Structure, Energetics, and Dynamics of Pedal-Like Motion in Stilbene

Simulations show the existence of pedal-like motion at higher temperatures in agreement with the recent X-ray diffraction measurements by Ogawa and ...
0 downloads 0 Views 347KB Size
J. Phys. Chem. B 2004, 108, 17403-17411

17403

Structure, Energetics, and Dynamics of Pedal-Like Motion in Stilbene from Molecular Simulation and ab Initio Calculations N. Arul Murugan and S. Yashonath* Solid state and Structural Chemistry Unit, Indian Institute of Science, Bangalore, India 560 012 ReceiVed: July 14, 2004; In Final Form: August 17, 2004

Molecular simulations of solid stilbene in the isothermal-isobaric ensemble with variable-shape simulation are reported. The structure has been characterized by means of lattice parameters and radial distribution functions. Simulations show the existence of pedal-like motion at higher temperatures in agreement with the recent X-ray diffraction measurements by Ogawa and co-workers and several others previously. The difference in energy between the major and minor conformers and the barrier to conformational change at both the crystallographic sites have been calculated. The temperature dependence of the equilibrium constant between the two conformers as well as the rate of conversion between the conformers at the two sites have also been calculated. These are in agreement with the recent analysis by Harada and Ogawa of nonequilibrium states obtained by rapid cooling of stilbene. (Harada, J.; Ogawa, K. J. Am. Chem. Soc. 2004, 126, 3539.) An estimate of the activation energies for interconversion between the two conformers at the two sites is reported. The volume and the total intermolecular energy suggests the existence of two transitions in agreement with previous Raman phonon spectroscopic and calorimetric studies. They seem to be associated with change from order to disorder at the two sites. Ab initio calculations coupled with simulations suggest that the disorder accounts for only a small part of the observed shortening in ethylene bond length.

Introduction Understanding the structures of molecular solids and their relationship to intermolecular interactions has received considerable interest in recent times.1 The crystal structure of transstilbene was first reported by Robertson and Woodward2 in 1937 based on X-ray diffraction measurements. This was subsequently confirmed by Finder et al.3 The X-ray diffraction study of stilbene crystals by Bernstein4 at room temperature found the bond length of the bridging ethylene double bond at the two crystallographic sites to be 1.313 and 1.288 Å. This may be compared to the length of the double bond in the ethylene molecule which is 1.337 Å. These reports suggest that the ethylene bond interconnecting the two phenyl rings is shorter. The calculations based on atom-atom potentials by Bernstein and Mirsky5 on trans-stilbene at two different temperatures suggest that stilbene at site 2 is disordered. Additionally, they found that at room temperature only about 20% of the molecules are disordered in site 2, which is in agreement with 13% obtained from X-ray studies of Bouwstra et al. in the case of stilbene6 and 17% in the case of structurally similar azobenzene.7 More careful collection and analysis of the X-ray diffraction data by Ogawa and co-workers8,9 on stilbenes and stilbene derivatives with methyl substituents at the ortho, meta, and para positions confirm these findings. Interestingly, their study revealed that, although the X-ray structure of the solid suggests a shorter ethylene bond length in stilbene and their derivatives, the NMR and UV study of these compounds in solution showed no evidence in support of such bond shortening. They attributed the shortening of the bond from X-ray data to an artifact arising from the fit to the dynamic averaging of torsional motion around the C-Ph bond. This motion gives rise to a crankshaft or pedal* Corresponding author. E-mail: [email protected].

like motion of the stilbenes. Saito and Ikemoto10 carried out lattice dynamical calculations within the harmonic approximation of stilbene and found that the intramolecular torsional rotation responsible for the pedal-like motion is coupled to the overall rotation of the molecule. In a more recent investigation, Harada and Ogawa11 showed that the pedal-like motion is not restricted to stilbene at site 2. The onset of pedal-like motion at site 1 occurs at a somewhat higher temperature. Crystallographically, site 1 corresponds to the two molecules near the origin, site 2 corresponds to two other molecules displaced by (0, 1/2, 1/2). More recently, Harada and Ogawa12 have reported X-ray diffraction measurements of the major and minor conformer population at various temperatures on stilbene along with fast cooling experiments in a bid to capture the less probable conformations arising from pedal-like motion observed in reasonable number of molecules (about 15%) at site 2. There have been attempts to model theoretically the pedallike motion of stilbene. Galli et al.13 reported a molecular mechanics calculation of stilbene crystal using an MM3(92) force field.14-16 These calculations give valuable information regarding the nature of the potential energy during the pedallike motion of stilbene. However, the values obtained by them for the barrier to torsional rotation are unusually large (40.7 kJ/mol at site 1 and 63.4 kJ/mol at site 2). They are many times the kinetic energy at room temperature; in fact, these results suggest that pedal-like motion is unlikely at room temperature. Additionally, the calculations reported so far, namely, those of Bernstein and Mirsky,5 Saito and Ikemoto10 as well as those by Galli et al.,13 are essentially static calculations of dynamic phenomena. Recently, Tasumi and co-workers17,18 assigned the various vibrational modes of stilbene in ground and excited states in the frequency range 300-1700 cm-1 by the use of Raman measurements and calculations to assign the various bands. They

