- Email: [email protected]
exponential hybrid force field scaling scheme to the bond length alternation modes of polyacetylene
Chemical Physics Letters 432 (2006) 356–361 www.elsevier.com/locate/cplett
Application of the linear/exponential hybrid force field scaling scheme to the bond length alternation modes of polyacetylene Shujiang Yang, Miklos Kertesz
*
Department of Chemistry, Georgetown University, 37th and O Streets NW, Washington, DC 20057-1227, United States Received 20 August 2006; in final form 6 October 2006 Available online 19 October 2006
Abstract The two bond length alternation related backbone carbon–carbon stretching Raman active normal modes of polyacetylene are notoriously difficulty to predict theoretically. We apply our new linear/exponential scaled quantum mechanical force field scheme to tackle this problem by exponentially adjusting the decay of the coupling force constants between backbone stretchings based on their distance which extends over many neighbors. With transferable scaling parameters optimized by least squares fitting to the experimental vibrational frequencies of short oligoenes, the scaled frequencies of trans-polyacetylene and its isotopic analogs agree very well with experiments. The linear/exponential scaling scheme is also applicable to the cis-polyacetylene case. Ó 2006 Elsevier B.V. All rights reserved.
1. Introduction trans-Polyacetylene (trans-PA) is a key conjugated polymer. Vibrational spectroscopy is a fundamental tool in determining their structures and properties. However, theoretical studies [1–5] have failed to reproduce the observed two longitudinal optical (LO) modes related to bond length alternation (BLA). Fig. 1a shows these two modes with Ag symmetry and involve mainly backbone C@C and CAC stretchings and in-plane CCH bendings, observed in Raman spectra at 1457 and 1066 cm1, respectively [6]. The former is hereafter referred to as the in-phase C@C LO stretching mode and the latter as the in-phase CAC LO stretching mode in this letter. Zerbi et al. [7] used the BLA coordinate (Z in their notation) in order to semiempirically account for the conjugation length effect on the two LO modes with the concept of amplitude mode theory [8], or the so-called effective conjugation coordinate theory [9], suggesting that oligomers with more than 20 units are required to establish a more realistic force field for trans-PA. Takeuchi et al. [10] built a force field which fit well with frequency data of both trans-PA and short oli*
Corresponding author. Fax: +1 202 687 6209. E-mail address: [email protected] (M. Kertesz).
0009-2614/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.10.050
goenes. The fact that it was partly derived from experimental data on trans-PA and its isotopic analogs makes the ‘‘prediction’’ less convincing. Using density functional crystal orbital frequency calculations at B3LYP/6-31G* level with their own codes, Hirata et al. [11] accurately predicted the two LO modes at 1452 and 1070 cm1 after uniform scaling with 0.964. However, their result is puzzling. Using the oligomer extrapolation approach including up to C50H52 (HA(CH@CH)25AH) oligoenes with the Gaussian package, we reproduced all of Hirata’s reported frequencies for trans-PA except for the in-phase C@C and CAC LO stretching modes. In our reproduction of these calculations, we obtained the two modes at 1419 and 1017 cm1 at the same level of theory using the same uniform frequency scaling. Besides, this method gives large errors in predicting the LO mode frequencies for cis-polyacetylene (cis-PA), both in their report [11] and in our reproduction. Accurate prediction for the two BLA related LO modes of cis-PA has proved to be even more challenging than the trans case. Theoretical vibrational studies of cis-PA have been scarce [1,11,12] and agreement with experiments [13–16] is poor. In this Letter, we apply our recent linear/exponential scaled quantum mechanical force field scaling scheme (Lin/Exp) [17] to address this problem. This scheme was
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
a
H
m H
2 H
t-PA
b
1
H
H
H
H
H
H
H
H
357
i
3 H
H
H
j H
H
H
Fig. 2. Numbering of the backbone carbon–carbon stretching coordinates of trans-PA oligomers. Here, the number of neighboring separation between i and j is ji jj = 3.
m H
H
cis-PA
coordinates [21]. Multiple parameter scaling is applied to the elements of the Hessian in terms of chemically meaningful internal coordinates allowing transferring the scaling factors from related molecules. Linear scaling is performed on all force constants except for those off-diagonal force constants involving the couplings of different backbone CC stretchings, as in SQM, pffiffiffiffiffiffiffi F scaled si sj ; ð1Þ ¼ F unscaled i;j i;j
Fig. 1. Structures of the polymers and the normal vectors of the bond length alternation (BLA) related longitudinal optical modes of (a) transPA and (b) cis-PA.
where i = j, or at least one of them is not a backbone CC stretching. For the off-diagonal force constants involving the couplings of different backbone CC stretchings, we apply exponential scaling [17]:
established to tackle the BLA related LO modes of highly conjugated systems where neighboring backbone stretchings couple strongly with each other over a large region of 10–15 conjugated bonds based on our calculations. The Lin/Exp scaling scheme has been shown to work well for the simple but extremely delocalized case, polyyne and its oligomers. Here, we report the in-phase C@C and CAC LO stretching mode frequencies of trans-PA, cis-PA as well as their isotopic analogs using our Lin/Exp scaling scheme.
F scaled ¼ F unscaled cjijj ; i;j i;j
2. Theoretical methodology The B3-type hybrid DFT functionals have proven to be the most reliable theory for vibrational predictions among several alternatives studied, including MP2, QCISD, and various DFT at GGA or LDA levels [18,19]. However, there is no publicly tested program for solid state vibrational calculations at this level. Therefore, we used a combined molecular-polymeric theoretical approach along with the B3LYP/6-31G* theory. Hydrogen end-capped oligoenes were used as PA oligomers. The Gaussian program was used for all molecular quantum mechanical vibrational calculations. For such highly conjugated systems, the convergency of the frequencies and force field with respect to the size of oligomers is slow, varying with theoretical levels. Therefore, our molecular calculations included the longest possible oligoenes that we could afford. The calculated molecular force field was scaled with our new Lin/Exp scaling scheme [17]. It uses the concept of Pulay’s scaled quantum mechanical (SQM) force field scaling methodology [20,21]. The force constant matrix (Hessian) obtained in Cartesian coordinates from Gaussian is transformed into the representation of redundant internal
ð2Þ
i 6¼ j and both correspond to CC backbone stretchings. Here, Fi,j is the force constant involving internal coordinates i and j; si and sj are the respective scaling factors; c is the exponential scaling factor, and ji jj is the number of neighboring separation between backbone stretchings i and j, as illustrated in Fig. 2. Application of exponential scaling to the off-diagonal backbone stretching coupling force constants of the PA systems as in Eq. (2) can be justified by the intrinsic decay behavior of Fi,j as a function of ji jj. This is demonstrated in Fig. 3, where the backbone stretching coupling force constants are plotted for the longest trans-PA oligomer HA(CH@CH)25AH. Fi,j oscillates in sign and its amplitude decreases monotonically as ji jj increases. An approximate exponential relationship of jFi,jj vs. ji jj can be obtained at the B3LYP/6-31G* level. This was also proved semi-analytically [17] to be the case with the distance dependent Hu¨ckel theory of Longuet-Higgins and Salem (LHS) [22]. In this case, the coupling force constant Fi,j between bonds i and j is approximately proportional to the mutual bond–bond polarizability Pi,j. The calculated Pi,j with the LHS model also shows an exponential decay behavior as a function of ji jj, similar to those shown in Fig. 3 at B3LYP/6-31G* level. The deviation from the exponential decay mainly comes from localized r electrons. The leading and delocalized contributions from p electrons dominate all force constant couplings beyond the first neighbor. The exponential scaling is also illustrated in Fig. 3 for comparison, showing the effect of distance-dependent scaling. It appears that such a distance dependent modification of the coupling force constants can effectively correct the over-estimation of correlation effect with the
358
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
a
0.8
0.6
Fij (mdyn/Å)
Unscaled 0.4
Lin/Exp scaled
0.2
0.0 0
4
8
12
-0.2
-0.4
|i-j|
b
0
Unscaled: R2 = 0.9714
ln |Fij|
-2
-4 Lin/Exp scaled: R2 = 0.9902 -6
-8 0
4
8
12
|i-j| Fig. 3. The backbone stretching coupling force constants Fi,j in trans-PA as a function of the number of neighboring separation between backbone stretchings i and j, ji jj. The data are based on the B3LYP/6-31G* calculated Hessian of the HA(CH@CH)25AH oligomer. The central C@C stretching (i = 25) is used as the reference stretching coordinate which couples with other C@C and CAC stretchings at different neighbors. The optimized exponential scaling factor of c = 0.920 is used.
of PA, we included in the fitting set eight in-phase C@C and CAC LO mode frequencies taken from well-defined Raman spectra of four trans-oligoenes: trans-1,3-butadiene [23,24], trans-1,3,5-hexatriene [25,26], trans-1,3,5,7-octatetraene [27], and trans-1,3,5,7,9-decapentaene [28]. The root mean square (rms) values of the best fits are summarized in Table 1. The Lin/Exp scaling scheme is compared with six other linear scaling methods which utilize Eq. (1) only. In the first, SQM(0), we used only the transferable scaling factors derived by Pulay et al. [21]. In the other five SQM methods, we selectively optimized one or two scaling factors based on the same experimental dataset. The optimized scaling factors are listed in the second row of Table 1, including scaling factors for backbone C@C or CAC stretchings s(CC) when they are treated the same, backbone C@C stretching s(C@C) and CAC stretching s(CAC) when treated differently, backbone C@C or CAC couplings s(CC/CC), and for CACAH or C@CAH bendings s(CCH). Two of these fits appear particularly inferior. The largest rms is for SQM(0), which performs quite admirable, given the fact that it uses a universal set of parameters without further optimization. SQM(2) offers only a modest improvement since it has only one optimized scaling factor for CCH bending indicating that CC stretching is key for these modes. All other methods listed perform extremely well for the fitting dataset. These optimized scaling factors are then applied to all oligoenes. After scaling, the Hessian blocks corresponding to the central part of a model oligomer were extracted to construct the force field of PA. Periodic boundary conditions (PBC) vibrational calculations were done with the Polycalc program [17,29,30] to generate the normal mode frequencies and displacement vectors of PA. Phonon dispersion curves with respect to k are obtained, but only the k = 0 modes are IR or Raman active according to the selection rules [31]. To account for the convergency of the force field, extrapolations are made on these k = 0 mode frequencies as a function of the model oligomer size in order to obtain the converged normal mode frequencies of the polymer1. 3. Results and discussion 3.1. Predictions for trans-PA
B3LYP/6-31G* theory. By applying exponential scaling in Eq. (2), the systematic errors of ab initio calculations at a specific theoretical level can be accurately corrected with suitable scaling factors, as will be shown by the success of this method. In the application of the Lin/Exp scaling scheme to polyacetylene, we used a single scaling factor for all backbone stretching coordinates (s(CC) = s(C@C) = s(CAC)). This s(CC) parameter, together with the exponential scaling factor, c, are optimized by least squares fitting of the calculated frequencies to the available experimental data of trans-PA oligomers. All other scaling factors are taken from the transferable scaling factors of Pulay et al. [21]. Since we are aiming at tackling the BLA related LO modes
trans-Polyacetylene has four atoms in the unit cell; at k = 0 eight normal modes have non-zero frequency. In Fig. 4, these unscaled frequencies are obtained through the Polycalc program, with force fields extracted from the central part of the Hessian of different oligomers, HA(CH@CH)mAH. The convergence of the k = 0 frequencies with respect to 1/2m reflects the effect of two factors: number of neighboring couplings included in the PBC 1 An alternative method is to extrapolate the calculated frequencies of oligomers directly to infinity without PBC calculations. The two methods converge to the same frequencies, but the convergency using PBC is faster and is therefore preferable due to computational costs.
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
359
Table 1 Least squares fitting results of SQM and the Lin/Exp scaling schemes for oligomers of trans-polyacetylene Method
SQM(0)
SQM(1)
SQM(2)
SQM(3)
SQM(4)
SQM(5)
Lin/Exp
Optimized parameters
None
s(CC) = 0.893
s(CCH) = 0.910
rmsa (cm1)
13.1
3.5
8.5
s(C@C) = 0.891 s(CAC) = 0.897 3.4
s(CC) = 0.892 s(CC/CC) = 0.885 3.4
s(CC) = 0.886 s(CCH) = 0.955 2.9
s(CC) = 0.893 c = 0.920 3.1
a
The eight experimental frequencies used in the fitting procedures are from Ref. [23–28] as discussed in the text.
vibrational calculation, and the resemblance of the central part of the Hessian of a model oligomer to that of the polymer. It is apparent in Fig. 4 that six of the k = 0 modes show almost no dispersion with respect to the size of the model oligomers. The frequencies of the other two modes show strong dependence on oligomer size. These are the two BLA-related modes. Since the frequencies of these two modes converge slowly with respect to the size of the model oligomers, an extrapolation process is necessary. To make the extrapolations more reliable, we did a series of molecular vibrational calculations on various transPA oligomers, with 2m up to 50. Fig. 5 demonstrates the extrapolation process for the in-phase C@C LO stretching mode frequency of trans-PA, where a quadratic fit is applied in the large model oligomer region. The extrapolated value of 1457 cm1 at the 1/2m ! 0 limit is in excellent agreement with experiments [6]. Fig. 5 implies that the effect of the exponential scaling relative to the original 3300
2800
SQM(0) scaling is more pronounced in the longer model oligomer region, where the environment is closer to that of the polymer. The performance of the Lin/Exp scaling scheme on trans-PA as well as its two isotopic analogs is presented and statistically analyzed in Table 2. The results of the six SQM scaling methods listed in Table 1 are also included for comparison. Although the SQM scaling scheme has been very successful in the literature and in fact it also gives reasonable predictions for the other six localized k = 0 modes of trans-PA (rms = 17.1 cm1), its prediction for the two BLA related LO modes is exceptionally poor. Conversely, this implies the importance of the exponential scaling to correctly address the backbone stretching couplings that extend over a number of CC bonds. With Lin/Exp scaling, excellent agreement between theory and experiment was obtained for all isotopic cases. None of the five SQM methods with optimized scaling parameters can improve the predictions for trans-PA relative to the SQM(0) method of Pulay et al. [21]. The success of the Lin/Exp method indicates that in the prediction of the BLA related LO modes, only a distance dependent scaling is capable of correcting for the limitations of the quantum chemical method that underlies the Hessian used in the frequency calculations.
-1
Freq. (cm )
1600
2300 1580 1560
Fitting for Lin/Expsclaed results: y = 7295x 2 + 705.87x + 1457.3 R2 = 0.9979
1540 -1
Freq. (cm )
1800
1520 1500 1480
1300
1460
Fitting for SQM scaled results: y = 12301x 2 + 748.11x + 1436.5 R2 = 0.9995
1440
800 0.00
0.05
0.10
0.15
1/2m Fig. 4. Unscaled k = 0 frequency of trans-PA at the B3LYP/6-31G* level obtained from the PBC vibrational calculations using the middle part Hessian of various sized oligomers. Frequencies are presented as a function of the inverse size (1/2m) of the oligomer. The solid triangles represent the in-phase C@C and CAC LO stretching modes.
1420 0.00
0.02
0.04
0.06
0.08
0.10
1/2m
Fig. 5. The in-phase C@C LO stretching mode frequency of trans-PA predicted by the Lin/Exp and SQM(0) scaling schemes, as a function of the inverse size, 1/2m, of the model oligomers used in the periodic boundary conditions calculations. Extrapolations are based on quadratic fits as shown by the solid lines. The empty triangle represents the experimental frequency (Ref. [6]) which was not used in the fit.
360
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
Table 2 Scaling results (in cm1) of the in-phase CAC and C@C LO stretching modes of trans-PA, and its deuterated and
13
C-substitued analogs
System
SQM(0)
SQM(1)
SQM(2)
SQM(3)
SQM(4)
SQM(5)
Lin/Exp
Experiment
A(CH@CH)xA
1029 1437 836 1302 1001 1414 32.7
1015 1433 830 1291 987 1412 41.8
1027 1417 828 1296 999 1393 41.8
1016 1432 831 1290 988 1411 41.7
1018 1434 831 1293 990 1412 39.9
1013 1439 831 1291 984 1419 41.7
1067 1457 847 1334 1040 1431 6.2
1066a 1457a 848b 1340b 1045a 1434a –
A(CD@CD)xA A(13CH@13CH)xA rms (cm1) a b
Ref. [6]. Ref. [15].
Table 3 Linear/Exponential scaling results of the in-phase C@C and CAC stretching of cis-PA, and its deuterated analog System
SQM(0)
Lin/Exp scaled
Experiment
A(CH@CHACH@CH)xA
1235 1530 969 1440 17.7a (14.0)b
1235 1534 967 1444 15.9a (12.6)b
1250a 1540a 976a (978)b 1470a (1460)b –
A(CD@CDACD@CD)xA rms (cm1)
Results from the standard SQM scaling scheme, SQM(0), are also listed for comparison. a Values refer to Ref. [16]. b Values in parentheses refer to Ref. [14].
3.2. Predictions for cis-PA cis-Polyacetylene is the less stable form of PA with eight atoms in the unit cell. At k = 0, it has 20 non-zero frequency normal modes. Similar to trans-PA, all other normal modes are localized except for the two BLA related modes with Ag symmetry as shown in Fig. 1b. A series of molecular vibrational calculations on cis-PA oligomers, A(CH@CHACH@CH)mA, were performed, with 4m up to 44 carbon atoms. The Hessians of these model oligomers were scaled the same way and then extracted for PBC vibrational calculations to obtain the frequencies of cis-PA, which in turn were extrapolated to the infinite size limit. The predictions of the BLA related in-phase C@C and CAC LO stretching mode frequencies of cis-PA and its deuterated analog with the Lin/Exp scaling scheme are listed in Table 3, along with predictions from the SQM(0) scaling method for comparison. The agreement between Lin/Exp scaling and experiments in the cis-PA case is acceptable but not as good as in the trans-PA case, with the overall rms deviation of 15.9 or 12.6 cm1 relative to different experimental data. The improvement over the SQM(0) scaling (rms = 17.7 or 14.0 cm1 depending on the experimental dataset) is not as significant as in the trans-PA case. The reason might be that experiments on cis-PA are not well-defined relative to its trans counterpart [6,15]. As the less stable form of PA, cis-PA samples are inevitably contaminated with trans-PA or various cis-PA segments. In Table 3, it is apparent that the Lin/Exp predicted frequencies are smaller than the available experimental results. Based on the relationship of frequencies with respect to the size of oligomers, this might imply the
samples were predominantly composed of oligomers with intermediate conjugation length. The largest deviation of Lin/Exp scaling for cis-PA from experimental data comes from the in-phase C@C LO stretching mode of the deuterated cis-PA (1444 vs. 1470 or 1460 cm1). Our scaled calculated frequencies on oligomers, which are not given here, suggest that the 1470 cm1 Raman peak corresponds to a deuterated cis-PA oligomer with about 24 carbon atoms along the chain, and the 1460 cm1 peak corresponds to the oligomer with 32 carbon atoms, rather than to an infinite chain. In summary, the linear/exponential force field scaling scheme has been applied to the bond length alternation related in-phase C@C and CAC LO stretching modes of polyacetylene. It offers a new way to correct for the deficiencies of quantum mechanically obtained Hessians based on the exponential distance dependence of the coupling force constants. In combination with B3LYP/6-31G* DFT calculations on very long oligomers, the Lin/Exp scaling scheme provides excellent agreement with experimental spectra of trans-PA and its isotopic analogs. The predictions for cis-PA are acceptable although our calculated frequencies suggest a new interpretation for the observed Raman peaks of cis-PA, especially for the deuterated samples, which are likely due to finite cis-PA segments of the order of 24–32 CH units. Acknowledgements Support from the US National Science Foundation (Grant No. DMR-0331710) is gratefully acknowledged. We are indebted to Professor Peter Pulay for valuable dis-
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
cussions on scaled quantum mechanical vibrational calculations. References [1] H. Teramae, T. Yamabe, A. Imamura, J. Chem. Phys. 81 (1984) 3564. [2] M. Kofranek, H. Lischka, A. Karpfen, J. Chem. Phys. 96 (1992) 982. [3] S. Hirata, H. Torii, M. Tasumi, J. Chem. Phys. 103 (1995) 8964, The MP2 predictions for PA were limited by using small oligomers up to HA(CH@CH)5AH. In our repeat using oligomers up to HA(C2H2)9AH, results for the BLA modes are much worse. [4] C.Q. Wu, R.T. Fu, Z.Q. Li, Y. Kawazoe, J. Phys.: Condens. Matter 9 (1997) L351. [5] J. Miao, C.Q. Wu, X. Sun, R.T. Fu, Z.Q. Li, Y. Kawazoe, Synth. Met. 101 (1999) 314. [6] H. Takeuchi, T. Arakawa1, Y. Furukawa, I. Harada, H. Shirakawa, J. Mol. Struct. 158 (1987) 179. [7] C. Castiglioni, J.T. Lopez Navarrete, G. Zerbi, M. Gussoni, Solid State Commun. 65 (1988) 625. [8] E. Ehrenfreund, Z. Vardeny, O. Brafman, B. Horovitz, Phys. Rev. B 36 (1987) 1535. [9] M. Del Zoppo, C. Castiglioni, P. Zuliani, G. Zerbi, in: T.A. Skotheim, R.L. Elsenbaumer, J.R. Reynolds (Eds.), Handbook of Conducting Polymers, Marcel Dekker, New York, 1998, p. 765. [10] H. Takeuchi, Y. Furukawa, I. Harada, H. Shirakawa, J. Chem. Phys. 84 (1986) 2882. [11] S. Hirata, S. Iwata, J. Chem. Phys. 107 (1997) 10075.
361
[12] S. Hirata, H. Torii, M. Tasumi, Bull. Chem. Soc. Japan 69 (1996) 3089, We took the results from the so-called Force Field I. Force Field II is not converged with respect to the oligomer size. [13] S. Lefrant, L.S. Lichtmann, H. Temkin, D.B. Fitchen, D.C. Miller, G.E. Whitwell Jr., J.M. Burlitch, Solid State Commun. 29 (1979) 191. [14] H. Kuzmany, Phys. Status Solidi B 97 (1980) 521. [15] L.S. Lichtmann, A. Sarhangi, D.B. Fitchen, Solid State Commun. 36 (1980) 869. [16] L.S. Lichtmann, E.A. Imhoff, A. Sarhangi, D.B. Fitchen, J. Chem. Phys. 81 (1984) 168. [17] S. Yang, M. Kertesz, V. Zolyomi, J. Kurti, submitted for publication. [18] M.W. Wong, Chem. Phys. Lett. 256 (1996) 391. [19] A.P. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502. [20] P. Pulay, G. Fogarasi, G. Pongor, J.E. Boggs, A. Vargha, J. Am. Chem. Soc. 105 (1983) 7037. [21] J. Baker, A.A. Jarzecki, P. Pulay, J. Phys. Chem. A 102 (1998) 1412. [22] H.C. Longuet-Higgins, L. Salem, Proc. Roy. Soc. A 251 (1959) 172. [23] Y.N. Panchenko, Spectrochim. Acta A 31 (1975) 1201. [24] Y. Furukawa, H. Takeuchi, I. Harada, M. Tasumi, Bull. Chem. Soc. Jpn. 56 (1983) 392. [25] R. McDiarmid, A. Sabljic´, J. Phys. Chem. 91 (1987) 276. [26] H. Yoshida, Y. Furukawa, M. Tasumi, J. Mol. Struct. 194 (1989) 279. [27] H. Yoshida, M. Tasumi, J. Chem. Phys. 89 (1988) 2803. [28] S. Hirata, H. Yoshida, H. Torii, M. Tasumi, J. Chem. Phys. 103 (1995) 8955. [29] C.X. Cui, M. Kertesz, J. Chem. Phys. 93 (1990) 5257. [30] C.H. Choi, Polycalc version 0.9.5.06, 1997. [31] P.W. Higgs, Proc. R. Soc. A 220 (1953) 472.
Application of the linear/exponential hybrid force field scaling scheme to the bond length alternation modes of polyacetylene Shujiang Yang, Miklos Kertesz
*
Department of Chemistry, Georgetown University, 37th and O Streets NW, Washington, DC 20057-1227, United States Received 20 August 2006; in final form 6 October 2006 Available online 19 October 2006
Abstract The two bond length alternation related backbone carbon–carbon stretching Raman active normal modes of polyacetylene are notoriously difficulty to predict theoretically. We apply our new linear/exponential scaled quantum mechanical force field scheme to tackle this problem by exponentially adjusting the decay of the coupling force constants between backbone stretchings based on their distance which extends over many neighbors. With transferable scaling parameters optimized by least squares fitting to the experimental vibrational frequencies of short oligoenes, the scaled frequencies of trans-polyacetylene and its isotopic analogs agree very well with experiments. The linear/exponential scaling scheme is also applicable to the cis-polyacetylene case. Ó 2006 Elsevier B.V. All rights reserved.
1. Introduction trans-Polyacetylene (trans-PA) is a key conjugated polymer. Vibrational spectroscopy is a fundamental tool in determining their structures and properties. However, theoretical studies [1–5] have failed to reproduce the observed two longitudinal optical (LO) modes related to bond length alternation (BLA). Fig. 1a shows these two modes with Ag symmetry and involve mainly backbone C@C and CAC stretchings and in-plane CCH bendings, observed in Raman spectra at 1457 and 1066 cm1, respectively [6]. The former is hereafter referred to as the in-phase C@C LO stretching mode and the latter as the in-phase CAC LO stretching mode in this letter. Zerbi et al. [7] used the BLA coordinate (Z in their notation) in order to semiempirically account for the conjugation length effect on the two LO modes with the concept of amplitude mode theory [8], or the so-called effective conjugation coordinate theory [9], suggesting that oligomers with more than 20 units are required to establish a more realistic force field for trans-PA. Takeuchi et al. [10] built a force field which fit well with frequency data of both trans-PA and short oli*
Corresponding author. Fax: +1 202 687 6209. E-mail address: [email protected] (M. Kertesz).
0009-2614/$ - see front matter Ó 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cplett.2006.10.050
goenes. The fact that it was partly derived from experimental data on trans-PA and its isotopic analogs makes the ‘‘prediction’’ less convincing. Using density functional crystal orbital frequency calculations at B3LYP/6-31G* level with their own codes, Hirata et al. [11] accurately predicted the two LO modes at 1452 and 1070 cm1 after uniform scaling with 0.964. However, their result is puzzling. Using the oligomer extrapolation approach including up to C50H52 (HA(CH@CH)25AH) oligoenes with the Gaussian package, we reproduced all of Hirata’s reported frequencies for trans-PA except for the in-phase C@C and CAC LO stretching modes. In our reproduction of these calculations, we obtained the two modes at 1419 and 1017 cm1 at the same level of theory using the same uniform frequency scaling. Besides, this method gives large errors in predicting the LO mode frequencies for cis-polyacetylene (cis-PA), both in their report [11] and in our reproduction. Accurate prediction for the two BLA related LO modes of cis-PA has proved to be even more challenging than the trans case. Theoretical vibrational studies of cis-PA have been scarce [1,11,12] and agreement with experiments [13–16] is poor. In this Letter, we apply our recent linear/exponential scaled quantum mechanical force field scaling scheme (Lin/Exp) [17] to address this problem. This scheme was
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
a
H
m H
2 H
t-PA
b
1
H
H
H
H
H
H
H
H
357
i
3 H
H
H
j H
H
H
Fig. 2. Numbering of the backbone carbon–carbon stretching coordinates of trans-PA oligomers. Here, the number of neighboring separation between i and j is ji jj = 3.
m H
H
cis-PA
coordinates [21]. Multiple parameter scaling is applied to the elements of the Hessian in terms of chemically meaningful internal coordinates allowing transferring the scaling factors from related molecules. Linear scaling is performed on all force constants except for those off-diagonal force constants involving the couplings of different backbone CC stretchings, as in SQM, pffiffiffiffiffiffiffi F scaled si sj ; ð1Þ ¼ F unscaled i;j i;j
Fig. 1. Structures of the polymers and the normal vectors of the bond length alternation (BLA) related longitudinal optical modes of (a) transPA and (b) cis-PA.
where i = j, or at least one of them is not a backbone CC stretching. For the off-diagonal force constants involving the couplings of different backbone CC stretchings, we apply exponential scaling [17]:
established to tackle the BLA related LO modes of highly conjugated systems where neighboring backbone stretchings couple strongly with each other over a large region of 10–15 conjugated bonds based on our calculations. The Lin/Exp scaling scheme has been shown to work well for the simple but extremely delocalized case, polyyne and its oligomers. Here, we report the in-phase C@C and CAC LO stretching mode frequencies of trans-PA, cis-PA as well as their isotopic analogs using our Lin/Exp scaling scheme.
F scaled ¼ F unscaled cjijj ; i;j i;j
2. Theoretical methodology The B3-type hybrid DFT functionals have proven to be the most reliable theory for vibrational predictions among several alternatives studied, including MP2, QCISD, and various DFT at GGA or LDA levels [18,19]. However, there is no publicly tested program for solid state vibrational calculations at this level. Therefore, we used a combined molecular-polymeric theoretical approach along with the B3LYP/6-31G* theory. Hydrogen end-capped oligoenes were used as PA oligomers. The Gaussian program was used for all molecular quantum mechanical vibrational calculations. For such highly conjugated systems, the convergency of the frequencies and force field with respect to the size of oligomers is slow, varying with theoretical levels. Therefore, our molecular calculations included the longest possible oligoenes that we could afford. The calculated molecular force field was scaled with our new Lin/Exp scaling scheme [17]. It uses the concept of Pulay’s scaled quantum mechanical (SQM) force field scaling methodology [20,21]. The force constant matrix (Hessian) obtained in Cartesian coordinates from Gaussian is transformed into the representation of redundant internal
ð2Þ
i 6¼ j and both correspond to CC backbone stretchings. Here, Fi,j is the force constant involving internal coordinates i and j; si and sj are the respective scaling factors; c is the exponential scaling factor, and ji jj is the number of neighboring separation between backbone stretchings i and j, as illustrated in Fig. 2. Application of exponential scaling to the off-diagonal backbone stretching coupling force constants of the PA systems as in Eq. (2) can be justified by the intrinsic decay behavior of Fi,j as a function of ji jj. This is demonstrated in Fig. 3, where the backbone stretching coupling force constants are plotted for the longest trans-PA oligomer HA(CH@CH)25AH. Fi,j oscillates in sign and its amplitude decreases monotonically as ji jj increases. An approximate exponential relationship of jFi,jj vs. ji jj can be obtained at the B3LYP/6-31G* level. This was also proved semi-analytically [17] to be the case with the distance dependent Hu¨ckel theory of Longuet-Higgins and Salem (LHS) [22]. In this case, the coupling force constant Fi,j between bonds i and j is approximately proportional to the mutual bond–bond polarizability Pi,j. The calculated Pi,j with the LHS model also shows an exponential decay behavior as a function of ji jj, similar to those shown in Fig. 3 at B3LYP/6-31G* level. The deviation from the exponential decay mainly comes from localized r electrons. The leading and delocalized contributions from p electrons dominate all force constant couplings beyond the first neighbor. The exponential scaling is also illustrated in Fig. 3 for comparison, showing the effect of distance-dependent scaling. It appears that such a distance dependent modification of the coupling force constants can effectively correct the over-estimation of correlation effect with the
358
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
a
0.8
0.6
Fij (mdyn/Å)
Unscaled 0.4
Lin/Exp scaled
0.2
0.0 0
4
8
12
-0.2
-0.4
|i-j|
b
0
Unscaled: R2 = 0.9714
ln |Fij|
-2
-4 Lin/Exp scaled: R2 = 0.9902 -6
-8 0
4
8
12
|i-j| Fig. 3. The backbone stretching coupling force constants Fi,j in trans-PA as a function of the number of neighboring separation between backbone stretchings i and j, ji jj. The data are based on the B3LYP/6-31G* calculated Hessian of the HA(CH@CH)25AH oligomer. The central C@C stretching (i = 25) is used as the reference stretching coordinate which couples with other C@C and CAC stretchings at different neighbors. The optimized exponential scaling factor of c = 0.920 is used.
of PA, we included in the fitting set eight in-phase C@C and CAC LO mode frequencies taken from well-defined Raman spectra of four trans-oligoenes: trans-1,3-butadiene [23,24], trans-1,3,5-hexatriene [25,26], trans-1,3,5,7-octatetraene [27], and trans-1,3,5,7,9-decapentaene [28]. The root mean square (rms) values of the best fits are summarized in Table 1. The Lin/Exp scaling scheme is compared with six other linear scaling methods which utilize Eq. (1) only. In the first, SQM(0), we used only the transferable scaling factors derived by Pulay et al. [21]. In the other five SQM methods, we selectively optimized one or two scaling factors based on the same experimental dataset. The optimized scaling factors are listed in the second row of Table 1, including scaling factors for backbone C@C or CAC stretchings s(CC) when they are treated the same, backbone C@C stretching s(C@C) and CAC stretching s(CAC) when treated differently, backbone C@C or CAC couplings s(CC/CC), and for CACAH or C@CAH bendings s(CCH). Two of these fits appear particularly inferior. The largest rms is for SQM(0), which performs quite admirable, given the fact that it uses a universal set of parameters without further optimization. SQM(2) offers only a modest improvement since it has only one optimized scaling factor for CCH bending indicating that CC stretching is key for these modes. All other methods listed perform extremely well for the fitting dataset. These optimized scaling factors are then applied to all oligoenes. After scaling, the Hessian blocks corresponding to the central part of a model oligomer were extracted to construct the force field of PA. Periodic boundary conditions (PBC) vibrational calculations were done with the Polycalc program [17,29,30] to generate the normal mode frequencies and displacement vectors of PA. Phonon dispersion curves with respect to k are obtained, but only the k = 0 modes are IR or Raman active according to the selection rules [31]. To account for the convergency of the force field, extrapolations are made on these k = 0 mode frequencies as a function of the model oligomer size in order to obtain the converged normal mode frequencies of the polymer1. 3. Results and discussion 3.1. Predictions for trans-PA
B3LYP/6-31G* theory. By applying exponential scaling in Eq. (2), the systematic errors of ab initio calculations at a specific theoretical level can be accurately corrected with suitable scaling factors, as will be shown by the success of this method. In the application of the Lin/Exp scaling scheme to polyacetylene, we used a single scaling factor for all backbone stretching coordinates (s(CC) = s(C@C) = s(CAC)). This s(CC) parameter, together with the exponential scaling factor, c, are optimized by least squares fitting of the calculated frequencies to the available experimental data of trans-PA oligomers. All other scaling factors are taken from the transferable scaling factors of Pulay et al. [21]. Since we are aiming at tackling the BLA related LO modes
trans-Polyacetylene has four atoms in the unit cell; at k = 0 eight normal modes have non-zero frequency. In Fig. 4, these unscaled frequencies are obtained through the Polycalc program, with force fields extracted from the central part of the Hessian of different oligomers, HA(CH@CH)mAH. The convergence of the k = 0 frequencies with respect to 1/2m reflects the effect of two factors: number of neighboring couplings included in the PBC 1 An alternative method is to extrapolate the calculated frequencies of oligomers directly to infinity without PBC calculations. The two methods converge to the same frequencies, but the convergency using PBC is faster and is therefore preferable due to computational costs.
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
359
Table 1 Least squares fitting results of SQM and the Lin/Exp scaling schemes for oligomers of trans-polyacetylene Method
SQM(0)
SQM(1)
SQM(2)
SQM(3)
SQM(4)
SQM(5)
Lin/Exp
Optimized parameters
None
s(CC) = 0.893
s(CCH) = 0.910
rmsa (cm1)
13.1
3.5
8.5
s(C@C) = 0.891 s(CAC) = 0.897 3.4
s(CC) = 0.892 s(CC/CC) = 0.885 3.4
s(CC) = 0.886 s(CCH) = 0.955 2.9
s(CC) = 0.893 c = 0.920 3.1
a
The eight experimental frequencies used in the fitting procedures are from Ref. [23–28] as discussed in the text.
vibrational calculation, and the resemblance of the central part of the Hessian of a model oligomer to that of the polymer. It is apparent in Fig. 4 that six of the k = 0 modes show almost no dispersion with respect to the size of the model oligomers. The frequencies of the other two modes show strong dependence on oligomer size. These are the two BLA-related modes. Since the frequencies of these two modes converge slowly with respect to the size of the model oligomers, an extrapolation process is necessary. To make the extrapolations more reliable, we did a series of molecular vibrational calculations on various transPA oligomers, with 2m up to 50. Fig. 5 demonstrates the extrapolation process for the in-phase C@C LO stretching mode frequency of trans-PA, where a quadratic fit is applied in the large model oligomer region. The extrapolated value of 1457 cm1 at the 1/2m ! 0 limit is in excellent agreement with experiments [6]. Fig. 5 implies that the effect of the exponential scaling relative to the original 3300
2800
SQM(0) scaling is more pronounced in the longer model oligomer region, where the environment is closer to that of the polymer. The performance of the Lin/Exp scaling scheme on trans-PA as well as its two isotopic analogs is presented and statistically analyzed in Table 2. The results of the six SQM scaling methods listed in Table 1 are also included for comparison. Although the SQM scaling scheme has been very successful in the literature and in fact it also gives reasonable predictions for the other six localized k = 0 modes of trans-PA (rms = 17.1 cm1), its prediction for the two BLA related LO modes is exceptionally poor. Conversely, this implies the importance of the exponential scaling to correctly address the backbone stretching couplings that extend over a number of CC bonds. With Lin/Exp scaling, excellent agreement between theory and experiment was obtained for all isotopic cases. None of the five SQM methods with optimized scaling parameters can improve the predictions for trans-PA relative to the SQM(0) method of Pulay et al. [21]. The success of the Lin/Exp method indicates that in the prediction of the BLA related LO modes, only a distance dependent scaling is capable of correcting for the limitations of the quantum chemical method that underlies the Hessian used in the frequency calculations.
-1
Freq. (cm )
1600
2300 1580 1560
Fitting for Lin/Expsclaed results: y = 7295x 2 + 705.87x + 1457.3 R2 = 0.9979
1540 -1
Freq. (cm )
1800
1520 1500 1480
1300
1460
Fitting for SQM scaled results: y = 12301x 2 + 748.11x + 1436.5 R2 = 0.9995
1440
800 0.00
0.05
0.10
0.15
1/2m Fig. 4. Unscaled k = 0 frequency of trans-PA at the B3LYP/6-31G* level obtained from the PBC vibrational calculations using the middle part Hessian of various sized oligomers. Frequencies are presented as a function of the inverse size (1/2m) of the oligomer. The solid triangles represent the in-phase C@C and CAC LO stretching modes.
1420 0.00
0.02
0.04
0.06
0.08
0.10
1/2m
Fig. 5. The in-phase C@C LO stretching mode frequency of trans-PA predicted by the Lin/Exp and SQM(0) scaling schemes, as a function of the inverse size, 1/2m, of the model oligomers used in the periodic boundary conditions calculations. Extrapolations are based on quadratic fits as shown by the solid lines. The empty triangle represents the experimental frequency (Ref. [6]) which was not used in the fit.
360
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
Table 2 Scaling results (in cm1) of the in-phase CAC and C@C LO stretching modes of trans-PA, and its deuterated and
13
C-substitued analogs
System
SQM(0)
SQM(1)
SQM(2)
SQM(3)
SQM(4)
SQM(5)
Lin/Exp
Experiment
A(CH@CH)xA
1029 1437 836 1302 1001 1414 32.7
1015 1433 830 1291 987 1412 41.8
1027 1417 828 1296 999 1393 41.8
1016 1432 831 1290 988 1411 41.7
1018 1434 831 1293 990 1412 39.9
1013 1439 831 1291 984 1419 41.7
1067 1457 847 1334 1040 1431 6.2
1066a 1457a 848b 1340b 1045a 1434a –
A(CD@CD)xA A(13CH@13CH)xA rms (cm1) a b
Ref. [6]. Ref. [15].
Table 3 Linear/Exponential scaling results of the in-phase C@C and CAC stretching of cis-PA, and its deuterated analog System
SQM(0)
Lin/Exp scaled
Experiment
A(CH@CHACH@CH)xA
1235 1530 969 1440 17.7a (14.0)b
1235 1534 967 1444 15.9a (12.6)b
1250a 1540a 976a (978)b 1470a (1460)b –
A(CD@CDACD@CD)xA rms (cm1)
Results from the standard SQM scaling scheme, SQM(0), are also listed for comparison. a Values refer to Ref. [16]. b Values in parentheses refer to Ref. [14].
3.2. Predictions for cis-PA cis-Polyacetylene is the less stable form of PA with eight atoms in the unit cell. At k = 0, it has 20 non-zero frequency normal modes. Similar to trans-PA, all other normal modes are localized except for the two BLA related modes with Ag symmetry as shown in Fig. 1b. A series of molecular vibrational calculations on cis-PA oligomers, A(CH@CHACH@CH)mA, were performed, with 4m up to 44 carbon atoms. The Hessians of these model oligomers were scaled the same way and then extracted for PBC vibrational calculations to obtain the frequencies of cis-PA, which in turn were extrapolated to the infinite size limit. The predictions of the BLA related in-phase C@C and CAC LO stretching mode frequencies of cis-PA and its deuterated analog with the Lin/Exp scaling scheme are listed in Table 3, along with predictions from the SQM(0) scaling method for comparison. The agreement between Lin/Exp scaling and experiments in the cis-PA case is acceptable but not as good as in the trans-PA case, with the overall rms deviation of 15.9 or 12.6 cm1 relative to different experimental data. The improvement over the SQM(0) scaling (rms = 17.7 or 14.0 cm1 depending on the experimental dataset) is not as significant as in the trans-PA case. The reason might be that experiments on cis-PA are not well-defined relative to its trans counterpart [6,15]. As the less stable form of PA, cis-PA samples are inevitably contaminated with trans-PA or various cis-PA segments. In Table 3, it is apparent that the Lin/Exp predicted frequencies are smaller than the available experimental results. Based on the relationship of frequencies with respect to the size of oligomers, this might imply the
samples were predominantly composed of oligomers with intermediate conjugation length. The largest deviation of Lin/Exp scaling for cis-PA from experimental data comes from the in-phase C@C LO stretching mode of the deuterated cis-PA (1444 vs. 1470 or 1460 cm1). Our scaled calculated frequencies on oligomers, which are not given here, suggest that the 1470 cm1 Raman peak corresponds to a deuterated cis-PA oligomer with about 24 carbon atoms along the chain, and the 1460 cm1 peak corresponds to the oligomer with 32 carbon atoms, rather than to an infinite chain. In summary, the linear/exponential force field scaling scheme has been applied to the bond length alternation related in-phase C@C and CAC LO stretching modes of polyacetylene. It offers a new way to correct for the deficiencies of quantum mechanically obtained Hessians based on the exponential distance dependence of the coupling force constants. In combination with B3LYP/6-31G* DFT calculations on very long oligomers, the Lin/Exp scaling scheme provides excellent agreement with experimental spectra of trans-PA and its isotopic analogs. The predictions for cis-PA are acceptable although our calculated frequencies suggest a new interpretation for the observed Raman peaks of cis-PA, especially for the deuterated samples, which are likely due to finite cis-PA segments of the order of 24–32 CH units. Acknowledgements Support from the US National Science Foundation (Grant No. DMR-0331710) is gratefully acknowledged. We are indebted to Professor Peter Pulay for valuable dis-
S. Yang, M. Kertesz / Chemical Physics Letters 432 (2006) 356–361
cussions on scaled quantum mechanical vibrational calculations. References [1] H. Teramae, T. Yamabe, A. Imamura, J. Chem. Phys. 81 (1984) 3564. [2] M. Kofranek, H. Lischka, A. Karpfen, J. Chem. Phys. 96 (1992) 982. [3] S. Hirata, H. Torii, M. Tasumi, J. Chem. Phys. 103 (1995) 8964, The MP2 predictions for PA were limited by using small oligomers up to HA(CH@CH)5AH. In our repeat using oligomers up to HA(C2H2)9AH, results for the BLA modes are much worse. [4] C.Q. Wu, R.T. Fu, Z.Q. Li, Y. Kawazoe, J. Phys.: Condens. Matter 9 (1997) L351. [5] J. Miao, C.Q. Wu, X. Sun, R.T. Fu, Z.Q. Li, Y. Kawazoe, Synth. Met. 101 (1999) 314. [6] H. Takeuchi, T. Arakawa1, Y. Furukawa, I. Harada, H. Shirakawa, J. Mol. Struct. 158 (1987) 179. [7] C. Castiglioni, J.T. Lopez Navarrete, G. Zerbi, M. Gussoni, Solid State Commun. 65 (1988) 625. [8] E. Ehrenfreund, Z. Vardeny, O. Brafman, B. Horovitz, Phys. Rev. B 36 (1987) 1535. [9] M. Del Zoppo, C. Castiglioni, P. Zuliani, G. Zerbi, in: T.A. Skotheim, R.L. Elsenbaumer, J.R. Reynolds (Eds.), Handbook of Conducting Polymers, Marcel Dekker, New York, 1998, p. 765. [10] H. Takeuchi, Y. Furukawa, I. Harada, H. Shirakawa, J. Chem. Phys. 84 (1986) 2882. [11] S. Hirata, S. Iwata, J. Chem. Phys. 107 (1997) 10075.
361
[12] S. Hirata, H. Torii, M. Tasumi, Bull. Chem. Soc. Japan 69 (1996) 3089, We took the results from the so-called Force Field I. Force Field II is not converged with respect to the oligomer size. [13] S. Lefrant, L.S. Lichtmann, H. Temkin, D.B. Fitchen, D.C. Miller, G.E. Whitwell Jr., J.M. Burlitch, Solid State Commun. 29 (1979) 191. [14] H. Kuzmany, Phys. Status Solidi B 97 (1980) 521. [15] L.S. Lichtmann, A. Sarhangi, D.B. Fitchen, Solid State Commun. 36 (1980) 869. [16] L.S. Lichtmann, E.A. Imhoff, A. Sarhangi, D.B. Fitchen, J. Chem. Phys. 81 (1984) 168. [17] S. Yang, M. Kertesz, V. Zolyomi, J. Kurti, submitted for publication. [18] M.W. Wong, Chem. Phys. Lett. 256 (1996) 391. [19] A.P. Scott, L. Radom, J. Phys. Chem. 100 (1996) 16502. [20] P. Pulay, G. Fogarasi, G. Pongor, J.E. Boggs, A. Vargha, J. Am. Chem. Soc. 105 (1983) 7037. [21] J. Baker, A.A. Jarzecki, P. Pulay, J. Phys. Chem. A 102 (1998) 1412. [22] H.C. Longuet-Higgins, L. Salem, Proc. Roy. Soc. A 251 (1959) 172. [23] Y.N. Panchenko, Spectrochim. Acta A 31 (1975) 1201. [24] Y. Furukawa, H. Takeuchi, I. Harada, M. Tasumi, Bull. Chem. Soc. Jpn. 56 (1983) 392. [25] R. McDiarmid, A. Sabljic´, J. Phys. Chem. 91 (1987) 276. [26] H. Yoshida, Y. Furukawa, M. Tasumi, J. Mol. Struct. 194 (1989) 279. [27] H. Yoshida, M. Tasumi, J. Chem. Phys. 89 (1988) 2803. [28] S. Hirata, H. Yoshida, H. Torii, M. Tasumi, J. Chem. Phys. 103 (1995) 8955. [29] C.X. Cui, M. Kertesz, J. Chem. Phys. 93 (1990) 5257. [30] C.H. Choi, Polycalc version 0.9.5.06, 1997. [31] P.W. Higgs, Proc. R. Soc. A 220 (1953) 472.