Stochastic model for photoinduced anisotropy - Wiley Online Library
FULL PAPER
http://WWW.C-CHEM.ORG
Stochastic Model for Photoinduced Anisotropy Valentina Cantatore,[a] Giovanni Granucci,[a] and Maurizio Persico∗[a] We present a stochastic model for the kinetics of photoinduced anisotropy in a sample of molecular chromophores that may undergo photoisomerization. It is assumed that the chromophores do not interact among them, but are embedded in a medium that slows down the rotational diffusion. The model makes use of data about the photoinduced reorientation of the single chromophore, its photoisomerization and its rotational diffusion, that are made available by molecular dynamics
Introduction It is well-known that a sample of randomly oriented molecular chromophores can develop anisotropy when irradiated with linearly polarized light. This phenomenon is particularly important when the spontaneous thermal reorientation (“rotational diffusion”) of the molecules is slowed down because of the properties of the environment. Another condition that may greatly enhance the anisotropic response of the sample is the presence of mesogenic interactions between the chromophores. In such cases, a variety of technological applications have been devised.[1–8] In this article, we consider the relatively simple case of light absorbing molecules, that is, chromophores, that do not interact with each other but form a dilute solution in a medium that hinders the rotational diffusion (a viscous solvent, a low temperature glass or a polymeric matrix). As photoisomerization is often competing or cooperating with photo-orientation, we must take it into account. In a previous article, we presented a simple model for the photo-orientation of axial molecules, that translates into an explicit set of differential equations for the time evolution of the angular distribution.[9] Although that formal treatment is very convenient to show some general features of the problem, it does not consider the photoisomerization and the details of the nonadiabatic dynamics, such as the flexibility of the molecule and the excited state lifetimes. We note that other important phenomena related to molecular rotation in condensed phases can be interpreted on the basis of a small number of statistical parameters; for instance, the depolarization of the fluorescence of rigid molecules depends on the rotational diffusion coefficients, that can be assumed to be the same for the ground and the excited states, especially if the fluorophore is bound to a large molecule (protein or the like).[10] The photo-orientation process is instead due to major geometrical deformations of a molecule that has absorbed light and interacts with its chemical environment, so its description needs specific information concerning the molecular photodynamics. Nowadays, the nonadiabatic classical trajectory methods have made available a full description of the excited state http://onlinelibrary.wiley.com
simulations. For the first time such molecular scale processes are computationally connected to the development of anisotropy in a large sample and on a long time scale. A test on azobenzene shows the potentiality of the method and the interplay between photoinduced anisotropy and photoisomerization. © 2012 Wiley Periodicals, Inc. DOI: 10.1002/jcc.22931
molecular dynamics in condensed phase, for example, for the azo chromophore.[11–15] The simulation results contain the necessary information about the reorientation undergone by a single chromophore interacting with a specific chemical environment, as a consequence of electronic excitation, geometrical relaxation first in the excited and then in the ground state, and possibly also isomerization. The reorientation data take the form of distributions of rotation matrices (see next section), that depend on the specific compound and isomer, on the chemical environment, and sometimes also on the targeted absorption band, as in the case of azobenzene. The model described in this work makes use of this information for the first time, to describe the development of the anisotropy of a sample of many chromophores, including the case of two or more interconverting isomers. It is based on a stochastic description of the photochemical and thermal events that change the molecular orientation. It can be used to predict the mutual influence of photo-orientation and photoisomerization, and their dependence on the properties of the chromophores, of the medium and of the exciting light. It also allows to simulate the changes in light absorption that can be used experimentally to monitor the anisotropy and the isomeric ratio. An application to azobenzene is shortly outlined.
Model and Connection with Single Chromophore Molecular Dynamics The stochastic description of the photo-orientation phenomenon considers a sample of molecular chromophores not interacting among them, but embedded in a medium which is assumed to be isotropic at all times. The chromophores are normally in the ground state and undergo rotational diffusion. From time to time each of them may adsorb light and go through a [a] V. Cantatore, G. Granucci, M. Persico Dipartimento di Chimica e Chimica Industriale, Università di Pisa, v. Risorgimento 35, I-56126 Pisa, Italy E-mail : [email protected] Contract/grant sponsor: University of Pisa © 2012 Wiley Periodicals, Inc.
Journal of Computational Chemistry 2012, 33, 1015–1022
1015
FULL PAPER
http://WWW.C-CHEM.ORG
sequence of photochemical and photophysical processes that result in a photoinduced reorientation. As the probability of exciting a chromophore with linearly polarized light depends on the direction of the transition dipole vector and therefore on the molecular orientation, the process induces anisotropy in the sample. The orientation of each molecule is defined by the molecular frame (ˆx , yˆ , zˆ ), that can be defined in various ways, for instance as the principal axis system or with reference to some specific
cos α cos β cos γ − sin α sin γ R() = sin α cos β cos γ + cos α sin γ − sin β cos γ
(3)
In this way the orientation matrix is propagated in time under the influence of photochemical or thermal processes, as specified in the next paragraphs. We shall consider a large but finite number of molecules, NM , that are numbered with the index m. If the compound can be photoconverted to different isomers, at any time t we will have NI (t) molecules in the isomeric form I, with I NI = NM . The isomeric form of molecule m will be denoted by the index Im . Each molecule is also characterized by its orientation m , that corresponds to the rotation matrix Rm . Moreover, a molecule can evolve either photochemically or thermally, that is, at a given time it can be engaged in a photochemical or in a thermal trajectory. The former goes from the instant of photon absorption to a final time that is determined when running the nonadiabatic molecular dynamics simulation; typically the stop conditions require that the molecule is back in the ground electronic state, is near the equilibrium geometry of one of the isomers, and has lost a substantial fraction of the photon energy, so that no further reactive events are expected.[15, 17] When not undergoing a photochemical process, the molecule evolves according to a thermal trajectory, that only implies rotational diffusion. We assume that thermal isomerization is so slow as to be negligible in the time scale of our simulations, but this limitation can be easily dropped. The information defining the photochemical and thermal trajectories is obtained as a by-product of molecular dynamics simulations concerning one chromophore embedded in the chosen medium (solvent, glass or whatever else). The photochemistry is simulated by swarms of nonadiabatic trajectories, as can be obtained for instance by the surface hopping method.[11–15] A specific swarm of trajectories will represent the photochemical events following excitation of each isomer. 1016
Journal of Computational Chemistry 2012, 33, 1015–1022
ˆ Z) ˆ R() (ˆx , yˆ , zˆ ) = (Xˆ , Y,
(1)
With , we indicate the collection of the Euler angles α, β, and γ as defined for instance in Zare (Notice that in this book the Euler angles are called φ, θ and χ instead of α, β and γ ).[16] R() is the rotation matrix
− cos α cos β sin γ − sin α cos γ − sin α cos β sin γ + cos α cos γ sin β sin γ
If R((t)) ≡ R(t) is the orientation of a molecule at time t and R(tx ) is the further rotation that the molecule undergoes during the time interval tx , the final orientation is R(t + tx ) = R(t) R(tx )
atoms or groups. The relationship between the (ˆx , yˆ , zˆ ) axes and ˆ Z), ˆ is the laboratory ones, (Xˆ , Y,
cos α sin β sin α sin β cos β
(2)
A trajectory k in the swarm of isomer I starts from suitably sampled initial conditions (nuclear coordinates and momenta), with a vertical transition from the ground state to one of the excited states. The “internal time” tx of each photochemical trajectory is counted from the initial excitation and goes up to T (I,k) , the duration of the photochemical event (not necessarily the same for all trajectories). At any time tx in the range [0, T (I,k) ] we need to know: the orientation of the molecular chromophore, expressed as R(I,k) (tx ); the transition energies EL (tx ) − E0 (tx ) and dipole moments µ 0,L(I,k) (tx ) between the ground and the excited states; and the electronic state L(I,k) (tx ) where the trajectory is running, that can change in time because of nonadiabatic transitions. Notice that most trajectory simulations of excited state dynamics are nowadays run with direct methods, that is, by computing quantum chemically the electronic energies and wavefunctions at each time step: it is then straightforward to compute also the transition dipole moments. At time tx the reorientation with respect to the excitation time is expressed by the matrix R(tx ) = [R(I,k) (0)]−1 R(I,k) (tx )
(4)
that is used in eq. (3) to propagate in time the orientation of the molecule. The thermal rotational diffusion of each isomer I is represented by a time dependent rotation matrix R(I,rd) (trd ). The R(I,rd) (trd ) matrices might be obtained by theories of rotational diffusion that consider shape, charge distribution and other properties of the solute and of the medium.[18, 19] We prefer to extract R(I,rd) (trd ) from simulations of the ground state molecular dynamics, for consistency with the photochemical dynamics. In this way, we take full account of the coupling of the overall rotation with the internal degrees of freedom: this is important, because the flexibility offered by single bond torsions can greatly enhance the rate of rotational diffusion. It is computationally more efficient to run in parallel NS trajectories of length TS , each starting with different initial conditions, than a single one of equivalent length NS · TS . Two subtrajectories in the space R1(I,rd) (t) and R2(I,rd) (t) can be joined by taking http://WWW.CHEMISTRYVIEWS.COM
FULL PAPER
http://WWW.C-CHEM.ORG
R(I,rd) (trd ) = R1(I,rd) (trd )
−1 R1(I,rd) (TS ) R2(I,rd) (0) R2(I,rd) (trd − TS )
have for 0 < trd < TS for TS < trd < 2TS
(5)
In a similar way one can chain up all the NS subtrajectories in random order, and reuse them several times to generate an arbitrarily long rotational diffusion trajectory R(I,rd) (trd ). When a molecule must start a thermal trajectory, a time trd is chosen at random as the beginning of the trajectory. The reorientation undergone after a time tx spent in the thermal trajectory is then expressed by the matrix R(tx ) = [R(I,rd) (trd )]−1 R(I,rd) (trd + tx )
(6)
to be used in eq. (3) as in the case of a photochemical trajectory. When a molecule is in the ground state it can absorb a photon. We shall neglect multiphoton phenomena, including the absorption or stimulated emission from excited states. The rate of excitation to state L depends on the transition dipole vector µ 0,L . If µx , µy , and µz are the components of µ 0,L in the molecular frame, those in the laboratory frame are:
µX µx µY = R() µy µZ µz
Aexc,Z () =
1000 ln(10) εI (νexc ) F 3µ2Z ln(10) εI (νexc ) F 3µ2Z = 2 NA hνexc µ Eexc µ2 (8)
Here, εI (νexc ) is the extinction coefficient of isomer I in mol/L/cm, F is the irradiance in W/cm2 , NA is Avogadro’s number and Eexc is the excitation energy NA hνexc in kJ/mol. Note that equally simple formulas can be used for other polarizations, including the circular or elliptic ones. Moreover, the irradiance F can be a function of time, so one may easily simulate the effect of sequences of pulses of different polarizations (see our previous work with a simpler and more formal model for an example[9] ). During a time step a molecule can be excited only once, with the probability Pexc = Aexc,Z () texc ≤ Aexc,Z () t
Here, texc is the fraction of t during which the molecule is not already engaged in a photochemical trajectory (most of the times texc = t). In this way, the possibility that two photons are absorbed in the same time step is neglected, an approximation which is warranted only if the time step is short enough:
(7)
If the exciting light is linearly polarized along the Zˆ axis, the photon absorption rate Aexc () is proportional to µ2Z . In isotropic samples, the orientational average of µ2Z is µ2Z = µ2 /3, but for a single molecule the dependence of µ2Z on the Euler angles is contained in eqs. (2) and (7). It can be seen that µ2Z and the excitation rate only depend on the β and γ angles. Therefore, the angular distribution of molecules irradiated only with Zpolarized light will depend on β and γ , but not on α (i.e., the rotation angle around Zˆ itself ).
t
Eexc ln(10) εI (νexc ) F
According to our model, each molecule evolves in time independently, but we update the parameters describing the status of every molecule at regular intervals of length t, to compute the ensemble averages that describe anisotropy, isomerization, power adsorption, and other quantities. In this section, we first adopt the point of view of the single molecule to describe the key choices and approximations we have applied, and then we outline the workflow of the simulation method. In principle, the excitation rate can be computed on the basis of transition energies and dipoles determined by quantum chemistry calculations. However, a more realistic and accurate simulation of the excitation process is obtained by combining the experimental absorption intensity and bandshape data with the computed direction of the transition dipole vector. With ˆ we monochromatic light of frequency νexc , polarized along Z, http://onlinelibrary.wiley.com
(10)
In field-free conditions, that is, when simulating the decay of the anisotropy caused by the rotational diffusion, the time step can be arbitrarily long. For simplicity, we also disregard excitations that might take place during a photochemical trajectory, even in its final part that takes place in the ground state (however, the contribution of this part of the trajectory to the optical anisotropy is taken into account, see next section). Therefore, the irradiance F should not be too high, so that ln(10) εI (νexc ) F (I) T 1 Eexc
Implementation of the Time Evolution
(9)
(11)
Here, T (I) is the average duration of a photochemical trajectory starting with isomer I. At time t0 = 0, we distribute the molecules isotropically, using the algorithm described in Appendix to produce a set of cells of equal volume in the α, β, γ space, and placing a molecule in the center of each cell. Then, at each time step beginning at tn = nt(n = 0, 1, 2 . . .), we go through the following steps: A. Cycle over all molecules B. At the beginning of the step, if molecule m is engaged in a photochemical trajectory, no action is taken. Otherwise (molecule engaged in a thermal trajectory or tn = 0)… B1 Select at random from the swarm of isomer Im a photochemical trajectory k to be provisionally associated with the molecule m. Journal of Computational Chemistry 2012, 33, 1015–1022
1017
FULL PAPER
http://WWW.C-CHEM.ORG
B2 Compute the excitation probability Pexc by eqs. (8) and (9), using texc = t and the transition dipole µ associated with the initial excitation of trajectory k. B3 Generate a random number R ∈ [0, 1]. B4.1 If Pexc ≥ R then start the photochemical trajectory k. B4.2 If Pexc < R and tn > 0 then continue the thermal trajectory already under way. B4.3 If Pexc < R and tn = 0 then start a thermal trajectory, by choosing at random the starting time trd of eq. (6). C. If a photochemical trajectory k (started during the same time step or earlier) ends before the end of the time step… C1 Update the isomeric tag for molecule m. C2 Update the orientation of molecule m by applying eqs. (3) and (4). C3 Try a new excitation as in steps B1–B4, but setting texc equal to the residual time in the step; if Pexc < R, start a new thermal trajectory as in B4.3. D. At the end of the step… D1 If a photochemical trajectory is running, update the isomeric tag and the orientation as in C1 and C2. D2 If a thermal trajectory is running, update the orientation using eqs. (3) and (6). D3 Compute the contribution of molecule m to the densities and statistical averages (see next section).
Statistical Averages At each time step several collective properties are computed. The isomeric fractions are defined as
(β) ρI (β, t)
=
(γ )
(12)
2π
dα
ρI (γ , t) =
dγ ρI (, t)
(16)
0
0 2π
π
dα
sin β dβ ρI (, t)
(17)
0
0
All the density functions, that is, ρI and the reduced ones, can be computed and visualized as histograms (for ρI we resort to the partition method illustrated in Appendix). However, it is more convenient to expand them in truncated Fourier series, which is a way to approximate the discrete distributions as continuous functions. We shall use real Fourier expansions for (γ ) (β) ρI(α) and ρI , and a Fourier–Legendre expansion for ρI :
ρI(α) (α, t)
(α) Mmax
=
(β)
ρI (β, t) = (γ )
ρI (γ , t) =
AI,M (t) FM (α)
(α) M=−Mmax (β) Lmax
BI,L (t) PL (cosβ)
L=0 (γ ) Mmax
CI,M (t) FM (γ )
(18)
(19)
(20)
(γ ) M=−Mmax
The FM basis functions are normalized sines and cosines and the PL are normalized Legendre polynomials: −1/2 cos(Mα) for M < 0 π FM (α) = (2π)−1/2 for M = 0 π −1/2 sin(Mα) for M > 0 PL (cos β) =
NI fI = NM
2π
2L + 1 2
(21)
1/2 PL (cos β)
(22)
The Fourier coefficients are then: The anisotropy of the sample is fully expressed by the angular distributions ρI (, t), one for each isomer I. Each distribution is normalized to 1, to facilitate comparisons:
ρI (, t) d =
ρI (α, β, γ , t) sin β dα dβ dγ = 1
(13)
The actual density of isomer I is therefore NI ρI (, t). As the distribution is discrete, the function ρI (, t) is a combination of Dirac’s δ functions: ρI (, t) =
1 δI,Im (sin βm )−1 δ(α − αm ) δ(β − βm ) δ(γ − γm ) NI m (14)
To facilitate the visualization of the results, it is convenient to define three reduced densities: ρI(α) (α, t) =
π 0
1018
sinβdβ
2π
dγ ρI (, t)
0
Journal of Computational Chemistry 2012, 33, 1015–1022
(15)
AI,M (t) =
NM 1 δI,I FM (αm ) NI m=1 m
(23)
BI,L (t) =
NM 1 δI,I P (βm ) NI m=1 m L
(24)
CI,M (t) =
NM 1 δI,I FM (γm ) NI m=1 m
(25)
The coefficient of the constant function in each expansion is fixed, because of the normalization √ of the reduced density: AI,0 = CI,0 = (2π)−1/2 and BI,0 = 1/ 2. At the beginning of every time step we also evaluate the absorption rates Aprobe,X , Aprobe,Y , and Aprobe,Z for light of frequency νprobe and polarizations along the three axes, by averaging over the NM molecules of the sample. The formulas are identical to eq. (8), with the obvious replacements. For simplicity, we have assumed νprobe = νexc in the present version of the code. In computing the contribution of molecule m to the averages, if the molecule is evolving thermally we use in eq. (8) the transition dipole vector for the initial excitation of http://WWW.CHEMISTRYVIEWS.COM
FULL PAPER
http://WWW.C-CHEM.ORG
Figure 1. trans- and cis-azobenzene and definition of the molecular axes.
a randomly chosen photochemical trajectory, that is, the same that enters the calculation of the excitation probability for the time step. If the molecule m is evolving photochemically in an excited state, we neglect its contribution, because normally the necessary information about higher lying excited states is missing. The fraction of trajectories in excited states at any time must anyway be small (otherwise two-photon processes would not be negligible). Moreover, excited state absorption and stimulated emission partially cancel out. Once the photochemical trajectory has reverted to the ground state, the dipole vectors µ 0,L(I,k) for transitions to the excited states are available. If more than one excited state has been considered in the nonadiabatic dynamics, the choice of the 0 → L transition is dictated by the two criteria of minimum detuning EL − E0 − hνprobe and largest strength of the transition, which 2 . Assuming lorentzian band shapes is proportional to µ (I,k) 0,L with HWHM = hνprobe /5, we choose the state L for which the 2
/ (EL − E0 − hνprobe )2 + (hνprobe /5)2 is maxiquantity µ (I,k) 0,L mum. Once the probe absorption rates are known, the optical anisotropy is defined as Rtotal =
A − A⊥ A + 2A⊥
(26)
where we identify A with Aprobe,Z and A⊥ with (Aprobe,X + Aprobe,Y )/2. The averages Aprobe,X , Aprobe,Y , and Aprobe,Z are also computed separately for each isomer, whence the isomer specific anisotropies RI , that are related to Rtotal by Rtotal = fI RI (27)
reversible photochemical interconvertion. The computed quantum yields for n → π ∗ excitation in ethylene glycol were t→c = 0.369 and c→t = 0.486.[15] At λ = 436 nm, near the maximum of the n → π ∗ band,[20] the absorption coefficients in a protic solvent (methanol) are εtrans = 490 and εcis = 1140 mol/L/cm. The molecular axes were defined as shown in Figure 1, so that xˆ coincides with the N–N axis and zˆ is the C2 symmetry axis for both isomers. Four simulations were run with different values of the irradiance, up to F = 5 GW/cm2 . At t = 0, we assumed only trans-azobenzene to be present, as it is the more stable isomer. The number of molecules was 3,732,480, corresponding to angular grids with 180 mesh points for each angle (see Appendix). The dependence on time of the isomeric ratio fcis /ftrans and of the optical anisotropy R is shown in Figures 2 and 3. In plotting the simulation results, it is convenient to introduce an adimensional reduced time[9] τ : τ = Atrans t =
3 ln(10) εtrans (νexc ) F t Eexc
(28)
where Atrans is the absorption rate for a trans-azobenzene molecule when its transition dipole is parallel to the polarization axis. In practice, one reduced time unit is equivalent to 81 ps, with an irradiance of 1 GW/cm2 .
I
A Test Simulation: Azobenzene in Ethylene Glycol To test the model and the code, we ran simulations of the photoinduced anisotropy of azobenzene in a moderately viscous solvent, ethylene glycol (η = 16.1 mPa/s). The trajectory data were taken from our previous study of the photoisomerization of a single azobenzene chromophore in a cluster of about 470 solvent molecules, with excitation in the n → π ∗ band.[15] Azobenzene has two isomers, trans and cis, that undergo a http://onlinelibrary.wiley.com
Figure 2. Isomeric ratio fcis /ftrans as a function of time, for different values of the irradiance F.The horizontal lines show the asymptotic values.
Journal of Computational Chemistry 2012, 33, 1015–1022
1019
FULL PAPER
http://WWW.C-CHEM.ORG
Table 1. Asymptotic values of the isomeric ratio fcis /ftrans and of the optical anisotropy. fcis /ftrans
Irradiance F
Figure 3. Optical anisotropy ratio Rtotal as a function of time,for different values of the irradiance F. The horizontal lines show the asymptotic values. For F = 5 GW/cm2 we also plot Rtrans and Rcis .
At the beginning of each simulation, the trans isomer develops a negative anisotropy Rtrans because the selective excitation and photoisomerization of molecules with a large Z component of the transition dipole decreases A . On the contrary, at the beginning of the irradiation the small number of cis molecules freshly produced have a positive anisotropy: apparently the direction of the transition dipole of the trans molecules from which they derive is roughly conserved through the photoisomerization process. In general, the anisotropy RI of each isomer I is affected by two processes: the excitation of the isomer itself, that decreases RI through the photoinduced reorientation and isomerization, and the photoisomerization of the other isomer, that tends to increase RI . Which process prevails depends, among other factors, on the relative abundance of the two isomers. This is why the anisotropy is not necessarily a monotonic function of time, as shown most clearly by the simulation with F = 0.5 GW/cm2 . For long irradiation times (20–30 reduced time units) a stationary state is reached. Given the values of the absorption coefficients and of the quantum yields, in isotropic conditions the isomeric ratio would be
fcis ftrans
= iso
εtrans t→c = 0.3258 εcis c→t
(β)
Figure 4. Reduced densities ρI
1020
(29)
(γ )
and ρI
Journal of Computational Chemistry 2012, 33, 1015–1022
Optical anisotropy
(GW/cm2 )
Simulation
eq. (30)
Rtrans
Rcis
0.5 1.0 2.0 5.0
0.3096 0.2996 0.2931 0.2802
0.3114 0.3030 0.2968 0.2861
−0.0383 −0.0632 −0.0849 −0.1184
−0.0170 −0.0304 −0.0443 −0.0654
Rtotal −0.0292 −0.0498 −0.0683 −0.0971
This asymptotic isomeric ratio is not expected to change with the excitation wavelength λexc within the n → π ∗ absorption band, because the ratios of absorption coefficients and the quantum yields are approximately constant. Higher yields of the cis isomer are obtained by π → π ∗ excitation at λexc = 300–350 nm.[21] When the sample is anisotropic, the absorption rate of linearly polarized light by isomer I is altered by the factor (1 + 2RI )/3, so
fcis ftrans
= aniso
fcis ftrans
iso
1 + 2Rtrans 1 + 2Rcis
(30)
Table 1 shows the asymptotic values of the isomeric ratio fcis /ftrans and of the optical anisotropy. We see that eq. (30) justifies almost perfectly the deviation of the isomeric ratio from the value expected in isotropic conditions, which is a good test for the correct functioning of the stochastic algorithms. Qualitatively, we can say that the isomer which develops the larger asymptotic anisotropy (the trans one, in our case) “avoids to be excited” more effectively, and its concentration is therefore higher than in the isotropic case. The difference increases with the F value, so in practice with higher irradiances one gets a smaller isomerization yield, but of course the time needed to approach the final isomeric ratio decreases with F (in comparing the plots of Fig. 2 remember that the reduced time unit is inversely proportional to F). Equation (30) reduces to eq. (29) when extrapolating to low F values, that correspond to vanishing anisotropies. The small remaining discrepancy between the simulation results and the predictions of eq. (30), that is, the second and
, obtained at the stationary state with the irradiance F = 5 GW/cm2 .
http://WWW.CHEMISTRYVIEWS.COM
FULL PAPER
http://WWW.C-CHEM.ORG
third columns of Table 1, is due to the fact that the photoisomerization is not instantaneous. In deriving eq. (30), as in most kinetic models,[9, 22] the excited state lifetimes are neglected. Actually, during the nonadiabatic dynamics that follows the excitation, the further absorption (or emission) of a photon would be most likely ineffective in increasing the photoisomerization rate, and in our model we have chosen to ignore completely this possibility (see “Model and Connection with Single Chromophore Molecular Dynamics” section). Therefore, the small fraction of molecules that are undergoing a photochemical process cannot start a new one (this fraction is ∼2% with the highest irradiance, F = 5 GW/cm2 ). As at molecular level the trans → cis photoisomerization is in the average slower than the cis → trans one (about 1.9 ps versus 1.0 ps), the simulation yields a slightly smaller fcis /ftrans ratio than eq. (30). This discussion shows that certain simplifying assumptions of our model affect the results in a physically sound way and however to a very small extent. Figure 4 shows the reduced densities obtained with the irradiance F = 5 GW/cm2 at the stationary state (1000 ps, about 60 reduced time units). As already observed in “Model and Connection with Single Chromophore Molecular Dynamics” section there is no dependence on the α angle if only Z-polarized light is used. To understand how the density depends on β and γ , one needs to know the direction of the n → π ∗ transition dipole in the molecular frame. In trans-azobenzene, the n → π ∗ transition is symmetry forbidden at the equilibrium C2h geometry, but it borrows intensity from the lowest π → π ∗ state. The orientation of the transition dipole is then well defined: it lies in the xy plane, almost parallel to the long axis of the molecule, and it makes an angle of about 50◦ with the xˆ axis.[23, 24] At β = 0◦ and 180◦ , the xy plane is orthogonal to the light polarˆ therefore at these angles we find the maximum ization axis Z, (β) of the ρI density. For intermediate β angles, the orthogonality of the transition dipole with respect to Zˆ can be achieved by a rotation around the molecular zˆ axis, to γ 130◦ or 310◦ (γ ) (i.e., −50◦ ), where we find the maxima of ρI . The angular distribution of cis-azobenzene is less easily interpreted, because of two reasons: first, as cis is the less abundant isomer, it is strongly influenced by the photoisomerization of trans (this is (β) why the ρI density of cis is somewhat complementary to that of trans); second, the direction of its n → π ∗ transition dipole fluctuates because of the interconversion of two conformational enantiomers.[24–26] Although the investigation of the molecular mechanisms of azobenzene photo-orientation are out of the scope of this work, this short discussion of the angular distributions suffices to show that our model can yield a detailed description of the anisotropy, beyond the mere optical parameter.
Conclusions The model we have presented describes the kinetics of photoinduced anisotropy in a sample of molecular chromophores, taking into account their photoisomerization. The model makes use in a stochastic way of the statistical data that can be obtained by simulations of the dynamics of a single molecule, especially concerning the photoinduced reorientation (with or without photoisomerization) and the rotational diffusion. Such http://onlinelibrary.wiley.com
molecular scale processes are connected for the first time to the development of anisotropy in the bulk. The main limitation of the model is that it can only be applied to diluted solutions in viscous liquids or amorphous isotropic solids, because the chromophores are assumed not to interact among them. There are no limitations as to the time scale of the simulation, which is physically related to the rate of rotational diffusion. The effect of light with different intensities and polarizations, and of sequences of pulses, could be easily simulated. A test on azobenzene shows that even in such a simple system the photoinduced anisotropy, the photoisomerization, and the other molecular processes are coupled in a non trivial way. In particular, we find that the photoisomerization of molecules selected by the exciting light according to their orientation can affect in opposite directions the anisotropy of the two isomers, the reactant and the product. As a result, the time dependence of the anisotropy can be non monotonic.
Appendix Here, we describe an algorithm to partition the space spanned by the Euler angles α, β, and γ into a number of threedimensional (3D) cells of equal volume. Our algorithm is similar to one recently put forward by Leopardi,[27] but is simpler and limited to three dimensions, while Leopardi deals with hyperspheres in an arbitrary number of dimensions. The partition we propose is more flexible, in the sense that the number and therefore the width of the cells in each of the three coordinates can be arbitrarily chosen, while Leopardi minimizes the Euclidean diameter of each cell. The angle γ can be treated independently. Its range is [0, 2π ] and it can be partitioned into Nγ segments of length γ = 2π/Nγ , numbered with the index K = 0, . . . Nγ − 1. We shall place the center of the K -th segment at γ = K γ . The problem is thus reduced to the definition of 2D cells in the {α, β} space, that is, the tessellation of a spherical surface. α and β are the azimuthal and zenithal angles, more often called φ and θ, respectively. Each 2D cell will have a trapezoidal shape and its boundaries will be segments of parallels and meridians. An exception are the two polar cells, each one being a spherical cap with center at β = 0 or π. Besides the polar caps, we have Nβ belts, each enclosed between two parallels and containing a variable number of cells. The J-th belt has a latitudinal width βJ and its NJ cells have a longitudinal width αJ = 2π/NJ . The dependence of the βJ width on the belt index J is minimized by our algorithm, but cannot be totally eliminated because of the constraint of equal surface for all the 2D cells. We shall adopt the following definitions: A. Nβ ≡ number of latitudinal belts, not including the polar caps. The latter are numbered J = 0 and J = Nβ + 1. The partitions of the north and south hemispheres are equivalent. If Nβ is odd, the equator is the midline of the central belt, otherwise the equator is the boundary between two belts. B. Nα ≡ number of cells in the equatorial belts (those having the equator as their midline or as a boundary). Journal of Computational Chemistry 2012, 33, 1015–1022
1021
FULL PAPER
http://WWW.C-CHEM.ORG
C. NC ≡ total number of 2D cells in the α, β space.
Acknowledgments
D. σ = 4π/NC ≡ surface of a 2D cell.
The authors are grateful to Piet Van Leuven and Christian Silvio Pomelli for helpful discussions.
E. βP = boundary of the north polar cap (symmetrically, the boundary of the southern cap is π − βP ). F. σcap (βP ) = 2π [1 − cos(βP )] = surface of the polar caps. G. βN,J ≡ northern (lower β) boundary of the J-th belt. H. βS,J ≡ southern (upper β) boundary of the J-th belt. I. σbelt,J = 2π [cos(βS,J ) − cos(βN,J )] ≡ surface of the J-th belt cos(β )+cos(β ) N,J S,J J. βC,J = arccos ≡ baricentric β of the J-th belt 2 (in this way we have the same surface in the two halves of the belt, above and below βC,J ). For the north cap βC,0 = 0 and for the south one βC,Nβ +1 = π . Nα and Nβ , as well as Nγ , are input data of the algorithm. In many problems where latitude is more important than longitude, one may choose Nβ > Nα . The algorithm first computes an approximate partition of the sphere in belts, the number of cells NJ for each belt, the total number of 2D cells NC and the area of the cells σ . Then, it proceeds to define the exact boundaries of the belts that contain NJ cells of the wanted surface. In detail, we go through the following steps: 1. The width βS,J −βN,J of each belt is approximated as π/(Nβ +1), and the radius of the polar caps as half of it: βP π/(2Nβ + 2).
Keywords: photo-orientation • stochastic chemistry models • photoinduced anisotropy • photoisomerization • azobenzene How to cite this article: V. Cantatore, G. Granucci, M. Persico, J. Comput. Chem. 2012, 33, 1015–1022. DOI: 10.1002/jcc.22931 [1] C. Cojocariu, P. Rochon, Pure. Appl. Chem 2004, 76, 1479. [2] K. G. Yager, C. J. Barrett, J. Photochem. Photobiol. A 2006, 182, 250. [3] I. Vecchi, A. Arcioni, C. Bacchiocchi, G. Tiberio, P. Zanirato, C. Zannoni, J. Phys. Chem. B 2007, 111, 3355. [4] D. Rais, Y. Zakrevskyy, J. Stumpe, S. Nešpøurek, Z. Sedláková, Opt. Mater 2008, 30, 1335. [5] K. Oki, Y. Nagasaka, Colloids. Surf. A Physicochem. Eng. Asp. 2009, 333, 182. ˇ [6] V. Domenici, G. Ambrožic, M. Copic, A. Lebar, I. Drevenšek-Olenik, P. Umek, B. Zalar, B. Zupanˇciˇc, M. Žigon, Polymer 2009, 50, 4837. [7] M. Devetak, B. Zupanˇciˇc, A. Lebar, P. Umek, B. Zalar, V. Domenici, G. ˇ Ambrožic, M. Žigon, M. Copic, I. Drevenšek-Olenik, Phys. Rev. E 2009, 80, 050701. [8] T. Yoshino, M. Kondo, J. Mamiya, M. Kinoshita, Y. Yu, T. Ikeda, Adv. Mater 2010, 22, 1361. [9] P. Van Leuven, C. Cantatore, M. Persico, Phys. Chem. Chem. Phys, in press, DOI: 10.1039/c2cp23128j. [10] J. R. Lakowicz, Principles of Fluorescence Spectroscopy; Springer: New York, 2006. [11] P. Cattaneo, M. Persico, J. Am. Chem. Soc. 2001, 123, 7638. [12] M. Böckmann, N. L. Doltsinis, D. Marx, Phys. Rev. E 2008, 78, 036101. [13] M. Böckmann, N. L. Doltsinis, D. Marx, J. Phys. Chem. A 2010, 114, 745.
2. The approximate central β of the J-th belt is then set to βC,J = πJ/(Nβ + 1).
[14] G. Tiberio, L. Muccioli, R. Berardi, C. Zannoni, ChemPhysChem 2010, 11, 1018. [15] T. Cusati, G. Granucci, M. Persico, J. Am. Chem. Soc 2011, 133, 5109.
3. The number of cells contained in the J-th belt is definitively set to NJ = Nα sin(βC,J ), by rounding to the nearest integer, because the surface σbelt,J is approximately proportional to sin(βC,J ). For the polar caps N0 = NNβ +1 = 1. 4. The total number of 2D cells is then NC =
1022
Nβ +1 J=0
NJ .
[16] R. N. Zare, Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics, J. Wiley (1988). [17] L. Creatini, T. Cusati, G. Granucci, M. Persico, Chem. Phys. 2008, 347, 492. [18] B. J. Berne, J. Chem. Phys. 1975, 62, 1154. [19] G. Lévi, J. P. Marsault, F. Marsault-Hérail, R. E. D. McClung, J. Chem. Phys. 1975, 63, 3543. [20] G. Gauglitz, S. Hubig, J. Photochem. 1985, 30, 121.
5. The surface of a 2D cell is computed as σ = 4π/NC .
[21] H. Rau, In Photochromism: Molecules and Systems; H. Dürr, H. BouasLaurent Eds. 2003.
6. The radius of the polar caps, that is, the southern boundary of the north cap, is computed as βP = βS,0 = arccos(1 − σ/2π).
[22] M. Persico, P. Van Leuven, Phys. Chem. Chem. Phys. 2009, 11, 8433.
7. Starting from J = 1 and going incrementally to J = Nβ : - we put βN,J = βS,J−1 ; - we compute σbelt,J = NJ σ ; - we compute βS,J = arccos[cos(βN,J ) − σbelt,J /2π ], as follows from definition (I); - we compute βC,J from definition (J).
[24] T. Cusati, G. Granucci, E. Martínez-Núñez, F. Martini, M. Persico, S. Vázquez, J. Phys. Chem. A, in press, DOI: 10.1021/jp208574q.
A 3D cell is identified with the indexes I, J, and K of the angles α, β, and γ . The center of the cell is put at αI = 2πI/NJ , βJ = βC,J , and γK = K γ . The volume of a 3D cell is V = σ γ = 8π 2 /(NC Nγ ).
Received: 17 October 2011 Revised: 27 November 2011 Accepted: 16 December 2011 Published online on 13 February 2012
Journal of Computational Chemistry 2012, 33, 1015–1022
[23] T. Cusati, G. Granucci, M. Persico, G. Spighi, J. Chem. Phys. 2008, 128, 194312.
[25] Y. Ootani, K. Satoh, A. Nakayama, T. Noro, T. Taketsugu, J. Chem. Phys. 2009, 131, 194306. [26] O. Weingart, Z. Lan, A. Koslowski, W. Thiel, J. Phys. Chem. Lett. 2011, 2, 1506. [27] P. Leopardi, Electron Trans. Numer. Anal. 2006, 25, 309.
http://WWW.CHEMISTRYVIEWS.COM
http://WWW.C-CHEM.ORG
Stochastic Model for Photoinduced Anisotropy Valentina Cantatore,[a] Giovanni Granucci,[a] and Maurizio Persico∗[a] We present a stochastic model for the kinetics of photoinduced anisotropy in a sample of molecular chromophores that may undergo photoisomerization. It is assumed that the chromophores do not interact among them, but are embedded in a medium that slows down the rotational diffusion. The model makes use of data about the photoinduced reorientation of the single chromophore, its photoisomerization and its rotational diffusion, that are made available by molecular dynamics
Introduction It is well-known that a sample of randomly oriented molecular chromophores can develop anisotropy when irradiated with linearly polarized light. This phenomenon is particularly important when the spontaneous thermal reorientation (“rotational diffusion”) of the molecules is slowed down because of the properties of the environment. Another condition that may greatly enhance the anisotropic response of the sample is the presence of mesogenic interactions between the chromophores. In such cases, a variety of technological applications have been devised.[1–8] In this article, we consider the relatively simple case of light absorbing molecules, that is, chromophores, that do not interact with each other but form a dilute solution in a medium that hinders the rotational diffusion (a viscous solvent, a low temperature glass or a polymeric matrix). As photoisomerization is often competing or cooperating with photo-orientation, we must take it into account. In a previous article, we presented a simple model for the photo-orientation of axial molecules, that translates into an explicit set of differential equations for the time evolution of the angular distribution.[9] Although that formal treatment is very convenient to show some general features of the problem, it does not consider the photoisomerization and the details of the nonadiabatic dynamics, such as the flexibility of the molecule and the excited state lifetimes. We note that other important phenomena related to molecular rotation in condensed phases can be interpreted on the basis of a small number of statistical parameters; for instance, the depolarization of the fluorescence of rigid molecules depends on the rotational diffusion coefficients, that can be assumed to be the same for the ground and the excited states, especially if the fluorophore is bound to a large molecule (protein or the like).[10] The photo-orientation process is instead due to major geometrical deformations of a molecule that has absorbed light and interacts with its chemical environment, so its description needs specific information concerning the molecular photodynamics. Nowadays, the nonadiabatic classical trajectory methods have made available a full description of the excited state http://onlinelibrary.wiley.com
simulations. For the first time such molecular scale processes are computationally connected to the development of anisotropy in a large sample and on a long time scale. A test on azobenzene shows the potentiality of the method and the interplay between photoinduced anisotropy and photoisomerization. © 2012 Wiley Periodicals, Inc. DOI: 10.1002/jcc.22931
molecular dynamics in condensed phase, for example, for the azo chromophore.[11–15] The simulation results contain the necessary information about the reorientation undergone by a single chromophore interacting with a specific chemical environment, as a consequence of electronic excitation, geometrical relaxation first in the excited and then in the ground state, and possibly also isomerization. The reorientation data take the form of distributions of rotation matrices (see next section), that depend on the specific compound and isomer, on the chemical environment, and sometimes also on the targeted absorption band, as in the case of azobenzene. The model described in this work makes use of this information for the first time, to describe the development of the anisotropy of a sample of many chromophores, including the case of two or more interconverting isomers. It is based on a stochastic description of the photochemical and thermal events that change the molecular orientation. It can be used to predict the mutual influence of photo-orientation and photoisomerization, and their dependence on the properties of the chromophores, of the medium and of the exciting light. It also allows to simulate the changes in light absorption that can be used experimentally to monitor the anisotropy and the isomeric ratio. An application to azobenzene is shortly outlined.
Model and Connection with Single Chromophore Molecular Dynamics The stochastic description of the photo-orientation phenomenon considers a sample of molecular chromophores not interacting among them, but embedded in a medium which is assumed to be isotropic at all times. The chromophores are normally in the ground state and undergo rotational diffusion. From time to time each of them may adsorb light and go through a [a] V. Cantatore, G. Granucci, M. Persico Dipartimento di Chimica e Chimica Industriale, Università di Pisa, v. Risorgimento 35, I-56126 Pisa, Italy E-mail : [email protected] Contract/grant sponsor: University of Pisa © 2012 Wiley Periodicals, Inc.
Journal of Computational Chemistry 2012, 33, 1015–1022
1015
FULL PAPER
http://WWW.C-CHEM.ORG
sequence of photochemical and photophysical processes that result in a photoinduced reorientation. As the probability of exciting a chromophore with linearly polarized light depends on the direction of the transition dipole vector and therefore on the molecular orientation, the process induces anisotropy in the sample. The orientation of each molecule is defined by the molecular frame (ˆx , yˆ , zˆ ), that can be defined in various ways, for instance as the principal axis system or with reference to some specific
cos α cos β cos γ − sin α sin γ R() = sin α cos β cos γ + cos α sin γ − sin β cos γ
(3)
In this way the orientation matrix is propagated in time under the influence of photochemical or thermal processes, as specified in the next paragraphs. We shall consider a large but finite number of molecules, NM , that are numbered with the index m. If the compound can be photoconverted to different isomers, at any time t we will have NI (t) molecules in the isomeric form I, with I NI = NM . The isomeric form of molecule m will be denoted by the index Im . Each molecule is also characterized by its orientation m , that corresponds to the rotation matrix Rm . Moreover, a molecule can evolve either photochemically or thermally, that is, at a given time it can be engaged in a photochemical or in a thermal trajectory. The former goes from the instant of photon absorption to a final time that is determined when running the nonadiabatic molecular dynamics simulation; typically the stop conditions require that the molecule is back in the ground electronic state, is near the equilibrium geometry of one of the isomers, and has lost a substantial fraction of the photon energy, so that no further reactive events are expected.[15, 17] When not undergoing a photochemical process, the molecule evolves according to a thermal trajectory, that only implies rotational diffusion. We assume that thermal isomerization is so slow as to be negligible in the time scale of our simulations, but this limitation can be easily dropped. The information defining the photochemical and thermal trajectories is obtained as a by-product of molecular dynamics simulations concerning one chromophore embedded in the chosen medium (solvent, glass or whatever else). The photochemistry is simulated by swarms of nonadiabatic trajectories, as can be obtained for instance by the surface hopping method.[11–15] A specific swarm of trajectories will represent the photochemical events following excitation of each isomer. 1016
Journal of Computational Chemistry 2012, 33, 1015–1022
ˆ Z) ˆ R() (ˆx , yˆ , zˆ ) = (Xˆ , Y,
(1)
With , we indicate the collection of the Euler angles α, β, and γ as defined for instance in Zare (Notice that in this book the Euler angles are called φ, θ and χ instead of α, β and γ ).[16] R() is the rotation matrix
− cos α cos β sin γ − sin α cos γ − sin α cos β sin γ + cos α cos γ sin β sin γ
If R((t)) ≡ R(t) is the orientation of a molecule at time t and R(tx ) is the further rotation that the molecule undergoes during the time interval tx , the final orientation is R(t + tx ) = R(t) R(tx )
atoms or groups. The relationship between the (ˆx , yˆ , zˆ ) axes and ˆ Z), ˆ is the laboratory ones, (Xˆ , Y,
cos α sin β sin α sin β cos β
(2)
A trajectory k in the swarm of isomer I starts from suitably sampled initial conditions (nuclear coordinates and momenta), with a vertical transition from the ground state to one of the excited states. The “internal time” tx of each photochemical trajectory is counted from the initial excitation and goes up to T (I,k) , the duration of the photochemical event (not necessarily the same for all trajectories). At any time tx in the range [0, T (I,k) ] we need to know: the orientation of the molecular chromophore, expressed as R(I,k) (tx ); the transition energies EL (tx ) − E0 (tx ) and dipole moments µ 0,L(I,k) (tx ) between the ground and the excited states; and the electronic state L(I,k) (tx ) where the trajectory is running, that can change in time because of nonadiabatic transitions. Notice that most trajectory simulations of excited state dynamics are nowadays run with direct methods, that is, by computing quantum chemically the electronic energies and wavefunctions at each time step: it is then straightforward to compute also the transition dipole moments. At time tx the reorientation with respect to the excitation time is expressed by the matrix R(tx ) = [R(I,k) (0)]−1 R(I,k) (tx )
(4)
that is used in eq. (3) to propagate in time the orientation of the molecule. The thermal rotational diffusion of each isomer I is represented by a time dependent rotation matrix R(I,rd) (trd ). The R(I,rd) (trd ) matrices might be obtained by theories of rotational diffusion that consider shape, charge distribution and other properties of the solute and of the medium.[18, 19] We prefer to extract R(I,rd) (trd ) from simulations of the ground state molecular dynamics, for consistency with the photochemical dynamics. In this way, we take full account of the coupling of the overall rotation with the internal degrees of freedom: this is important, because the flexibility offered by single bond torsions can greatly enhance the rate of rotational diffusion. It is computationally more efficient to run in parallel NS trajectories of length TS , each starting with different initial conditions, than a single one of equivalent length NS · TS . Two subtrajectories in the space R1(I,rd) (t) and R2(I,rd) (t) can be joined by taking http://WWW.CHEMISTRYVIEWS.COM
FULL PAPER
http://WWW.C-CHEM.ORG
R(I,rd) (trd ) = R1(I,rd) (trd )
−1 R1(I,rd) (TS ) R2(I,rd) (0) R2(I,rd) (trd − TS )
have for 0 < trd < TS for TS < trd < 2TS
(5)
In a similar way one can chain up all the NS subtrajectories in random order, and reuse them several times to generate an arbitrarily long rotational diffusion trajectory R(I,rd) (trd ). When a molecule must start a thermal trajectory, a time trd is chosen at random as the beginning of the trajectory. The reorientation undergone after a time tx spent in the thermal trajectory is then expressed by the matrix R(tx ) = [R(I,rd) (trd )]−1 R(I,rd) (trd + tx )
(6)
to be used in eq. (3) as in the case of a photochemical trajectory. When a molecule is in the ground state it can absorb a photon. We shall neglect multiphoton phenomena, including the absorption or stimulated emission from excited states. The rate of excitation to state L depends on the transition dipole vector µ 0,L . If µx , µy , and µz are the components of µ 0,L in the molecular frame, those in the laboratory frame are:
µX µx µY = R() µy µZ µz
Aexc,Z () =
1000 ln(10) εI (νexc ) F 3µ2Z ln(10) εI (νexc ) F 3µ2Z = 2 NA hνexc µ Eexc µ2 (8)
Here, εI (νexc ) is the extinction coefficient of isomer I in mol/L/cm, F is the irradiance in W/cm2 , NA is Avogadro’s number and Eexc is the excitation energy NA hνexc in kJ/mol. Note that equally simple formulas can be used for other polarizations, including the circular or elliptic ones. Moreover, the irradiance F can be a function of time, so one may easily simulate the effect of sequences of pulses of different polarizations (see our previous work with a simpler and more formal model for an example[9] ). During a time step a molecule can be excited only once, with the probability Pexc = Aexc,Z () texc ≤ Aexc,Z () t
Here, texc is the fraction of t during which the molecule is not already engaged in a photochemical trajectory (most of the times texc = t). In this way, the possibility that two photons are absorbed in the same time step is neglected, an approximation which is warranted only if the time step is short enough:
(7)
If the exciting light is linearly polarized along the Zˆ axis, the photon absorption rate Aexc () is proportional to µ2Z . In isotropic samples, the orientational average of µ2Z is µ2Z = µ2 /3, but for a single molecule the dependence of µ2Z on the Euler angles is contained in eqs. (2) and (7). It can be seen that µ2Z and the excitation rate only depend on the β and γ angles. Therefore, the angular distribution of molecules irradiated only with Zpolarized light will depend on β and γ , but not on α (i.e., the rotation angle around Zˆ itself ).
t
Eexc ln(10) εI (νexc ) F
According to our model, each molecule evolves in time independently, but we update the parameters describing the status of every molecule at regular intervals of length t, to compute the ensemble averages that describe anisotropy, isomerization, power adsorption, and other quantities. In this section, we first adopt the point of view of the single molecule to describe the key choices and approximations we have applied, and then we outline the workflow of the simulation method. In principle, the excitation rate can be computed on the basis of transition energies and dipoles determined by quantum chemistry calculations. However, a more realistic and accurate simulation of the excitation process is obtained by combining the experimental absorption intensity and bandshape data with the computed direction of the transition dipole vector. With ˆ we monochromatic light of frequency νexc , polarized along Z, http://onlinelibrary.wiley.com
(10)
In field-free conditions, that is, when simulating the decay of the anisotropy caused by the rotational diffusion, the time step can be arbitrarily long. For simplicity, we also disregard excitations that might take place during a photochemical trajectory, even in its final part that takes place in the ground state (however, the contribution of this part of the trajectory to the optical anisotropy is taken into account, see next section). Therefore, the irradiance F should not be too high, so that ln(10) εI (νexc ) F (I) T 1 Eexc
Implementation of the Time Evolution
(9)
(11)
Here, T (I) is the average duration of a photochemical trajectory starting with isomer I. At time t0 = 0, we distribute the molecules isotropically, using the algorithm described in Appendix to produce a set of cells of equal volume in the α, β, γ space, and placing a molecule in the center of each cell. Then, at each time step beginning at tn = nt(n = 0, 1, 2 . . .), we go through the following steps: A. Cycle over all molecules B. At the beginning of the step, if molecule m is engaged in a photochemical trajectory, no action is taken. Otherwise (molecule engaged in a thermal trajectory or tn = 0)… B1 Select at random from the swarm of isomer Im a photochemical trajectory k to be provisionally associated with the molecule m. Journal of Computational Chemistry 2012, 33, 1015–1022
1017
FULL PAPER
http://WWW.C-CHEM.ORG
B2 Compute the excitation probability Pexc by eqs. (8) and (9), using texc = t and the transition dipole µ associated with the initial excitation of trajectory k. B3 Generate a random number R ∈ [0, 1]. B4.1 If Pexc ≥ R then start the photochemical trajectory k. B4.2 If Pexc < R and tn > 0 then continue the thermal trajectory already under way. B4.3 If Pexc < R and tn = 0 then start a thermal trajectory, by choosing at random the starting time trd of eq. (6). C. If a photochemical trajectory k (started during the same time step or earlier) ends before the end of the time step… C1 Update the isomeric tag for molecule m. C2 Update the orientation of molecule m by applying eqs. (3) and (4). C3 Try a new excitation as in steps B1–B4, but setting texc equal to the residual time in the step; if Pexc < R, start a new thermal trajectory as in B4.3. D. At the end of the step… D1 If a photochemical trajectory is running, update the isomeric tag and the orientation as in C1 and C2. D2 If a thermal trajectory is running, update the orientation using eqs. (3) and (6). D3 Compute the contribution of molecule m to the densities and statistical averages (see next section).
Statistical Averages At each time step several collective properties are computed. The isomeric fractions are defined as
(β) ρI (β, t)
=
(γ )
(12)
2π
dα
ρI (γ , t) =
dγ ρI (, t)
(16)
0
0 2π
π
dα
sin β dβ ρI (, t)
(17)
0
0
All the density functions, that is, ρI and the reduced ones, can be computed and visualized as histograms (for ρI we resort to the partition method illustrated in Appendix). However, it is more convenient to expand them in truncated Fourier series, which is a way to approximate the discrete distributions as continuous functions. We shall use real Fourier expansions for (γ ) (β) ρI(α) and ρI , and a Fourier–Legendre expansion for ρI :
ρI(α) (α, t)
(α) Mmax
=
(β)
ρI (β, t) = (γ )
ρI (γ , t) =
AI,M (t) FM (α)
(α) M=−Mmax (β) Lmax
BI,L (t) PL (cosβ)
L=0 (γ ) Mmax
CI,M (t) FM (γ )
(18)
(19)
(20)
(γ ) M=−Mmax
The FM basis functions are normalized sines and cosines and the PL are normalized Legendre polynomials: −1/2 cos(Mα) for M < 0 π FM (α) = (2π)−1/2 for M = 0 π −1/2 sin(Mα) for M > 0 PL (cos β) =
NI fI = NM
2π
2L + 1 2
(21)
1/2 PL (cos β)
(22)
The Fourier coefficients are then: The anisotropy of the sample is fully expressed by the angular distributions ρI (, t), one for each isomer I. Each distribution is normalized to 1, to facilitate comparisons:
ρI (, t) d =
ρI (α, β, γ , t) sin β dα dβ dγ = 1
(13)
The actual density of isomer I is therefore NI ρI (, t). As the distribution is discrete, the function ρI (, t) is a combination of Dirac’s δ functions: ρI (, t) =
1 δI,Im (sin βm )−1 δ(α − αm ) δ(β − βm ) δ(γ − γm ) NI m (14)
To facilitate the visualization of the results, it is convenient to define three reduced densities: ρI(α) (α, t) =
π 0
1018
sinβdβ
2π
dγ ρI (, t)
0
Journal of Computational Chemistry 2012, 33, 1015–1022
(15)
AI,M (t) =
NM 1 δI,I FM (αm ) NI m=1 m
(23)
BI,L (t) =
NM 1 δI,I P (βm ) NI m=1 m L
(24)
CI,M (t) =
NM 1 δI,I FM (γm ) NI m=1 m
(25)
The coefficient of the constant function in each expansion is fixed, because of the normalization √ of the reduced density: AI,0 = CI,0 = (2π)−1/2 and BI,0 = 1/ 2. At the beginning of every time step we also evaluate the absorption rates Aprobe,X , Aprobe,Y , and Aprobe,Z for light of frequency νprobe and polarizations along the three axes, by averaging over the NM molecules of the sample. The formulas are identical to eq. (8), with the obvious replacements. For simplicity, we have assumed νprobe = νexc in the present version of the code. In computing the contribution of molecule m to the averages, if the molecule is evolving thermally we use in eq. (8) the transition dipole vector for the initial excitation of http://WWW.CHEMISTRYVIEWS.COM
FULL PAPER
http://WWW.C-CHEM.ORG
Figure 1. trans- and cis-azobenzene and definition of the molecular axes.
a randomly chosen photochemical trajectory, that is, the same that enters the calculation of the excitation probability for the time step. If the molecule m is evolving photochemically in an excited state, we neglect its contribution, because normally the necessary information about higher lying excited states is missing. The fraction of trajectories in excited states at any time must anyway be small (otherwise two-photon processes would not be negligible). Moreover, excited state absorption and stimulated emission partially cancel out. Once the photochemical trajectory has reverted to the ground state, the dipole vectors µ 0,L(I,k) for transitions to the excited states are available. If more than one excited state has been considered in the nonadiabatic dynamics, the choice of the 0 → L transition is dictated by the two criteria of minimum detuning EL − E0 − hνprobe and largest strength of the transition, which 2 . Assuming lorentzian band shapes is proportional to µ (I,k) 0,L with HWHM = hνprobe /5, we choose the state L for which the 2
/ (EL − E0 − hνprobe )2 + (hνprobe /5)2 is maxiquantity µ (I,k) 0,L mum. Once the probe absorption rates are known, the optical anisotropy is defined as Rtotal =
A − A⊥ A + 2A⊥
(26)
where we identify A with Aprobe,Z and A⊥ with (Aprobe,X + Aprobe,Y )/2. The averages Aprobe,X , Aprobe,Y , and Aprobe,Z are also computed separately for each isomer, whence the isomer specific anisotropies RI , that are related to Rtotal by Rtotal = fI RI (27)
reversible photochemical interconvertion. The computed quantum yields for n → π ∗ excitation in ethylene glycol were t→c = 0.369 and c→t = 0.486.[15] At λ = 436 nm, near the maximum of the n → π ∗ band,[20] the absorption coefficients in a protic solvent (methanol) are εtrans = 490 and εcis = 1140 mol/L/cm. The molecular axes were defined as shown in Figure 1, so that xˆ coincides with the N–N axis and zˆ is the C2 symmetry axis for both isomers. Four simulations were run with different values of the irradiance, up to F = 5 GW/cm2 . At t = 0, we assumed only trans-azobenzene to be present, as it is the more stable isomer. The number of molecules was 3,732,480, corresponding to angular grids with 180 mesh points for each angle (see Appendix). The dependence on time of the isomeric ratio fcis /ftrans and of the optical anisotropy R is shown in Figures 2 and 3. In plotting the simulation results, it is convenient to introduce an adimensional reduced time[9] τ : τ = Atrans t =
3 ln(10) εtrans (νexc ) F t Eexc
(28)
where Atrans is the absorption rate for a trans-azobenzene molecule when its transition dipole is parallel to the polarization axis. In practice, one reduced time unit is equivalent to 81 ps, with an irradiance of 1 GW/cm2 .
I
A Test Simulation: Azobenzene in Ethylene Glycol To test the model and the code, we ran simulations of the photoinduced anisotropy of azobenzene in a moderately viscous solvent, ethylene glycol (η = 16.1 mPa/s). The trajectory data were taken from our previous study of the photoisomerization of a single azobenzene chromophore in a cluster of about 470 solvent molecules, with excitation in the n → π ∗ band.[15] Azobenzene has two isomers, trans and cis, that undergo a http://onlinelibrary.wiley.com
Figure 2. Isomeric ratio fcis /ftrans as a function of time, for different values of the irradiance F.The horizontal lines show the asymptotic values.
Journal of Computational Chemistry 2012, 33, 1015–1022
1019
FULL PAPER
http://WWW.C-CHEM.ORG
Table 1. Asymptotic values of the isomeric ratio fcis /ftrans and of the optical anisotropy. fcis /ftrans
Irradiance F
Figure 3. Optical anisotropy ratio Rtotal as a function of time,for different values of the irradiance F. The horizontal lines show the asymptotic values. For F = 5 GW/cm2 we also plot Rtrans and Rcis .
At the beginning of each simulation, the trans isomer develops a negative anisotropy Rtrans because the selective excitation and photoisomerization of molecules with a large Z component of the transition dipole decreases A . On the contrary, at the beginning of the irradiation the small number of cis molecules freshly produced have a positive anisotropy: apparently the direction of the transition dipole of the trans molecules from which they derive is roughly conserved through the photoisomerization process. In general, the anisotropy RI of each isomer I is affected by two processes: the excitation of the isomer itself, that decreases RI through the photoinduced reorientation and isomerization, and the photoisomerization of the other isomer, that tends to increase RI . Which process prevails depends, among other factors, on the relative abundance of the two isomers. This is why the anisotropy is not necessarily a monotonic function of time, as shown most clearly by the simulation with F = 0.5 GW/cm2 . For long irradiation times (20–30 reduced time units) a stationary state is reached. Given the values of the absorption coefficients and of the quantum yields, in isotropic conditions the isomeric ratio would be
fcis ftrans
= iso
εtrans t→c = 0.3258 εcis c→t
(β)
Figure 4. Reduced densities ρI
1020
(29)
(γ )
and ρI
Journal of Computational Chemistry 2012, 33, 1015–1022
Optical anisotropy
(GW/cm2 )
Simulation
eq. (30)
Rtrans
Rcis
0.5 1.0 2.0 5.0
0.3096 0.2996 0.2931 0.2802
0.3114 0.3030 0.2968 0.2861
−0.0383 −0.0632 −0.0849 −0.1184
−0.0170 −0.0304 −0.0443 −0.0654
Rtotal −0.0292 −0.0498 −0.0683 −0.0971
This asymptotic isomeric ratio is not expected to change with the excitation wavelength λexc within the n → π ∗ absorption band, because the ratios of absorption coefficients and the quantum yields are approximately constant. Higher yields of the cis isomer are obtained by π → π ∗ excitation at λexc = 300–350 nm.[21] When the sample is anisotropic, the absorption rate of linearly polarized light by isomer I is altered by the factor (1 + 2RI )/3, so
fcis ftrans
= aniso
fcis ftrans
iso
1 + 2Rtrans 1 + 2Rcis
(30)
Table 1 shows the asymptotic values of the isomeric ratio fcis /ftrans and of the optical anisotropy. We see that eq. (30) justifies almost perfectly the deviation of the isomeric ratio from the value expected in isotropic conditions, which is a good test for the correct functioning of the stochastic algorithms. Qualitatively, we can say that the isomer which develops the larger asymptotic anisotropy (the trans one, in our case) “avoids to be excited” more effectively, and its concentration is therefore higher than in the isotropic case. The difference increases with the F value, so in practice with higher irradiances one gets a smaller isomerization yield, but of course the time needed to approach the final isomeric ratio decreases with F (in comparing the plots of Fig. 2 remember that the reduced time unit is inversely proportional to F). Equation (30) reduces to eq. (29) when extrapolating to low F values, that correspond to vanishing anisotropies. The small remaining discrepancy between the simulation results and the predictions of eq. (30), that is, the second and
, obtained at the stationary state with the irradiance F = 5 GW/cm2 .
http://WWW.CHEMISTRYVIEWS.COM
FULL PAPER
http://WWW.C-CHEM.ORG
third columns of Table 1, is due to the fact that the photoisomerization is not instantaneous. In deriving eq. (30), as in most kinetic models,[9, 22] the excited state lifetimes are neglected. Actually, during the nonadiabatic dynamics that follows the excitation, the further absorption (or emission) of a photon would be most likely ineffective in increasing the photoisomerization rate, and in our model we have chosen to ignore completely this possibility (see “Model and Connection with Single Chromophore Molecular Dynamics” section). Therefore, the small fraction of molecules that are undergoing a photochemical process cannot start a new one (this fraction is ∼2% with the highest irradiance, F = 5 GW/cm2 ). As at molecular level the trans → cis photoisomerization is in the average slower than the cis → trans one (about 1.9 ps versus 1.0 ps), the simulation yields a slightly smaller fcis /ftrans ratio than eq. (30). This discussion shows that certain simplifying assumptions of our model affect the results in a physically sound way and however to a very small extent. Figure 4 shows the reduced densities obtained with the irradiance F = 5 GW/cm2 at the stationary state (1000 ps, about 60 reduced time units). As already observed in “Model and Connection with Single Chromophore Molecular Dynamics” section there is no dependence on the α angle if only Z-polarized light is used. To understand how the density depends on β and γ , one needs to know the direction of the n → π ∗ transition dipole in the molecular frame. In trans-azobenzene, the n → π ∗ transition is symmetry forbidden at the equilibrium C2h geometry, but it borrows intensity from the lowest π → π ∗ state. The orientation of the transition dipole is then well defined: it lies in the xy plane, almost parallel to the long axis of the molecule, and it makes an angle of about 50◦ with the xˆ axis.[23, 24] At β = 0◦ and 180◦ , the xy plane is orthogonal to the light polarˆ therefore at these angles we find the maximum ization axis Z, (β) of the ρI density. For intermediate β angles, the orthogonality of the transition dipole with respect to Zˆ can be achieved by a rotation around the molecular zˆ axis, to γ 130◦ or 310◦ (γ ) (i.e., −50◦ ), where we find the maxima of ρI . The angular distribution of cis-azobenzene is less easily interpreted, because of two reasons: first, as cis is the less abundant isomer, it is strongly influenced by the photoisomerization of trans (this is (β) why the ρI density of cis is somewhat complementary to that of trans); second, the direction of its n → π ∗ transition dipole fluctuates because of the interconversion of two conformational enantiomers.[24–26] Although the investigation of the molecular mechanisms of azobenzene photo-orientation are out of the scope of this work, this short discussion of the angular distributions suffices to show that our model can yield a detailed description of the anisotropy, beyond the mere optical parameter.
Conclusions The model we have presented describes the kinetics of photoinduced anisotropy in a sample of molecular chromophores, taking into account their photoisomerization. The model makes use in a stochastic way of the statistical data that can be obtained by simulations of the dynamics of a single molecule, especially concerning the photoinduced reorientation (with or without photoisomerization) and the rotational diffusion. Such http://onlinelibrary.wiley.com
molecular scale processes are connected for the first time to the development of anisotropy in the bulk. The main limitation of the model is that it can only be applied to diluted solutions in viscous liquids or amorphous isotropic solids, because the chromophores are assumed not to interact among them. There are no limitations as to the time scale of the simulation, which is physically related to the rate of rotational diffusion. The effect of light with different intensities and polarizations, and of sequences of pulses, could be easily simulated. A test on azobenzene shows that even in such a simple system the photoinduced anisotropy, the photoisomerization, and the other molecular processes are coupled in a non trivial way. In particular, we find that the photoisomerization of molecules selected by the exciting light according to their orientation can affect in opposite directions the anisotropy of the two isomers, the reactant and the product. As a result, the time dependence of the anisotropy can be non monotonic.
Appendix Here, we describe an algorithm to partition the space spanned by the Euler angles α, β, and γ into a number of threedimensional (3D) cells of equal volume. Our algorithm is similar to one recently put forward by Leopardi,[27] but is simpler and limited to three dimensions, while Leopardi deals with hyperspheres in an arbitrary number of dimensions. The partition we propose is more flexible, in the sense that the number and therefore the width of the cells in each of the three coordinates can be arbitrarily chosen, while Leopardi minimizes the Euclidean diameter of each cell. The angle γ can be treated independently. Its range is [0, 2π ] and it can be partitioned into Nγ segments of length γ = 2π/Nγ , numbered with the index K = 0, . . . Nγ − 1. We shall place the center of the K -th segment at γ = K γ . The problem is thus reduced to the definition of 2D cells in the {α, β} space, that is, the tessellation of a spherical surface. α and β are the azimuthal and zenithal angles, more often called φ and θ, respectively. Each 2D cell will have a trapezoidal shape and its boundaries will be segments of parallels and meridians. An exception are the two polar cells, each one being a spherical cap with center at β = 0 or π. Besides the polar caps, we have Nβ belts, each enclosed between two parallels and containing a variable number of cells. The J-th belt has a latitudinal width βJ and its NJ cells have a longitudinal width αJ = 2π/NJ . The dependence of the βJ width on the belt index J is minimized by our algorithm, but cannot be totally eliminated because of the constraint of equal surface for all the 2D cells. We shall adopt the following definitions: A. Nβ ≡ number of latitudinal belts, not including the polar caps. The latter are numbered J = 0 and J = Nβ + 1. The partitions of the north and south hemispheres are equivalent. If Nβ is odd, the equator is the midline of the central belt, otherwise the equator is the boundary between two belts. B. Nα ≡ number of cells in the equatorial belts (those having the equator as their midline or as a boundary). Journal of Computational Chemistry 2012, 33, 1015–1022
1021
FULL PAPER
http://WWW.C-CHEM.ORG
C. NC ≡ total number of 2D cells in the α, β space.
Acknowledgments
D. σ = 4π/NC ≡ surface of a 2D cell.
The authors are grateful to Piet Van Leuven and Christian Silvio Pomelli for helpful discussions.
E. βP = boundary of the north polar cap (symmetrically, the boundary of the southern cap is π − βP ). F. σcap (βP ) = 2π [1 − cos(βP )] = surface of the polar caps. G. βN,J ≡ northern (lower β) boundary of the J-th belt. H. βS,J ≡ southern (upper β) boundary of the J-th belt. I. σbelt,J = 2π [cos(βS,J ) − cos(βN,J )] ≡ surface of the J-th belt cos(β )+cos(β ) N,J S,J J. βC,J = arccos ≡ baricentric β of the J-th belt 2 (in this way we have the same surface in the two halves of the belt, above and below βC,J ). For the north cap βC,0 = 0 and for the south one βC,Nβ +1 = π . Nα and Nβ , as well as Nγ , are input data of the algorithm. In many problems where latitude is more important than longitude, one may choose Nβ > Nα . The algorithm first computes an approximate partition of the sphere in belts, the number of cells NJ for each belt, the total number of 2D cells NC and the area of the cells σ . Then, it proceeds to define the exact boundaries of the belts that contain NJ cells of the wanted surface. In detail, we go through the following steps: 1. The width βS,J −βN,J of each belt is approximated as π/(Nβ +1), and the radius of the polar caps as half of it: βP π/(2Nβ + 2).
Keywords: photo-orientation • stochastic chemistry models • photoinduced anisotropy • photoisomerization • azobenzene How to cite this article: V. Cantatore, G. Granucci, M. Persico, J. Comput. Chem. 2012, 33, 1015–1022. DOI: 10.1002/jcc.22931 [1] C. Cojocariu, P. Rochon, Pure. Appl. Chem 2004, 76, 1479. [2] K. G. Yager, C. J. Barrett, J. Photochem. Photobiol. A 2006, 182, 250. [3] I. Vecchi, A. Arcioni, C. Bacchiocchi, G. Tiberio, P. Zanirato, C. Zannoni, J. Phys. Chem. B 2007, 111, 3355. [4] D. Rais, Y. Zakrevskyy, J. Stumpe, S. Nešpøurek, Z. Sedláková, Opt. Mater 2008, 30, 1335. [5] K. Oki, Y. Nagasaka, Colloids. Surf. A Physicochem. Eng. Asp. 2009, 333, 182. ˇ [6] V. Domenici, G. Ambrožic, M. Copic, A. Lebar, I. Drevenšek-Olenik, P. Umek, B. Zalar, B. Zupanˇciˇc, M. Žigon, Polymer 2009, 50, 4837. [7] M. Devetak, B. Zupanˇciˇc, A. Lebar, P. Umek, B. Zalar, V. Domenici, G. ˇ Ambrožic, M. Žigon, M. Copic, I. Drevenšek-Olenik, Phys. Rev. E 2009, 80, 050701. [8] T. Yoshino, M. Kondo, J. Mamiya, M. Kinoshita, Y. Yu, T. Ikeda, Adv. Mater 2010, 22, 1361. [9] P. Van Leuven, C. Cantatore, M. Persico, Phys. Chem. Chem. Phys, in press, DOI: 10.1039/c2cp23128j. [10] J. R. Lakowicz, Principles of Fluorescence Spectroscopy; Springer: New York, 2006. [11] P. Cattaneo, M. Persico, J. Am. Chem. Soc. 2001, 123, 7638. [12] M. Böckmann, N. L. Doltsinis, D. Marx, Phys. Rev. E 2008, 78, 036101. [13] M. Böckmann, N. L. Doltsinis, D. Marx, J. Phys. Chem. A 2010, 114, 745.
2. The approximate central β of the J-th belt is then set to βC,J = πJ/(Nβ + 1).
[14] G. Tiberio, L. Muccioli, R. Berardi, C. Zannoni, ChemPhysChem 2010, 11, 1018. [15] T. Cusati, G. Granucci, M. Persico, J. Am. Chem. Soc 2011, 133, 5109.
3. The number of cells contained in the J-th belt is definitively set to NJ = Nα sin(βC,J ), by rounding to the nearest integer, because the surface σbelt,J is approximately proportional to sin(βC,J ). For the polar caps N0 = NNβ +1 = 1. 4. The total number of 2D cells is then NC =
1022
Nβ +1 J=0
NJ .
[16] R. N. Zare, Angular Momentum: Understanding Spatial Aspects in Chemistry and Physics, J. Wiley (1988). [17] L. Creatini, T. Cusati, G. Granucci, M. Persico, Chem. Phys. 2008, 347, 492. [18] B. J. Berne, J. Chem. Phys. 1975, 62, 1154. [19] G. Lévi, J. P. Marsault, F. Marsault-Hérail, R. E. D. McClung, J. Chem. Phys. 1975, 63, 3543. [20] G. Gauglitz, S. Hubig, J. Photochem. 1985, 30, 121.
5. The surface of a 2D cell is computed as σ = 4π/NC .
[21] H. Rau, In Photochromism: Molecules and Systems; H. Dürr, H. BouasLaurent Eds. 2003.
6. The radius of the polar caps, that is, the southern boundary of the north cap, is computed as βP = βS,0 = arccos(1 − σ/2π).
[22] M. Persico, P. Van Leuven, Phys. Chem. Chem. Phys. 2009, 11, 8433.
7. Starting from J = 1 and going incrementally to J = Nβ : - we put βN,J = βS,J−1 ; - we compute σbelt,J = NJ σ ; - we compute βS,J = arccos[cos(βN,J ) − σbelt,J /2π ], as follows from definition (I); - we compute βC,J from definition (J).
[24] T. Cusati, G. Granucci, E. Martínez-Núñez, F. Martini, M. Persico, S. Vázquez, J. Phys. Chem. A, in press, DOI: 10.1021/jp208574q.
A 3D cell is identified with the indexes I, J, and K of the angles α, β, and γ . The center of the cell is put at αI = 2πI/NJ , βJ = βC,J , and γK = K γ . The volume of a 3D cell is V = σ γ = 8π 2 /(NC Nγ ).
Received: 17 October 2011 Revised: 27 November 2011 Accepted: 16 December 2011 Published online on 13 February 2012
Journal of Computational Chemistry 2012, 33, 1015–1022
[23] T. Cusati, G. Granucci, M. Persico, G. Spighi, J. Chem. Phys. 2008, 128, 194312.
[25] Y. Ootani, K. Satoh, A. Nakayama, T. Noro, T. Taketsugu, J. Chem. Phys. 2009, 131, 194306. [26] O. Weingart, Z. Lan, A. Koslowski, W. Thiel, J. Phys. Chem. Lett. 2011, 2, 1506. [27] P. Leopardi, Electron Trans. Numer. Anal. 2006, 25, 309.
http://WWW.CHEMISTRYVIEWS.COM