10.1021/jp0468859 CCC: $27.50 © 2004 American Chemical Society Published on Web 10/15/2004

17404 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Murugan and Yashonath

also investigated the distortion or torsional rotation around the CdC and C-Ph bonds in the solution phase of stilbene with benzene and n-hexane as solvents by measuring the Raman intensity of the in-phase ethylenic CH wag. More recently, Takeshita et al.19 investigated phase transitions in stilbene crystals as a function of hydrostatic pressure. Luminescence studies on stilbene by Ghoshal et al.20 and Raman phonon spectroscopic studies by Chakrabarti and Misra21 suggest the possibility of two transitions in stilbene, one at 105115 K and another at a higher temperature in the range of 215225 K. However, because of the large experimental error associated with the width of the Raman bands, these are somewhat inconclusive. Heat capacity measurements by Saito et al.22 suggest an anomalous behavior at 170 K. They concluded that the glass transition arises from the freezing of the pedallike motion. Evidence from these studies is suggestive of one or maybe two weak second-order transitions in stilbene. Many of the above-mentioned experimental methods do provide insight into the nature of molecular conformations in crystalline stilbene, although they remain indirect at best. The limitations in experimental techniques to provide a detailed microscopic view and the difficulty in understanding the relationship between the changes at the molecular level with the macroscopic changes suggest necessity of alternative methods of investigation. Additionally, the inadequacy of the molecular mechanics calculations to model an essentially dynamic phenomena have led us to carry out detailed simulations in the isothermal-isobaric ensemble (NPT) using a variable-shape simulation cell. Molecular mechanics has several serious limitations, for example, changes in lattice parameters and cooperative translational and reorientational movement of first shell of neighbors might actually lead to facile pedal-like motion. In fact, as we shall see, our results based on a realistic simulation suggest that the barriers to pedal-like motion are significantly lower than indicated by the molecular mechanics calculations. The isothermal-isobaric ensemble variable-shape simulation was first proposed by Parrinello and Rahman. 23 This along with the subsequent modification by Yashonath and Rao24 permits a realistic calculation of stilbene molecule. The beautiful aspect of these simulations is that the lattice parameters and motion at the molecular level evolve with the external conditions, viz., temperature and pressure. These simulations provide macroscopic details of the crystal structure and energetics with temperature as well as microscopic details such as molecular conformation, rotation, and translation and the associated changes in enthalpy, intra and intermolecular energy, and volume. These provide valuable insight into the nature of motion and the accompanying changes in the energetics which are usually obtainable only indirectly through experiments. Methods Intermolecular Potential. The atom-atom potential of Williams25 has been used to model intermolecular interactions in stilbene. The potential is based on experimental crystal structures of 134 hydrocarbon molecules having 5-16 carbon atoms and 8 heats of sublimation. The potential is of 6-exp type with Na ) 26 interaction sites per stilbene molecule (eq 1).

φ(rij,µν) ) -

Aµν (rij,µν)6

+ Bµν exp(-Cµνrij,µν)

(1)

Here, i and j refer to molecular indices, and µ and ν refer to atomic indices. The well depth of the potential has been reduced

Figure 1. Molecular structure of stilbene with numbering.

TABLE 1: Intermolecular and Intramolecular Potential Parameters Used in the Simulation of Stilbene intermolecular potential type

A kJ/mol Å6

B kJ/mol

C Å-1

C-C C-H H-H

1314.17 538.68 220.81

19 0861.6 44 696.64 10 467.2

3.60 3.58 3.56

intramolecular potential c0 kJ/mol

c1 kJ/mol

c2 kJ/mol

c3 kJ/mol

0.04168

1.19446

-0.36037

-0.01116

to obtain a better agreement of the simulated results with experiment. The well depth was lowered by a factor of 0.8. This is consistent with previous finding that the Williams potential leads to stronger binding and, therefore, slow dynamics.26 The modified potential parameters employed are shown in Table 1. Figure 1 shows the numbering employed for the stilbene molecule. Intramolecular Potential. We are not aware of any suitable intramolecular potential function derived experimentally or theoretically to model the pedal like motion of the stilbene molecule. We have, therefore, derived an intramolecular potential function based on ab initio calculations carried out at the Hartree-Fock level with the 6-31 g(d) basis set as implemented in Gaussian 98.27 This calculation involves a set of constrained optimizations, where the dihedral angles φ1 for rotation around C(1)-C(2) and φ2 for rotation around C(14)C(15) bond are constrained in such a way that φ1 ) -φ2. Here, the direction of the rotations is arranged such that the two phenyl rings are parallel to each other after both rotations. The ab initio energies were calculated as a function of the dihedral angle in steps of 10°. To obtain agreement between the temperature at which the onset of the pedal-like motion is observed and the experimentally observed onset temperature, we found it necessary to adjust the potential parameters. Additionally, ab initio energies generally overstimate the values of the energy barriers; therefore, the energies were scaled by a factor of 0.1. We found that this led to good agreement with the experimental temperature at which the pedal-like motion begins to occur. A function of the form given below (eq 2) was fit to these ab initio energies. The curve is shown in Figure 2. The constants c0, c1, c2, and c3 were determined by the fit. These are listed in Table 1.

