Calculation of the Vibrational Frequencies of Carbon Clusters and Fullerenes with Empirical Potentials

Vibrational frequencies for carbon clusters, fullerenes and nanotubes evaluated using empirical carbon-carbon potentials are presented. For linear and cyclic clusters, frequencies evaluated with the reactive empirical bond order (REBO) potential provide the closest agreement with experiment. The mean absolute deviation (MAD) between experiment and the calculated harmonic frequencies is 79 cm 1 for the bending modes and 76 cm 1 for the stretching modes. The effects of anharmonicity are included via second order vibrational perturbation theory and tend to increase the frequency of the bending modes while the stretching modes have negative shifts in the region of 20 60 cm 1, with larger shifts for the higher frequency modes. This results in MADs for the bending and stretching modes of 84 cm 1 and 58 cm 1, respectively. For the fullerene molecule C60, the high frequency modes are predicted to have harmonic frequencies that are significantly higher than experiment, and this is not corrected by accounting for anharmonicity. This overestimation of experimental observed frequencies is also evident in the calculated frequencies of the G band in nanotubes. This suggests that the REBO potential is not optimal for these larger systems and it is shown that adjustment of the parameters within the potential leads to closer agreement with experiment, particularly if higher and lower frequency modes are considered separately.


Introduction
The prediction of the infrared (IR) and Raman spectroscopy of carbon clusters and closed carbon cages is a problem of fundamental interest and a challenge for computational methods.Carbon clusters are often studied in relation to the chemistry of carbon stars, 1,2 and IR measurements have had a prominent role in the detection of C 60 fullerene and the possible detection C 70 in a young planetary nebula. 3Furthermore, IR and Raman spectroscopy are used to characterise the structure of nanotubes, 4,5 and IR spectroscopy has been used to study charge dynamics in graphene. 6The capability to compute the vibrational frequencies and associated spectra of these systems accurately can potentially aid the interpretation and identification of fullerene species in experimental measurements and allow the relationship between the molecular structure and the observed features to be explored.From a quantum chemical perspective within the Born-Oppenheimer approximation, the calculation of vibrational frequencies and the associated vibrational modes requires solutions of the nuclear Schrödinger equation Adopting the harmonic approximation wherein V (R) is assumed to be quadratic greatly simplifies this problem, and School of Chemistry, University of Nottingham, University Park, Nottingham, NG7 2RD.E-mail: nick.besley@nottingham.ac.uk the vibrational frequencies and normal modes are obtained through diagonalization of the mass-weighted Hessian matrix.
For small to medium sized molecules, Kohn-Sham density functional theory 7 (DFT) is most commonly used to evaluated the necessary derivatives.Within this approach the calculated frequencies are usually too high, and often a uniform scaling factor is applied. 8There are a number of well established methods for going beyond the harmonic approximation, such as second-order vibrational perturbation theory (VPT2) and vibrational configuration interaction (VCI).In general, these calculations are based upon a quartic force field, and the accuracy of VPT2 with various DFT exchange-correlation functionals has been assessed and show that hybrid functionals achieve an accuracy of about 30 cm 1 . 9The limitation of these methods compared to a harmonic analysis is their computational cost.For VPT2 calculations, the majority of the computational effort is in the evaluation of the necessary third and fourth derivatives by numerical methods.The number of energy and gradient calculations required increases rapidly as the size of the system increases.This increase in computational cost is compounded by the cost of the individual energy and gradient evaluations, which also increase as the system gets larger.
1][12] One example are methods that exploit factors such as the localized nature of the vibrational modes which can be applied to certain systems to reduce the cost of the calculations. 13,14This constitutes an approximation in the solution of the nuclear part of the problem, while using DFT to solve the electronic part of the problem and provide the force field.An alternative approach, which is the focus of this work is to simplify the description of the force field through the use of empirical potentials.The evaluation of energies and gradients with empirical potentials can be computationally trivial compared to DFT calculations, and their use can make the calculation of IR spectra for large systems, that would be too expensive with DFT forces fields, computationally tractable.Of course, for such methods to be of value then this decrease in computational cost cannot result in an unacceptable reduction in accuracy.
Empirical potentials describing the carbon-carbon interaction are some of the most highly developed and widely used in chemistry and materials science.Currently, the most widely used potential for carbon is probably the reactive empirical bond order (REBO) potential. 15This potential is a development of the original Brenner potential 16 which was based on the Tersoff potential 17 and is designed to account for changes in atomic hybridisation and allow for the breaking and forming of covalent bonds.The REBO potential is a relatively short ranged potential and the TLHT potential of Takai et al. 18 is an example of a longer range potential and further potentials for carbon have been developed by Murrell and co-workers. 19,20he focus of this work is an assessment of the accuracy of such empirical potentials in the prediction of the vibrational frequencies of carbon clusters and fullerenes.This is motivated by goal of exploiting accurate empirical potentials for the calculation of IR and Raman spectra of large fullerene and nanotube systems in a computationally tractable manner.Furthermore, potentials such as the REBO potential are widely used in molecular dynamics simulations of carbon based materials.It is an open question whether a potential that is developed and tested based upon relatively small molecules will be transferable and describe systems such as fullerenes and nanotubes accurately.The calculation of vibrational frequencies probes the curvature of the potential energy surface around a minimum.Consequently, if a potential does not accurately reproduce vibrational frequencies, then the curvature of the potential at the minima is not correct, which raises doubt over any quantitative analysis of, for example, molecular dynamics trajectories using the potential.
The study of the spectroscopy of small carbon clusters has been an area of research of considerable activity and extensive reviews of the subject can be found in the literature. 21,22arbon clusters (C n ) comprising less than 10 carbon atoms have low energy linear structures, with a1 S + g ground state for odd n and 3 S g ground state for even n. 22Although high level calculations suggest that the lowest energy structures are cyclic for some n < 10. 23 Larger clusters (n > 10), are believed to have cyclic structures owing to a reduction in ring strain.Halicioglu reported harmonic vibrational frequencies from both the Brenner and TLHT potentials for linear C n (n=2-5) clusters. 24The results showed the TLHT potential to provide the closest agreement with experiment, although the predicted stretching frequencies were considerably higher than the experimental values (by up to 500 cm 1 ).The reported frequencies for the Brenner potential were even worse, and tended to be too low by several hundred wavenumbers.
The harmonic vibrational frequencies of carbon clusters and C 60 computed with the empirical potential of Murrell and coworkers have also been reported. 20The stretching modes for small linear carbon clusters were high compared with experiment and a scaling factor of 0.615 was used for these modes.
For C 60 the low frequency (less than 1000 cm 1 ) were underestimated and the higher frequency (greater than 1000 cm 1 ) modes were overestimated and it was necessary to scale the modes.Small carbon clusters of this size are accessible to direct calculation of the harmonic frequencies with quantum chemical methods, including DFT, 25 coupled cluster theory 23 and multi reference based approaches. 26These studies have shown that overall results from DFT are reasonably accurate compared to those from coupled cluster theory, 25,26 although some discrepancies have been identified where the minimum energy coupled cluster theory structure is a saddle point according to DFT. 26 Furthermore, it is probably necessary to account for both dynamical and non-dynamical correlation to accurately treat these systems. 268][29][30] The most studied fullerene is C 60 , and the IR spectrum of C 60 has four infrared active modes (T 1u ) and ten Raman active modes (two A g and eight H g ).The four active IR modes give bands at 526, 577, 1180 and 1433 cm 1 , with a ratio of intensities of 1, 0.48, 0.45, and 0.378, respectively, 27 and there has been a considerable effort to assign all of the vibrational modes. 31The Raman spectrum has bands at 267, 431, 495, 711, 775, 1101, 1251, 1427, 1470 and 1576 cm 1 , with the most intense peaks at 267, 495 and 1470 cm 1 . 28The strong covalent bonds between carbon atoms lead vibrations with predominantly tangential displacements to have higher frequencies, while the lower frequency part of the spectrum consists mainly of radial modes. 32IR spectra have been reported for other fullerenes, including C 70 , 27 C 76 29 and C 84 , 30 and the lower symmetry of these fullerenes leads to more bands being evident in the IR and Raman spectra.Raman spectroscopy is also used to probe the structure of carbon nanotubes and graphene. 33The G band is a multiple peak feature at 1540-1595 cm 1 and is an important component of the Raman spectroscopy of these systems [34][35][36] Group theory pre-dicts the Raman-active G band in achiral nanotubes consists of A 1g , E 1g and E 2g modes. 35Furthermore, in carbon nanotubes there is a radial breathing mode at 100 -400 cm 1 , where the frequency of this mode is dependent on the diameter of the nanotubes. 37,382][43][44][45][46] The work of Hands et al. 40 is illustrative of force field based studies, and assessed several existing force fields and developed a new force field wherein the 13 force constants contained within the force field were fitted to the experimental values for the Raman active vibrational modes.This type of model will naturally reproduce the data to which it was fitted but was also able to account for the full spectrum.However, there is significant difference between this type of force field and the empirical potentials considered in this work.The potentials considered here are general carbon-carbon potentials and have not been fitted based using data for C 60 or to specifically describe vibrational frequencies.Fabian 39 also used a force constant based model with parameters fitted to experimental IR and Raman data to simulate the IR spectrum of C 60 .This work developed an approach to evaluate the intensities of the IR bands based upon the bond charge model. 47 good description of the IR spectrum was achieved and it was suggested that anharmonicity provides a possible mechanism for activating weak modes resolved in IR spectra of C 60 thin films and single crystals.Calculation of harmonic frequencies of C 60 and C 70 with DFT is computationally expensive, particularly if good quality basis sets are used.2][43] The aim of this work is predominantly to assign all of the vibrational modes, although quite recently these assignments have been updated based upon inelastic neutron scattering data and periodic DFT calculations. 46In this work the authors note that while calculations normally pertain to an isolated molecule, almost all of the available data are for the solid and therefore a periodic description is necessary.Recently, 45 a self-consistent charge density-functional tight-binding method was applied to study the vibrational frequencies of some fullerenes using a harmonic treatment of the vibrations and a root mean squared deviation of about 30 cm 1 with respect to BLYP/cc-pVTZ calculations was reported.Comparison with the experimental frequencies for the IR and Raman active modes of C 60 shows that the calculated frequencies tend to be too low by up to about 80 cm 1 in the worst case.Overwhelmingly, theoretical studies of the vibrational spectroscopy of these systems have been based upon harmonic frequencies, and there has been relatively few attempts to incorporate anharmonicity into the calculations. 39,48,49A potential energy surface has been developed to describe the anharmonic vibrational mo-tions of C 60 and used to calculate anharmonically corrected fundamentals frequencies within a vibrational self-consistent field approach. 49The anharmonic frequencies were found to be within about 10 cm 1 of the harmonic frequencies, and this small degree of anharmonicity was associated with the stiff carbon-carbon bonds in C 60 .In this paper, the harmonic and anharmonic calculations of the IR spectroscopy of small carbon clusters and fullerenes calculated with empirical potentials are presented, with the accuracy of the calculations assessed through comparison with DFT based calculations, where these are computationally feasible, and experiment.