Vintra(φ) ) c0 + c1(1.0 - cos(2φ)) + c2(1.0 - cos(4φ)) + c3(1.0 - cos(6φ)) (2) where φ is specified in radians. Note that the torsional motion is applied on one of the phenyl rings and that the other rotation is on the other phenyl ring. This coupled with overall rotation of the molecule as a whole leads to facile pedal-like motion in stilbene.

Pedal-Like Motion in Stilbene

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17405 reducing the 9 degrees of freedom to the essential 6. This has been done here by choosing vector a along the x axis, b in xy plane, and c in any direction. This leads to a total of 6 degrees of freedom. This permits the system to undergo solid-solid phase transformation from one space group to another. Thus, the simulation involves integration over (7N + 6) degrees of freedom. In systems such as biphenyl29 or stilbene, the amplitude of the torsional mode involving the torsional displacement of the two phenyl rings is comparable to the molecular translation or rotation modes. Hence, it is necessary to sample over the torsional degrees of freedom in the configurational phase space. Here, the total interaction energy of the stilbene crystal has two components. N

U ) Uintra + Uinter )

1

N

N

Na

Na

Vintra(φi) + ∑ ∑ ∑ ∑ φ(rij,µν) ∑ 2 i)1 j)1 µ)1 ν)1 i)1 (4)

Figure 2. Ab initio derived intramolecular potential to include pedallike motion.

Computational Details Intramolecular energy has two principal contributions: (i) conjugation energy, which stabilizes the planar conformation due to π overlap, and (ii) repulsion between ethylenic hydrogens and the ortho hydrogens, which destabilizes the planar conformation. The role of intermolecular energy in the resulting conformation of the stilbene is also important. The resulting conformation is due to these three competing interactions. Variable-Shape NPT-MC Simulations. Simulations have been carried out in isothermal-isobaric or NPT ensemble with variable-shape simulation23 cell using the Monte Carlo method and the important sampling algorithm of Metropolis.28 The simulation cell is represented by 9 degrees of freedom. Seven degrees are associated with the each molecule, three with translation, three with rotation, and one with the torsional motion involving the two dihedral angles, φ1 and φ2 (because φ1 ) -φ2). Thus, the average of any property a is obtained by integrating over these (7N + 9) variables 〈a〉 )

∫da ∫db ∫dc ∫dΩ ∫dφ ∫dr a(r , Ω , φ ) p(r , Ω , φ ) ∫da ∫db ∫dc ∫dΩ ∫dφ ∫dr p(r , Ω , φ ) N

N

N

N

N

N

N

N

N

N

N

N

N

N

N

(3) -βU(rN,ΩN,φN

where p(rN, ΩN, φN) ) e ) and a(rN, ΩN, φN) are the probability and the property, respectively, for each configuration specified by (rN, ΩN, φN). Here, rN values are the center of mass positions of the molecules, and ΩN values specify the orientations of molecules as a whole, whereas φN specifies the dihedral angles of the N molecular species. Note that r and Ω have three components each: r specifies the three Cartesian coordinates and Ω, the three Euler angles. Variable-shape simulation cell was employed because it provides the necessary degrees of freedom for the simulation cell to be in any one of the crystal systems. The simulation cell is represented by three vectors, the cell vectors a, b, and c. Parrinello and Rahman23 originally chose the three vectors without any constraints. The cell vectors contribute nine additional degrees of freedom. However, as Yashonath and Rao24 showed, subsequently, only 6 degrees of freedom are necessary to represent the simulation cell. These could be, for example, the typical cell parameters a, b, c, R, β, and γ. The three additional degrees of freedom in the original formulation of Parrinello and Rahman23 often lead to rotation of the simulation cell24 as a whole, especially in polyatomic systems, which is not desirable. This can be overcome by

Simulations have been carried out on 3 × 7 × 3 unit cells with 252 stilbene molecules. The initial configuration and the lattice parameters for the simulation at 150 K have been taken from the X-ray structure at 113 K reported by Hoekstra et al.30 (space group ) P21/a, Z ) 4.) Each Monte Carlo move (MCM) consists of a random displacement, a rotation around a randomly selected axis by a random amount, and a torsional rotation again chosen randomly around the C-Ph bonds. The way the temperature and pressure appear in the calculations is through the Boltzmann probability in eq 3. One Monte Carlo step (MCS) consists of an attempted MCM once on each of the N molecules and the attempted displacement of each of the three simulation cell vectors a, b, and c. Displacements have been adjusted to obtain an acceptance ratio of 0.40 to ensure that the sampling over the configurational phase space is efficient. For runs at temperatures greater than 150 K, we used the final configuration from the nearest lower temperature. Calculations have been carried out at temperatures of 150, 160, 170, 180, 190, 200, 225, 250, 275, 300, 325, 350, 365, and 375 K at 1 atm pressure. The simulations at a given pressure and temperature are carried out for 16 000 MCS (or about 4 × 106 MCMs), which include 4000 MCS for equilibration. A center of mass (com)-com cutoff of 13 Å has been employed. Longer runs have been made (24 000 MCS) for 250 K and higher temperatures. A set of constrained optimizations in HF/6-31 g(d) level have been carried out with various dihedral angles(φ1 and φ2, where φ1 ) - φ2) between 0-180°. Except for these torsional degrees of freedom (φ1 and φ2), all other degrees of freedom are allowed to relax during the optimization. The bond lengths are obtained from the optimized structure. Results and Discussion Structure. Table 2 lists the variation in the cell parameters with temperature along with the experimental values at 113 K obtained from the X-ray measurements of Hoekstra et al.30 Figure 3 shows a plot of the lattice parameter variation in the temperature range of 150-350 K. Also shown are the lattice parameters obtained from the more recent X-ray diffraction measurements of Ogawa and co-workers11 (obtained from the Supporting Information deposited at the American Chemical Society Publications website). The agreement of a, c, R, β, and γ is generally within 1-2% of the experimental value except for b, where the values obtained from the simulation are about

17406 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Murugan and Yashonath

Figure 5. Distribution of dihedral angle around the C-Ph bonds at various temperatures in crystalline stilbene.

Figure 3. Various lattice parameters as a function of temperature compared to the experimental (XRD) cell parameters.

Figure 4. com-com, C-C, C-H, and H-H pair correlation functions as a function of temperature.

TABLE 2: Cell Parameters, Volume, and Energies for Stilbene as a Function of Temperature at Atmospheric Pressure (Energies in kJ/mol) T K

〈Utotal〉

150.0 200.0 250.0 300.0 350.0 113a a

-59.233 -57.111 -55.126 -52.362 -48.920

‚〈Uinter〉 〈Uintra〉 -59.274 -57.185 -55.233 -52.550 -49.203

0.041 0.074 0.107 0.188 0.283

200 K, there are indications of some degree of disorder as observed by a gradual decrease in the well-defined features in the rdfs leading to broad, somewhat liquidlike rdfs. By 250 K, the well-defined peaks have almost completely disappeared. This suggests that the amount of disorder is increasing in the transstilbene crystals. In globular molecules such as methane or neopentane, similar behavior is observed. In these molecules, the presence of a plastic crystalline phase with orientational disorder is responsible for the observed disappearance of the well-defined features in the rdfs. Rotational motion of the stilbene molecule as a whole is unlikely due to the tight packing of the neighboring molecules within the trans-stilbene crystals. Nature of Disorder. Pedal Motion. Previous reports by several groups 4,7,8,11 suggest the presence of pedal-like motion. Figure 5 shows the distribution of the dihedral angle φ around the C-Ph bonds. At 150 K, these angles are found to vary over a reasonably wide range of (-80 to 80°). By 200 K, the range has increased to (-100 to 100°) degrees but the probability at (180° is still small. The presence of a nonzero population at φ ) (90° at 200 K suggests the onset of pedal-like motion. Ogawa and co-workers11 suggest that by 150 K stilbene exhibits a small amount of disorder (6% at site 2). Potential Energy Landscape and Pedal-Like Motion. We have investigated the variation of total potential energy of interaction 〈u〉 of a molecule with the rest of the lattice and its own intramolecular energy. This is defined as N

a Å

b Å

c Å

R deg

β deg

γ deg

12.37 12.46 12.52 12.59 12.61 12.29

6.15 6.21 6.31 6.39 6.53 5.66

15.58 15.64 15.65 15.77 15.94 15.48

89.78 89.67 89.83 90.10 90.17 90.0

112.11 112.35 112.69 112.86 112.86 112.0

90.12 90.25 90.25 90.11 90.36 90.0

Experimental.30

8-10% higher. The relatively high value for b is likely to affect the pedal-like motion, probably making it easier. A better agreement of b with the experimental value will ensure that the predicted behavior, including the activation energy for pedallike motion, will be more accurately predicted. Figure 4 shows the radial distribution functions (rdfs) for com-com, C-C, C-H, and H-H. At 150 K, the rdfs show well-defined features, suggesting absence of any disorder. By

u)

Na

Na

φ(rijµν) + Vintra(φ) ∑ ∑ ∑ j)1,j*i µ)1 ν)1

(5)

Here 〈〉 indicates an average over all molecules and over all MC steps after equilibration. These are plotted in Figure 6. At low temperatures, because the stilbenes are not performing pedal-like motion, there are no molecules at values of φ, for example, outside of (-80 to 80°) range. At higher temperatures, the curves extend over the whole range of φ, namely, (-180 to +180°). Thus, Figure 6 provides the underlying potential energy landscape or profile for a stilbene molecule performing pedallike motion. From the Figure, it is evident that an energetic barrier exists at φ ) 90° for pedal-like motion. The φ ) 90° conformation appears to be the transition state for pedal-like motion. Note that, at low temperatures, this barrier is different for clockwise and anticlockwise rotations. At 200 K, the energetic barrier ∆u+ a ) u(φ ) +90°) - u(φ ) +0°) and ∆ua ) u(φ )