Computational Details
DFT harmonic and anharmonic frequencies were computed following full geometry optimisation using the Q-Chem software package. 50The B3LYP 51 and B97-1 exchangecorrelation functional 52 were used in conjunction withe the 6-311G* basis set 53 unless stated otherwise.Vibrational frequencies for several empirical potentials for carbon have been investigated.The first is the Murrell-Mottram 54 potential for carbon parameterised by Eggen et al. 20 The potential (denoted MM here) has the following two-body and three-body terms where and r i j = (r i j r e )/r e (9)   with r i j the distance between atoms i and j.The potential was fitted to the phonon frequencies and elastic constants of diamond and the values for the parameters (D, r e , a 2 , a 3 and c 0 c 10 ) can be found elsewhere. 20The TLHT potential 18 also has two and three-body terms, with the two-body part given by The three-body part of the potential is represented by the angle dependent term where q i , q j and q k are the angles of the triangle formed by the three atoms i, j and k, and q 1 q 5 , Z, h, p and b are parameters.In the original Brenner potential 16 (denoted Brenner), the energy is expressed as where f c (r i j ) is a cutoff function and V R and V A and b i j are defined as with and The second generation REBO potential 15 is significantly more complex than the Brenner potential and is designed to describe changes in atomic hybridization.V R and V A are modified to and g i jk is a bond-bending spline function and the P RC i j term is relevant for radicals and is zero for closed shell systems.b DH i j is a dihedral bending function that depends on the local conjugation and involves the third nearest neighbour atoms.
where Q i jkl is the dihedral angle of four atoms and T i j describes the rotation about the bonds and depends whether the atoms are conjugated.There are many parameters incorporated in the REBO potential which are described in more detail elsewhere. 15ibrational frequencies were computed for the REBO, 15 THLT 18 and Murrell 19 potentials using our own code.Structures were optimised according to the empirical potentials with the conjugate gradient technique.In this software analytical first derivatives are available for the REBO potential, otherwise the derivatives were evaluated numerically with a step size of 0.005 Å in cartesian coordinates for the first and second derivatives and 0.1 bohr along the normal modes for the third and fourth derivatives used in the evaluation of the anharmonic correction.These values of step size are typically used within quantum chemistry codes. 50For the fullerene molecules, we found it necessary to use a larger step size of 0.5 bohr for the numerical third and fourth derivatives.The smaller step size did not give a sufficiently large change in energy and the resulting anharmonic shifts were not reliable.While this sensitivity of the anharmonic shifts for the fullerenes is undesirable, implementing analytical higher order derivatives for a potential as complex as the REBO potential is not practical, and we only discuss the anharmonic shifts for the fullerenes at a qualitative rather than quantitative level.Anharmonic corrections for the vibrational frequencies were computed according to VPT2 using the formula 55 where the h's represent derivatives of the energy with respect to the normal modes, w i is the harmonic vibrational frequency of normal mode q i and there are m normal modes.While evaluating the vibrational frequencies for the empirical potentials is relatively straightforward, determining the associated intensities is more problematic.In this work we have used two approaches to estimate the infrared intensities.In the first, atomic partial charges for the structure optimised according to the empirical potential were determined in a separate DFT calculation and the relative intensities of the modes were estimated through displacement of the partial charges with respect to the normal modes.Some of the molecules studied in this work are highly symmetric and the atomic partial charges are zero.In this case, dipole derivatives computed in a separate B97-1/STO-3G calculation were combined with the computed normal modes as given by the Hessian of the empirical potential to give the intensity.However, this approach is limited in that it introduces a significant computational expense in determining the dipole derivatives.Other approaches involving the introduction of bond charges have been used to evaluate intensities for C 60 and C 70 39,47 and provide a basis for computing intensities without introducing additional expense.Anharmonic corrections for the intensities were not considered.