Pedal-Like Motion in Stilbene

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17407

Figure 6. Potential energy profile of a stilbene molecule. Note the presence of barrier at φ ) (90°. The barriers are different at φ ) +90° and φ ) -90°.

TABLE 3: Barriers for Clockwise and Anticlockwise Torsional Rotations around C-Ph Bonds during Pedal Motion T K

∆u+ a

∆ua

∆uconf kJ/mol

200 300 365

6.681 6.972 6.719

8.093 7.111 6.876

6.654 3.746 3.230

-90°) - u(φ ) +0°) are 6.68 and 8.09 kJ/mol, respectively. The difference in energy between φ ) +90° and φ ) -90° conformation arises for reasons identical to those leading to difference in energy between the minor (φ ) 180°) and the major (φ ) 0°) conformers. In fact, any two conformers differing in torsional angle by 180° will have a difference in energy or enthalpy. The reason for this left-right asymmetry arises from crystal packing. In a free molecule, φ ) +90° and φ ) -90° molecular conformations or even the major and minor conformers are energetically or otherwise indistinguishable. The crystal consists of largely the major conformers, and their own packing is quite efficient optimizing the intermolecular interaction energy. However, a minor conformer is unable to pack itself efficiently among this packing of major conformers leading to a difference in energy. The barriers for clockwise and anticlockwise rotations during pedal motion and their variation with temperature are listed in Table 3. Figure 7 shows a snapshot of the stilbene molecules down the crystallographic b axis at 150 and 300 K. The molecular arrangement at 150 K clearly shows little disorder. By 300 K, the disorder is widespread with a reasonable fraction in minor conformation (φ ) 180°). A few molecules can be observed in the process of performing pedal motion with φ close to 90°. The set of first (and probably second) neighbors surrounding a molecule performing pedal-like motion slightly expand and reorient to permit the pedal-like motion of the central molecule. This was pointed out by Bernstein and Mirsky5 who found that lattice expansion plays an important role. This leads to significant lowering of the barrier for pedal like motion. This largely arises from a lowering of the potential energy of the φ ) 90° conformation and not the φ ) 0° conformation. By 350 K, the energetic barrier is almost the same for both φ ) (90° conformations, ∆ua ) 6.72 kJ/mol with no left-right asymmetry. The availability of greater free volume is responsible

Figure 7. Snapshot down the b axis at (a) 150 K and (b) 300 K, both at 1 atm pressure. (c) Molecule performing pedal-like motion (labeled as M), with its neighbors.

for absence of left-right asymmetry. For the same reason, rotation around the C-Ph bonds is relatively facile. Harada and Ogawa12 obtained van’t Hoff plots of K, where K is the ratio of population of the major to the minor conformer.

17408 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Murugan and Yashonath

Figure 8. Arrhenius plot of the (a) equilibrium constant, K, between the major and the minor conformers in crystalline stilbene. (b) Rate constant, k, for rate of conversion between the two conformers (in units of /MC step/molecule). These yield the difference in the enthalpy between the two conformers ∆H and activation energy Ea, for their interconversion (major H minor conformer).

Figure 9. (a) Variation of 〈Utotal〉 with temperature and (b) 〈Uintra〉 averaged over site 1 and 2 separately.

They obtained it by X-ray diffraction measurements at various temperatures. In our calculations, we computed K by taking the ratio

K)

n(φ < 90°) n(φ g 90°)

(6)

where n(φ < 90°) and n(φ g 90°) give the population of major and minor conformers, respectively. The van’t Hoff plot of K (where K ) K0 e-∆H/RT) as calculated from the above ratio is shown in Figure 8. The value of ∆H obtained from this is -15.59 kJ/mol with the interatomic potential employed in the calculation. We have also plotted, in Figure 8, an Arrhenius plot of the

k ) k0e-Ea/RT

(7)

against reciprocal temperature where k is the rate of conversion between major and minor conformers. Here, Ea is the activation energy for conversion between the two conformers. As pointed out by Truhlar,31 this is a phenomenological quantity and is different from the barrier obtained from the potential energy profile. The value we obtain from this is 12.46 kJ/mol. This value is considerably lower than the values obtained for either site 1 or site 2 by Galli et al.13 Site Specific Pedal-Like Motion. Starting with the early work of Bernstein and Mirsky,5 several groups, including the most recent work of Harada and Ogawa,11 found that the onset of disorder at site 2 occurs at a lower temperature than that at site 1. Site-specific properties have been calculated from our simulations. Figure 9 shows a plot of the variation of the total interaction energy Utotal. Stronger interaction of a given molecule implies closer neighbors, which, in turn, means a higher barrier for torsional and, therefore, pedal-like motion. The total interaction energy for site 2 is lower than that for site 1,

Figure 10. ln(k) vs 1/T for site 2.

suggesting that stilbene at site 2 has a stronger interaction with the rest of the lattice. This is surprising because it is at site 2 where the disorder is observed first, at a lower temperature. Figure 10 shows a plot of logarithm of k against reciprocal temperature. We obtained an activation energy of 10.38 kJ/mol for site 2. Because of very few jumps at site 1, the data available is not sufficient to estimate the activation energy at this site accurately; but, our estimate suggests that the activation energy is higher at site 1. The overall activation energy of 12.53 kJ/ mol discussed above is therefore also higher than obtained for site 2 alone. These values may be compared with the values obtained by Sironi and co-workers13 (40.7 and 63.4 kJ/mol, respectively, for site 1 and 2) by a molecular mechanics calculation. There are, unfortunately, no experimental estimates available. More importantly, the present study suggests that the activation energy for pedal motion is lower at site 2 in contrast to the results of Sironi and co-workers13 who found that Ea is lower for site 1. Second, the magnitude of the activation energy is significantly lower than those suggested by the MM calculation. The values of activation energy obtained by Sironi and co-workers13 are rather large, being several times the molecular kinetic energy at 300 K; in fact, our results suggest that pedallike motion is more likely to occur than indicated by MM calculations. Figure 11 shows the variation in the potential energy with the dihedral angle at different temperatures for stilbene during

Pedal-Like Motion in Stilbene

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17409

Figure 11. Potential energy profile during pedal-like motion for 〈u(φ)〉 for sites 1 and 2 at different temperatures (a) 150, (b) 200, (c) 250, and (d) 300 K.

pedal-like motion at site 1 and site 2. First, we noticed that at any given temperature, the molecule at site 2 performs a larger amplitude motion than the one at site 1, although the molecule at site 2 is more strongly interacting with the remaining molecules in the crystal; that is, the magnitude of u(φ) is larger for site 2. This, as we have seen, is due to the lower barrier to rotation encountered by molecules at site 2. Additionally, as already pointed out, the clockwise and anticlockwise rotation lead to different values for u(φ) at low temperatures, that is, the curves are not symmetric around φ ) 0°. This is related to the fact that the two conformers do not have the same interaction energy. The values of the barriers for rotation for site 1 and 2 at φ ) -90° are 7.3620 and 6.8597 kJ/mol, respectively, and φ ) 90° are 7.3171 and 6.6266 kJ/mol, respectively. The differences in barrier heights between molecules at site 1 and 2 are 0.5023 kJ/mol at φ ) - 90° and 0.6905 kJ/mol at φ ) 90°. This corresponds to a temperature difference of 60 (0.5023/R) and 82.065 K, respectively. In other words, if the molecules in site 2 begin to perform pedal-like motion around a temperature say, T onset , then such a motion will occur in molecules at site 1 2 ) T onset + 60 K. at T onset 1 2 The reason for the lower values obtained in the present work could arise from the fact that, in these simulations, lattice relaxation and motion, translational as well as rotational, of neighboring molecules during pedal-like motion are permitted; in fact, the translational and rotational motion of the molecules neighboring a central molecule performing pedal-like motion probably occur cooperatively. Such cooperative motion could reduce the magnitude of the barrier to pedal-like motion significantly; in fact, the torsional motion associated with the pedal-like motion leads to change in the orientation of the two phenyl rings with respect to the neighboring stilbene molecules. To offset such changes in the orientation of phenyl rings resulting from pedal-like motion, reorientation of the molecule as a whole is necessary. This explanation is consistent with the finding by Saito and Ikemoto, who found that the intramolecular torsional rotation is coupled with overall molecular rotation.10 Thus, it appears that overall molecular rotation (apart from torsional rotation) might play an important role during the pedallike motion. Role of Pedal-Like Motion in the Thermal Expansion of the Crystal. Two different simulations have been carried out. In one set of simulations, referred to as FB, we performed the

Figure 12. (a) Difference in volume between normal stilbene capable of pedal-like motion and those of stilbene without pedal-like motion at temperatures 200, 250, 300, and 350 K at atmospheric pressure. (b) Excess volume per minor conformer at temperatures 200, 250, 300, and 350 K and at atmospheric pressure.

simulations as described earlier, including molecular translation and rotation and torsional motion leading to pedal-like motion. In another set of simulations, referred to as RB, the calculations are akin to those of FB except that torsional motion and, therefore, pedal-like motion are not permitted. This was done to obtain the contribution to the expansion in volume of the unit cell arising exclusively from the pedal-like motion. Figure 12 shows a plot of the difference in unit cell volume, (VFB VRB) and the excess volume per minor conformer, (VFB - VRB)/ nminor, where VFB, VRB, and nminor are the volumes of the unit cell from FB and RB simulations and nminor is the number of minor conformations at that temperature. At low temperatures, the increase in unit cell volume due to pedal-like motion is significant, but this is not the case at high temperatures. Phase Transition in Stilbene. We have plotted (Figure 13) the variation of the derivative of the volume with respect to temperature obtained from the NPT simulations. Also plotted is the derivative of the total energy with respect to temperature. Both of the curves show two anomalies with temperature, one is at 170 K and the other is at 250 K. Although the error associated with the two derivatives may be larger than the magnitude of the observed anomaly present, the results seem to suggest possible existence of two transitions. These results are in agreement with the results obtained by Saito et al.20 They measured heat capacity of stilbene and found that a slight change of slope in the plot of CP/T versus T around 170 K. More importantly, on cooling, stilbene showed a well-defined freezing transition near 170 K. This has been attributed to freezing of disorder below 170 K, leading to glassy crystal. Because they restricted their measurements to the range of 90-240 K, they could not have found the second transition at 250 K. Luminescence intensity measurements by Ghoshal et al.20 and Raman measurements by Chakrabarti and Misra21 indicate the existence of two transitions, one near 110 K and another near 220 K. However, our results seem to agree with those of Saito et al.22 We associate the transition at 170 K with the onset of disorder at site 2 and the transition at 250 K with the onset of disorder at site 1.