Linear C n Molecules
Table 1 shows the harmonic frequencies for linear C n clusters computed with the second generation REBO potential, TLHT potential, MM potential and the Brenner potential for n=3 and 4. Also shown are frequencies computed using DFT with the B97-1 and B3LYP exchange-correlation functionals and 6-311G* basis set.For odd n these clusters have a 1 S + g ground state and for even n a 3 S g ground state.The variation of the predicted frequencies between these two functionals is small and the frequencies generally lie within 20 cm 1 of each other, and both lie within about 40 cm 1 of the reported CCSD(T) frequencies. 23The values from B97-1 are marginally closer to the available experimental data and we focus on this functional in the following discussion.Before considering the calculated frequencies we will briefly discuss the optimised structures.The minimum energy structures given by the empirical potentials are more uniform than those from DFT, with all of the bond lengths equal except for the two carbon-carbon bonds at the ends of cluster.For the REBO potential the bond lengths between the central carbons is 1.348 Å and 1.300 Å for the terminal bonds, with corresponding values of 1.226 Å and 1.210 Å for the TLHT potential.The predicted bond lengths for the MM potential are significantly larger, this is a consequence of the r e parameter that is set to 1.507 Å.In contrast to the REBO and THLT potentials the bond lengths are predicted to be longer in the centre of the cluster, with values of 1.529 Å for the central bonds and 1.517 Å for the terminal bonds.The DFT calculations are consistent with the REBO and TLHT potentials and also predicts the carbon-carbon bond lengths of the central carbons to be shorter than for the two end carbons.However, for the longer carbon chains the bond lengths of the central carbons are not uniform.Furthermore, the bond lengths of the end carbons are shorter for the odd n clusters ( 1 S + g states) than the even n clus-ters ( 3 S g states) of similar size, and for both states decreases as the length of the clusters increases.For the C 6 cluster the bond length between the central carbons is 1.276 Å and the 1.305 Å for the end carbons, which are closest to the values predicted by the REBO potential.We have also explored the Tersoff and Brenner potentials.The Tersoff potential does not predict linear structures for the small carbon clusters, while the Brenner potential predicted linear structures for C 3 and C 4 with very low frequencies for the bending modes (see Table 1) and non-linear structures for the larger clusters.These potentials have a closely related bending potential (eqn.17), and the results suggest that this does not describe these small clusters well.These findings are consistent with an earlier study, 24 and we do not consider these potentials further.For C 3 , the REBO potential predicts frequencies for the stretching modes that are very close to the values from DFT.The frequency for the bending mode is significantly higher than DFT, and comparison with experiment shows the value from DFT to be more accurate.The REBO calculated frequency for the symmetric stretching mode (s g ) lies within 10 cm 1 of the experimental value, while the predicted frequency for the antisymmetric stretching mode (s u ) is considerably higher than experiment.The calculated frequencies for the stretching modes are significantly higher with the TLHT potential while the bending mode frequency is slightly lower than the REBO value, but remains higher than experiment.
The frequencies for the MM potential are even higher, 3321 cm 1 for the s u mode, although the frequency of the bending mode is close to the experimental value.This is consistent with the finding of a previous study where these frequencies were scaled by 0.615 to achieve agreement with experiment. 20The higher stretching modes frequencies as given by the TLHT potential is consistent with its shorter bond length, although it is counterintuitive for that the MM potential has the longest bond length and the highest frequencies.
For the larger clusters, the bending frequencies for the REBO potential are in closer agreement with DFT and the large error for the bending mode observed for C 3 appears to be an exception.The stretching mode frequencies from the REBO potential are reasonably similar to the DFT values, however, the TLHT potential predicts frequencies that are too high, in some cases by over 700 cm 1 and the stretching frequencies for the MM potential remain considerably higher than experiment.For C 8 the order of the two highest frequency modes are interchanged by both empirical potentials relative to the DFT calculations.For the smaller clusters, our calculated frequencies for the TLHT potential are in agreement with earlier work. 24The results also demonstrate a considerable improvement between the original and second generation REBO potentials with the original potential not predicting the correct structures for these clusters while the newer version is reasonably accurate.Some of the discrepancy between the computed frequencies and experiment can be associated with the harmonic approximation.Table 2 shows the computed anharmonic corrections to the normal mode frequencies as given by VPT2 for DFT and the REBO potential.Two values for the anharmonic shift are shown.The first includes all vibrational modes within the evaluation of the anharmonic shift, while the second excludes the bending modes from the VPT2 calculation.Near degeneracy effects can often result in VPT2 giving erroneous predictions of anharmonic shifts.This can often be corrected by removing low frequency modes from the calculation resulting in more accurate anharmonic corrections for the higher frequency modes, although by removing vibrational modes from the calculation the anharmonic shift will tend to be underestimated.The DFT calculations predict a large shift to lower frequency (-65 cm 1 ) for the s u mode and a positive shift of (+20 cm 1 ) for the s g mode.These anharmonic shifts are consistent with earlier coupled cluster calculations that reported anharmonic shifts of -51 cm 1 and +18 cm 1 for the s u and s g modes, respectively. 56Anharmonic shifts at the coupled cluster level have also been reported for C 5 where values of -7, -17, -26, and -31 cm 1 were evaluated for the four stretching modes with the bending modes excluded from the anharmonic frequency calculation. 57Again, these results are consistent with our corresponding DFT calculation.For C 3 the REBO potential gives anharmonic shifts of -19 cm 1 and -69 cm 1 for the s g and s u modes.For the s u mode this is close to the value from DFT, but the shift in the s g mode has a different sign to DFT.However, combining these anharmonic corrections with the computed harmonic frequencies gives values of 1213 cm 1 and 2066 cm 1 for the s g and s u modes, which are close to the experimental values of 1225 cm 1 and 2040 cm 1 .Although, the inclusion of anharmonicity does not account for the deviation from experiment of the REBO potential for the bending mode.
For the larger clusters (n > 5) DFT anharmonic corrections are only reported with the bending modes excluded since the inclusion of the bending modes results in clearly unreliable frequencies which is likely to be a consequence of neardegeneracy effects.For the REBO potential this is less of a problem and anharmonic shifts are given with and without the bending modes included.The calculations for the larger linear clusters do reveal some general trends in the computed anharmonic shifts.Anharmonicity tends increase the frequency of the bending modes while the stretching modes have negative shifts in the region of 20 -60 cm 1 , with larger shifts for the higher frequency modes.The calculations also show that neglecting the bending modes in the evaluation of the anharmonic correction does result in an underestimation of the anharmonic shift compared to the full anharmonic calculation.The agreement between the calculated REBO frequencies and experiment is illustrated in Figure 1.For the bending modes the mean absolute deviation (MAD) between experiment and the calculated harmonic frequencies is 79 cm 1 , and this in- creases to 84 cm 1 with the inclusion of the VPT2 anharmonic correction.This represents a significantly poorer agreement with experiment than DFT for which the MAD for the bending modes is 20 cm 1 .For the stretching modes, the REBO potential is more accurate with a MAD for the harmonic frequencies of 76 cm 1 which decreases to 58 cm 1 with the inclusion of anharmonic effects.This compares to a MAD of 63 cm 1 for the harmonic DFT frequencies, although this MAD decreases to 21 cm 1 following scaling of the frequencies by a standard value of 0.96.For these clusters, the displacement of atomic charges along the normal modes allows for the modes with non-zero intensity of be identified.