17410 J. Phys. Chem. B, Vol. 108, No. 45, 2004

Figure 13. (a) Variation of V/V150 K as a function of temperature. (b) Variation of d(V/V150 K)/dT as a function of temperature. (c) Variation of d(〈Utotal〉)/dT as a function of temperature.

Because the number of molecules with disorder is small (in case of both site 1 and 2) and becasue the disorder develops gradually, there may be no thermodynamic signature of this. The previous experimental studies and the present study suggest that the transition is quasi-continuous over the temperature range of 215-225 K.21 They did not find any soft mode or any noticeable increase in phonon bandwidth when the transition is approached from below. They conclude that the transition is of order-disorder type. These conclusions are in good agreement with our findings here as well as those of Saito et al.22 who observed a slight change in the slope in the plot of specific heat with temperature. So far, no discontinuities have been found in any first-order derivative of free energy. Thus, it appears that these changes may only be detected by more sensitive measurements, such as Raman phonon bands. Ethylene Bond Length and Disorder. We have carried out ab initio electronic structure calculations to determine the effect of disorder on the ethylene bond length. The dependence of the ethylene bond length on the dihedral angle was first computed for a single molecule of stilbene by varying the two dihedral angles (with the constraint φ1 ) -φ2) so that the molecule was in various phases of the pedal motion. This was done by using Gaussian 98.27 This gave us the variation of bond length of ethylene as a function of φ. The resulting variation of bond length on the dihedral angle is shown in Figure 14. From the MC simulations, we obtained the dihedral angle for each of the 252 molecules of stilbene in the crystal at each MC step. We then obtained the corresponding bond length from the ab initio derived curve (Figure 14) for each molecule. An average over all of these bond lengths and N molecules and all MC steps yielded the average bond length of stilbene at a given temperature, which is also shown in Figure 14. The change in the bond length obtained from ab initio calculations is about 0.01 Å as compared to the experimentally observed value of 0.05 Å at site 2. If we take into account the number of molecules that are disordered within the crystal and take an average over all molecules in the lattice, then the actual observed reduction in length of the ethylene bond is still smaller

Murugan and Yashonath

Figure 14. Variation of (a) ethylenic and (b) C-Ph bond length as a function of dihedral angle obtained from ab initio calculations. Variation of (c) ethylenic and (d) C-Ph bond length as a function of temperature derived from distribution of dihedral angle obtained from MC simulation and bond length-dihedral angle correlation obtained from ab initio calculations.

(0.002 Å) in going from 150 to 300 K. These results are in good agreement with those of Ogawa et al.,9 whose study suggested that NMR and UV showed no significant shortening. They suggested that the observed shortening in the X-ray diffraction is due to an artifact of the disorder or pedal-like motion and has its origin in the fitting of structure to the measured diffraction pattern. Our results are in agreement with this view. Conclusions The results obtained here reproduce the structure of crystalline stilbene to a good accuracy. In agreement with experiment, our simulations show the onset of pedal-like motion at site 2 between 150 and 200 K. Pedal motion starts at site 1 at higher temperatures. We have obtained the enthalpy change and activation energy from van’t Hoff and Arrhenius plots, respectively. Our results suggest a larger ∆H for site 2 than indicated by experiment. The activation energy values that we obtained are significantly smaller than the values reported by previous MM calculations of Sironi and co-workers.13 An estimate of activation energy for site 2 is reported. These are expected to be more realistic than those obtained from molecular mechanics calculations. Our calculations also provide the energy profile associated with pedal motion. This suggests that stilbene at site 2 interacts relatively strongly with the remaining molecules of the lattice; however, they encounter lower barriers, which are responsible for the early onset of pedal motion at site 2. Within the accuracy afforded by the potential employed in the present study, it appears that the onset of pedal motion at sites 1 and 2 is related by T onset ) T onset + 60 K. A more accurate potential 1 needs to be developed so that these properties can be predicted more precisely. The results of the calculations indicate the possible existence of two transitions. Our analysis indicates that these are associated with the onset of disorder at the two different sites, the lower temperature transition with site 2. Additionally, it appears that

Pedal-Like Motion in Stilbene these can only be detected with more sensitive experimental techniques such as Raman bandwidth measurements. The results of Saito et al.,22 who observed a glass transition due to freezing of torsional disorder leading to a glassy crystal, indicate that such a transition is indeed possible. The Raman bandwidth measurements of Chakrabarti and Misra21 and the luminescence studies of Ghoshal et al.20 point to the existence of a higherorder transition. It is also observed that the pedal-like motion is coupled to the overall molecular rotation. This coupling appears to be required to ensure that the pedal-motion does not lead to unfavorable orientation of the molecule with its neighbors. This result is in agreement with the lattice dynamical calculations of Saito and Ikemoto,10 which suggest that such coupling indeed is observed in stilbene crystal with disorder. Gaussian calculations coupled with our molecular simulation found that the shortening of ethylenic bond length observed in X-ray diffraction is indeed an artifact. This conclusion is in agreement with the conclusions of Ogawa et al.,9 whose careful NMR and UV experiments on stilbene in the solution phase of stilbene showed no such shortening. References and Notes (1) Burgi, H. B. Faraday Discuss. 2002, 122, 41. (2) Robertson, J. M.; Woodward, I. Proc. R. Soc. London, Ser. A. 1937, 162, 568. (3) Finder, C. J.; Newton, M. G.; Allinger, N. L. Acta Crystallogr., Sect. B 1974, 31, 411. (4) Bernstein, J. Acta Crystallogr., Sect. B 1975, 31, 1268. (5) Bernstein, J.; Mirsky, K. Acta Crystallogr., Sect. A 1978, 34, 161. (6) Bouwstra, J. A.; Schouten, A.; Kroon, J. Acta Crystallogr., Sect. C 1984, 40, 428. (7) Bouwstra, J. A.; Schouten, A.; Kroon, J. Acta Crystallogr., Sect. C 1983, 39, 112. (8) Ogawa, K.; Harada, J.; Tomodo, S. Acta Crystallogr., Sect. B 1995, 51, 240. (9) Ogawa, K.; Sano, T.; Yoshimura, S.; Takeuchi, Y.; Toriumi, K. J. Am. Chem. Soc. 1992, 114, 1041.

J. Phys. Chem. B, Vol. 108, No. 45, 2004 17411 (10) Saito, K.; Ikemoto, I. Bull. Chem. Soc. Jpn. 1996, 69, 909. (11) Harada, J.; Ogawa, K. J. Am. Chem. Soc. 2001, 123, 10884. (12) Harada, J.; Ogawa, K. J. Am. Chem. Soc. 2004, 126, 3539. (13) Galli, S.; Mercandelli, P.; Sironi, A. J. Am. Chem. Soc. 1999, 121, 3767. (14) Allinger, N. L.; Yuh, Y. H.; Lii, J.-H. J. Am. Chem. Soc. 1989, 111, 8551. (15) Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 111, 8566. (16) Lii, J.-H.; Allinger, N. L. J. Am. Chem. Soc. 1989, 111, 8576. (17) Furuya, K.; Kawato, K.; Yokoyama, H.; Sakamoto, A.; Tasumi, M. J. Phys. Chem. A 2003, 107, 8251. (18) Watanabe, H.; Okamoto, Y.; Furuya, K.; Sakamoto, A.; Tasumi, M. J. Phys. Chem. A 2002, 106, 3318. (19) Takeshita, H.; Suzuki, Y.; Hirai, Y.; Nibu, Y.; Shimada, H.; Shimada, R. Bull. Chem. Soc. Jpn. 2004, 77, 477. (20) Ghoshal, S. K.; Sarkar, S. K.; Kastha, G. S. Mol. Cryst. Liq. Cryst. 1983, 91, 1. (21) Chakrabarti, S.; Misra, T. N. Bull. Chem. Soc. Jpn. 1991, 64, 2454. (22) Saito, K.; Yamamura, Y.; Kikuchi, K.; Ikemoto, I. J. Phys. Chem. Solids 1995, 56, 849. (23) Parrinello, M.; Rahman, A. Phys. ReV. Lett. 1980, 45, 1196. (24) Yashonath, S.; Rao, C. N. R. Chem. Phys. Lett. 1985, 119, 22. (25) Williams, D. E. J. Mol. Struct. 1999, 485-486, 321. (26) Yashonath, S.; Price, S. L.; McDonald, I. R. Mol. Phys. 1988, 64, 361. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A., Jr.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98, revision A.9; Gaussian, Inc.: Pittsburgh, PA, 1998. (28) Metropolis, N.; Rosenbluth, A. W.; Rosenbluth, M. N.; Teller, A. H.; Teller, E. J. Chem. Phys. 1953, 21, 1087. (29) Murugan, N. A.; Jha, P. C.; Yashonath, S.; Ramasesha, S. J. Phys. Chem. B 2004, 108, 4178. (30) Hoekstra, A.; Meertens, P.; Vos, A. Acta Crystallogr., Sect. B 1975, 31, 2813. (31) Truhlar, J. E. J. Chem. Educ. 1978, 55, 309.