C 16 , C 60 , C 70 and Nanotubes
The calculated frequencies for the linear carbon clusters show the REBO potential to be the significantly more accurate than the other empirical potentials considered.While the computed REBO frequencies are significantly less accurate than DFT for the bending modes, its accuracy for the stretching modes is comparable to unscaled harmonic DFT frequencies and overall the decrease in accuracy is modest in view of the vastly reduced computational cost.Consequently, the use of the REBO potential has the potential to allow relatively accurate calculations of the vibrational frequencies of much larger carbon systems and we focus on the use of this potential.In order to test The optimised structures of cyclic C 16 as given by the REBO potential is a highly symmetrical planar structure with all of the bond lengths equal to 1.356 Å.The computed harmonic and anharmonic frequencies are illustrated in Figure 2 and tabulated in Table 3.In this Figure, superimposed on the computed spectra are asterisks denoting the computed frequencies with zero intensities.Also shown are the harmonic frequencies computed at the B97-1/cc-pVDZ level.In general, the lower frequency modes (< 450 cm 1 ) correspond to vibrations involving motions out of the molecular plane and the higher frequency modes involve motions within the plane of the molecule.For both REBO and DFT spectra, there is only one (doubly-degenerate) vibrational mode that has nonzero intensity.This is computed to have an overall intensity of 34.1 km mol 1 from DFT and is predicted to be slightly weaker using the combination of the REBO normal modes and B97-1/STO-3G dipole derivatives at 19.4 km mol 1 .This difference may be attributed to using a smaller basis set in the calculation of the dipole derivatives.The effect of anharmonicity in the cyclic clusters is similar to the linear clusters.Most of the bending modes are shifted to higher frequency by about 5-15 cm 1 , while the stretching modes are shifted to lower frequency by typically 30 cm 1 .Comparison with the DFT calculations can provide some estimate of the accuracy of the REBO frequencies.For the mode with non-zero intensity, the REBO potential gives harmonic and anharmonic frequencies of 595 cm 1 and 589 cm 1 which is in reasonable agreement with the DFT value of 606 cm 1 .However, for the modes with no intensity there is a overestimation of the higher frequency modes and also some underestimation of the lower frequency modes.The calculated and experimental frequencies of the IR and Raman active modes of C 60 are shown in Table 4.These modes have been assigned based on a visual comparison of the animated normal modes using IQMOL (http://iqmol.org/) with those from a DFT calculation where the IR and Raman active modes can be determined directly through calculated intensities.Overall, the lower frequency modes tend to be underestimated and the high frequency modes are significantly overestimated, which is consistent with the findings for C 16 .Furthermore, the relative order of the vibrational modes is not consistent with experiment.For example, the lowest T 1u mode is predicted to have a lower frequency that the lowest A g mode.A positive shift in the low frequency modes arising from anharmonicity would be expected based upon the calculations the smaller clusters.REBO VPT2 calculations are consistent with this with shifts in the range +10 cm 1 to +30 cm 1 for the lower frequency modes.For the higher frequency modes (>1000 cm 1 ), VPT2 anharmonic shifts are predicted in the range -10 to -20 cm 1 with some higher frequency shifts of -30 cm 1 for some modes, for example the A g mode at 1666 cm 1 in C 60 .The magnitude of these shifts are larger than the values of about -10 cm 1 suggested previously, 49 but are too small to fully account for the deviations from experiment suggesting that the improvements in the potential can be made.Nanotubes are related systems and the G band in nanotubes has been well characterised by Raman spectroscopy and lies at 1540-1595 cm 1 .Within a non-periodic framework nanotubes can be modelled by finite open or capped nanotubes. 58,59For finite uncapped zigzag nanotubes of length 30.6 Å and varying diameters including (9,0), (10,0), (11,0), (12,0) and (13,0) nanotubes (illustrated in Figure 3) the vibrational modes comprising the G band are computed to lie in the range 1730-1780 cm 1 .This is approximately 200 cm 1 too high and is inline with the findings for C 60 .The calculated intensities of the IR active modes evaluated by combining dipole derivatives from a B97-1/STO-3G calculation with the REBO normal modes gives relative intensities of 1:0.01:0.86:0.5 compared to the experimental values of 1:0.48:0.45:0.4 which shows the intensity of the second T 1u mode to be vastly underestimated.However, the source of this discrepancy can be associated with the DFT calculation which also underestimates the intensity of this band.
Alternative parameterisations of the REBO potential have been reported in the literature.In particular, Lindsay and Broido modified the potentials to improve its description of structural data and to the in-plane phonon-dispersion data for graphite, and the resulting potential also gave better lattice thermal conductivity values in single-walled carbon nanotubes. 60In this work the re-parameterisation was limited to T 0 (equation 22) and the coefficients of the bond bending spline function.Also shown in Table 4 are the com-Fig.3 Model of the (13,0) nanotube with a component of the G band vibration indicated puted frequencies from this version of the potential, denoted REBO LD .The performance of this parameterisation relative to the original REBO potential is mixed.Overall, there is an improvement in the predicted frequencies as shown by a reduction in the root-mean square (RMS) deviation from experiment.However, the predicted frequencies of the high frequency modes remain significantly overestimated and also the order of the two lowest IR active T 1u modes are incorrectly interchanged.
We have also modified the potential to optimise the values of the predicted vibrational frequencies for these modes without a significant perturbation of the underlying potential.The two parameters in the potential that have been varied are the exponent a in the repulsive part of the potential (equation 18) and the bond angle term is modified to where a scaling factor z is introduced.Variation of the two parameters a and z leads to optimum values of a = 4.74 Å 1 (compared to the original value of a = 4.7465390606595 Å 1 ) and z = 1.10, which gives an RMS error of 79 cm 1 compared with experiment.The new parameters do reduce the deviation from experiment for the higher frequency modes but within the form of the potential it was not possible to achieve this without adversely affecting the low frequency modes and the resulting parameters represent a compromise in accuracy between the low and high frequency modes.The potential has been further optimised wherein the low frequency modes (<1200 cm 1 ) and high frequency modes are treated separately.In these optimisations the T i j term in the dihedral bending function (equation 22) is scaled by a factor h in addition to the variation of a and z .For the low frequency modes values of a = 4.76 Å 1 , z = 0.80 and h = 1.8 are optimal and for the high frequency modes the respective values are 4.74 Å 1 , 1.60 and 1.5.For these parameters there is a significant decrease in the RMS error to 45 cm 1 and an improvement in the predicted ordering of the vibrational modes.This error is significantly lower than that obtained through a linear scaling of the REBO frequencies, where an optimum scaling factor of 0.95 gives an RMS error of 85 cm 1 .The optimisation of the potentials was performed on the computed harmonic frequencies and subsequent inclusion of anharmonicity leads to a mixed picture and no overall improvement in the agreement with experiment.Figure 4 shows the experimental and simulated REBO opt2 IR spectra for C 70 .In the simulated spectrum the IR bands have been broadened by representing with gaussian functions with full width at half maximum of 7 cm 1 .The simulated spectrum predicts bands in reasonable agreement with the experiment.The intensity of the weaker bands is overestimated by the calculation but the most intense band is reproduced with a calculated frequency of 1459 cm 1 compared with the experimental band at 1427 cm 1 .The computed harmonic frequencies for the highest frequency component of the G band and the radial breathing mode (RBM) for a range of zigzag nanotubes calculated with the REBO opt2 parameterisation are shown in Table 5.The predicted values for G band are much closer to the experimentally observed range of 1540-1595 cm 1 compared with the original REBO potential.Similarly, the predicted frequencies for the RBM are within the 100 -400 cm 1 range that is seen experimentally.Furthermore, strong dependence of the RBM frequency on the diameter of the nanotube is reproduced by the calculations.

Conclusions
The calculation of the vibrational frequencies and associated spectroscopy of large carbon systems lies at the limits of the size of system that can be studied with standard quantum chemical approaches.This work has explored the use of empirical potentials in the calculation of harmonic and anharmonic frequencies of carbon clusters and fullerenes.For the linear carbons clusters C 3 to C 8 , the most accurate frequencies are predicted by the REBO potential, with a MAD from experiment of 79 cm 1 for the bending modes and 76 cm 1 for the stretching modes.Incorporating anharmonicity via VPT2 does not improve the predicted frequencies of the bending modes but does reduce the MAD for the stretching modes to 58 cm 1 .For C 60 and the G band in nanotubes the predicted frequencies for the higher frequency modes are consistently much too high, and for C 60 there is a large RMS deviation tween the calculated and experimental vibrational frequencies for the IR and Raman active modes.Modification of the potential to reproduce the experimental vibrational frequencies does result in a reduction in the error in the calculated frequencies but it is found that to achieve satisfactory agreement with experiment it is necessary to consider the low frequency (<1200 cm 1 ) and high frequency modes separately.This gives an RMS error of 45 cm 1 for the IR and Raman active modes of C 60 and leads to a reasonable agreement with experiment for the IR spectrum of C 70 and the G band and RBM of nanotubes.

Acknowledements
We thank the Engineering and Physical Sciences Research Council (EPSRC) for funding (Grant No. EP/I012303).

1 Fig. 1
Fig. 1 Correlation between the calculated REBO and experimental vibrational frequencies for the linear carbon clusters.The top panel shows the bending modes and the lower panel stretching modes.Calculated harmonic frequencies are shown in red and anharmonic frequencies are shown in blue

Fig. 4
Fig. 4 Experimental (upper panel) and simulated (lower panel) IR spectra of C 70 .Experimental data adapted from reference 27

Table 3
Computed harmonic and anharmonic DFT and REBO frequencies for cyclic carbon clusters, all modes are doubly degenerate except those denoted with an asterisk

Table 4
27mputed REBO harmonic frequencies in cm 1 for the IR and Raman active modes of C 60 with different parameterisations of the REBO potential (see text for details).Experimental data.27aroot mean square deviation from experiment.

Table 5
1omputed harmonic vibrational frequencies in cm1for the G band and radial breathing mode (RBM) in model nanotubes