Relative Raman Intensities in C_{6}H_{6}, C_{6}D_{6}, and C_{6}F_{6}: a Comparison of Different Computational Methods
Stephen D. Williams, Timothy J. Johnson*, Thomas P. Gibbons, and Christopher L. Kitchens
A. R. Smith Department of Chemistry, Appalachian State University, Boone, NC 28618
*Battelle Pacific Northwest National Laboratory, P.O. Box 999 Mailstop K888, Richland, WA 99352
Abstract: The accuracy of various computational methods (HartreeFock, MP2, CCSD, CASSCF, and several types of DFT) for predicting relative intensities in Raman spectra for C_{6}H_{6}, C_{6}D_{6}, and C_{6}F_{6} were compared. The predicted relative intensities for n_{1} and n_{2} were compared with measured relative intensities. While none of these methods excelled at this prediction, HartreeFock with a large basis set was most successful for C_{6}H_{6} and C_{6}D_{6}, while PW91PW91 was the most successful for C_{6}F_{6}.
Keywords: relative Raman intensity, ab inito, HartreeFock, MP2, DFT, CCSD, CASSCF
Introduction: The assignment of Raman spectra, even for small sized molecules, is not a simple task. Quantum mechanically computed vibrational frequencies can provide significant help in this task, since the typical errors and empirical corrections for these frequencies are well known [1]. However, frequencies may not be enough to assign the spectrum, for example when a weak fundamental is near a combination or overtone. Quantum mechanically computed Raman intensities are quite helpful in such cases, particularly if their reliabilities are known. There have been several studies of computed Raman intensities compared to experiment; these have focused on one of a few methods or basis sets, and have not attempted to survey several different types of computation methods. These studies include HartreeFock methods for ethylene [2], benzene [3], HFCs [4], and triptycene [5]; DFT methods for N_{2}, HF, and ethane [6], as well as water, ammonia, and methane [7]; and hybrid DFT methods for azulene [8], porphine [9], PAHs [10], and aniline and its radical cation [11]. While some of these (small molecule) studies compare the computed Raman intensities with measured absolute Raman activities or cross sections, most compare the computed intensity to an experimental Raman spectrum, and in most of these cases, the comparison is made to an intensityuncorrected spectrum. As McCreery [12] points out in his excellent review, not correcting for the spectrometer’s response can dramatically distort a Raman band’s intensity. It is thus difficult to gauge the theoretical to experimental intensity comparisons of some of these reports. Shinohara et al. [10] do make comparisons with corrected spectra, but they do not include the CH stretching bands.
We have chosen to compare computed relative Raman intensities with accurately measured relative intensities, using several methods including HartreeFock, correlated ab initio (MP2, CCSD, CASSCF), pure DFT (PW91PW91), and hybrid DFT (B3LYP, B3PW91, BB1K). In order to be as certain as possible that we were comparing intensities of correctly assigned modes, we considered the two a_{1g} modes of C_{6}H_{6}, C_{6}D_{6}, and C_{6}F_{6}. These modes are n_{1}, the totally symmetric CX stretching mode, and n_{2}, the totally symmetric ring CC stretch; they are readily assigned due to their strongly polarized nature. We assessed the quality of a computational method on how well it reproduces the experimental ratio of the relative intensities of these modes, I(n_{1})/I(n_{2}), since determining absolute Raman cross sections is difficult, and our hope is to provide practical advice to those who routinely use Raman spectroscopy in their work. We also show a comparison of a theoretically simulated spectrum with the corresponding corrected experimental spectrum.
Experimental: Fourier transform Raman spectra were measured with an FTRaman spectrometer using 1064 nm excitation and 2 cm^{1} resolution. 180° backscattered light was collected, typically averaging 512 interferograms. Polarized spectra were also collected, typically measured with 4 cm^{1} resolution. The spectra of the neat liquids were recorded using a Bruker FRA 106 Raman module coupled to the interferometer of a Bruker IFS 66v/S vacuum Fourier infrared spectrometer. A similar FTIR spectrometer has been previously described [13], using the methods suggested by IUPAC [14]. The 106 Raman accessory consists of a CW 1.064 μm Nd:YAG laser with a focused or unfocused laser spot and with either 90^{o} or 180^{o} backscatter spectral light collection.
As described by Simon et al. [15], a notch filter in the FRA 106 removes the Rayleigh scattered photons; the Ramanscattered light is modulated by the FTIR interferometer and focused onto a sensitive liquidnitrogen cooled Ge detector [16]. In the present case the interferometer’s wavelengths are referenced to those of the HeNe laser. These were calibrated in the nearIR by using a white light source and the 2 ← 0 overtone transition of ~22 Torr of neat carbon monoxide [17] in a 10 cm cell and adjusting the spectrometer’s measured wavelengths (via the HeNe frequency) to match those of the HITRAN database [18]. To calibrate the Raman Nd:YAG laser wavelength and thus the Raman scattered frequencies, a strong scatterer, elemental sulfur, was used and the laser wavelength empirically shifted such that the Stoke’s and antiStoke’s Raman shifts were equal to within 0.2 cm^{1} for a few of the strongest lines.
Calibration of the scattering intensity as function of wavelength is somewhat more complicated. The Raman module has a 3000 K white light source whose intensity is scattered into the spectrometer off a NISTtraceable Spectralon surface mounted at the sample position. The instrument response function as a function of wavelength is then generated by ratioing this spectrum relative to the intensity predicted by Planck’s law for a 3000 K source. Subsequent spectra are then multiplied by the (frequency dependant) scaling factor. The quality of the intensity correction was assessed by measuring the intensities in the corrected Stokes and antiStokes spectra; the relative intensities were within a few percent of the expected value ([19], equation 4.36). As McCreery has pointed out, [12] obtaining Raman intensity values that correlate even on a relative scale, i.e. that properly ratio out the instrumental response, is more challenging than first supposed.
For polarization experiments the natural polarization of the laser was used. A reducing Jacquinot stop and a calcite GlanTaylor analyzer (i.e. crosspolarizer) were placed after the Raman notch filter at the focal position of the interferometer’s collimator mirror. Perpendicular and parallel orientations were obtained by rotating the analyzer to the 0^{o} and 90^{o} positions. The method was verified by measuring the parallel and perpendicular spectra of CCl_{4} whose depolarization ratio is known [20]. The depolarization ratios of natural isotope benzene, particularly the symmetric CH stretching and ring breathing modes, have also been reported and were reproduced [20, 21].
Relative intensities of the two a_{1g} stretching modes (n_{1}, CX stretch and n_{2} , CC stretch) were determined by fitting the appropriate bands to a pseudoVoigt lineshape function [22]:
_{} (1)
In this four parameter function the peak area is A, the full width at half maximum is w, the peak center is n_{0}, and h is the fraction the Lorentzian contributes to the peak. Sigmaplot was used for the curve fitting. Ratios of the resulting areas were then compared to corresponding theoretical values for C_{6}H_{6}, C_{6}D_{6}, and C_{6}F_{6}. We also measured the spectrum of ^{13}C_{6}H_{6} but did not include this in the analysis due the extensive Fermi resonance perturbation of the CH stretching modes [23].
The computations were carried out with Gaussian98 [24] on an Intel PC and Gaussian03 [25] on an SGI altix computer. The default basis sets in Gaussian were used except for a few calculations with Truhlar’s MG3 [26] basis set. For each calculation the molecular geometry was optimized followed by a frequency calculation. The geometry was constrained to have D_{6h} symmetry during the optimization. For C_{6}F_{6} some MP2 constrained optimizations yielded structures with a single imaginary frequency associated with an outofplane motion of alternate C atoms. This problem was only observed with large Pople basis sets (6311+G(3df,2p)), and did not occur when Dunning correlation consistent basis sets [27] were used. For the HartreeFock, DFT, and MP2 calculations described here the automatic methods in Gaussian were used to compute the Raman activities (keyword: freq=Raman), but for other calculations (CCSD and CASSCF) this automatic method was not available. What appears in the Gaussian output as the Raman activity is _{}(following the notation of Long [19]). Here _{}is the square magnitude of derivative of the isotropic part of the polarizability tensor with respect to the k^{th} normal mode, and likewise _{} for the anisotropic part of the polarizability, where the derivatives are evaluated at zero displacement. (The factor of 7 is for 90° scattering; for 180° scattering a factor of 4 appears instead.) For the CCSD and CAS calculations the normal modes (included in the Gaussian freq=hpmodes output) for the two a_{1g} modes were used to generate a series of molecular structures displaced along the mass weighted normal mode (using the reduced mass reported by Gaussian). The polarizability tensor components were then computed for each of these and numerical differentiation was used to estimate the polarizability derivatives; these were then used to calculate the needed polarizability tensor invariants for the Raman activities, following standard methods [19]. This numerical method was also applied to a few HartreeFock calculations; the Raman activities derived from these pointbypoint polarizability curves differed by at most a few percent from the automatic ones from the same calculation.
The CASSCF calculations used a (6,6) CAS where the active space included the 6 p MO’s composed primarily of the 2p_{z} C atom AO’s. A series of single point calculations starting with STO3g and gradually increasing to the largest basis set was used to ensure that the correct orbitals were used in the highest level CAS; visualization of these orbitals with GaussView confirmed this.
Results:
There are two strong fundamental CH stretching modes in the Raman spectrum of
benzene; these are usually denoted by n_{1}
and n_{15} of a_{1g}
and e_{2g} symmetry [28], although other notations are sometimes used [29].
These bands overlap in spectra of room temperature samples, to the extent that
recent [29] absolute Raman scattering cross section measurements of benzene
vapor have reported only the sum of the cross sections for these modes. For
our intensity measurements we need the contribution to the total intensity for n_{1} only. The pseudoVoigt
curve fitting procedure is successful in resolving these two contributions, as
shown in Figure 1 below.
Figure 1
In this case the R^{2}
statistic was 0.9998 and the areas for n_{1}
and n_{15} are 279±1 and 105±1,
in arbitrary units. The fits for the other modes and species were not always
this good, but in all cases the R^{2} statistic was at least 0.993. A
particularly challenging fit was for n_{1}
in C_{6}F_{6}. This band is quite weak (weaker than some
overtones in the spectrum) and is complicated by the presence of an isotopic
band and a broad underlying peak. The area was found to be 2.7±0.5 with R^{2}
= 0.998. Note that the intensity for n_{1}
in C_{6}F_{6} is about 1% of that for n_{1} in C_{6}H_{6}. The fit for
this band is shown in Figure 2 below.
Figure 2
The
intensity of n_{1} in the
C_{6}F_{6} Raman spectrum is quite weak and is of similar
intensity to an overtone band that is somewhat close in frequency. Since both
of these are totally symmetric, they have the same polarization properties, so
polarization measurements did not help in confirming the assignment for n_{1}. We confirmed the
assignment of n_{1} by
noting that the nearby totally symmetric mode was at nearly exactly twice the
frequency of a forbidden (b_{2g}) fundamental (n_{7}), and that the observed isotopic shift in
frequency upon substitution of a single ^{13}C atom was in near perfect
agreement with the predicted isotopic shift for n_{1}
computed at the B3LYP/6311+g(3df) and MP2/ccpVTZ levels of theory. This
isotopic band is not observed in the spectrum of C_{6}H_{6};
its appearance in the C_{6}F_{6} spectrum is a consequence
much larger atomic mass of fluorine compared to hydrogen._{ } Our
assignment of the C_{6}F_{6} Raman spectrum is shown in Figure
3 below, and agrees well with previous assignments [30].
Figure 3
The Raman spectrum
for ^{13}C_{6}H_{6} was planned for this study but it
turned out not to be suitable for measuring intensities for n_{1}, due to Fermi resonance
with overtones. This resonance perturbs the intensities of the CH stretching
modes in a way that is not simple to correct. The Fermi resonance perturbed CH
stretching region is much more complicated for ^{13}C_{6}H_{6}
compared to C_{6}H_{6} as shown in Figure 4 below.
Figure 4
We considered C_{6}H_{6}
as a test case, since Ozkabak et al. [3] had good results with Raman
intensities computed at the HartreeFock level for this compound. Methods that
seemed relatively poor for C_{6}H_{6} were not used to predict
intensities for C_{6}D_{6} and C_{6}F_{6}. We
compared the experimental intensity ratio I(n_{1})/I(n_{2}) (ratio of areas from the
curve fitting) with the predicted intensity ratio from a given level of
theory. As mentioned previously, Gaussian does not predict the Raman
intensity, rather it predicts the Raman activity. As McCreery [12] has pointed
out, the intensity of a Raman band has contributions from three sources: a part
that is intrinsic to the scattering molecule (this is Gaussian’s Raman
activity), a part that depends on temperature (which affects the population of
the scattering vibrational state), and a part that is dependent on the exciting
laser frequency and the type of photometric detection in the Raman
spectrometer. In the notation of Long [19], this intensity for 180° scattering, Stokes shift, and power detection,
may be represented as [19, 3]:
_{} (2)
where N is the number of molecules excited by the laser, h is Planck’s
constant, _{}is the
exciting laser wavenumber, _{}is the
wavenumber of the vibrational mode, c is the speed of light, k is Boltzmann’s
constant, T is the temperature, e_{0}
is the permitivity of free space, P is the laser irradiance, _{}is the
squared derivative of the isotropic part of the polarizability tensor with
respect to the k^{th} normal mode, evaluated at zero displacement, and _{} is the
squared derivative of the anisotropic part of the polarizability tensor with
respect to the k^{th} normal mode, evaluated at zero displacement.
This is the appropriate formula for our instrument. Some common modifications
might be: if a photon counting detector is used, the _{}should
be replaced with _{}[12],
and if 90° scattering is measured, the
4 should be replaced with a 7 [19]. Since we have restricted our study to the
two a_{1g} stretching modes (since we can be very confident in their
correct assignment), the _{} term is
unimportant since it is zero for these totally symmetric vibrational modes;
scattering geometry is not a factor for our study. Hence for every method
studied, we converted the computed Raman activity into a corresponding
intensity appropriate to our experimental conditions, then calculated the I(n_{1})/I(n_{2}) intensity ratio.
The results for
this computed prediction of intensity ratio are presented in Table 1 below.
Table 1
Method 
Basis Set 
I(n_{1})/I(n_{2}) 
Experimental 

0.91±0.01 
RHF 
6311+g(3df,2p) 
1.037 
RHF 
6311++g(3df,2p) 
1.038 
RHF 
631++g(d) 
1.157 
RHF 
631+g(d) 
1.166 
MP2 
augccpvdz 
1.354 
B3PW91 
augccpvtz 
1.442 
CCSD 
augccpvdz 
1.447 
B3LYP 
6311+g(3df,2p) 
1.458 
BB1K [31] 
MG3 [26] 
1.467 
MP2 
631+g(d) 
1.522 
B3LYP 
631+g(d) 
1.583 
PW91PW91 
augccpvtz 
1.615 
MP2 
ccpvtz 
1.823 
CAS(6,6) p 
augccpvdz 
2.193 
CAS(6,6) p 
631g 
3.662 
There are three trends worth noting for the results in Table 1. One: for any given method, the results generally improve with larger basis sets, and the inclusion of diffuse functions is particularly important, as noted by Ozkabak [3]. Two: HartreeFock seems to be superior to methods that include in various ways the electron correlation energy. Three: all of the methods tested overestimate the intensity of the CH stretching mode compared to the CC stretching mode. It is also clear from Table 1 that there is no benefit for predicting Raman intensities from treating the basis functions that contribute to the 6 p MOs in a completely equivalent fashion as was done in the CAS(6,6) calculations.
Most of the
methods in Table 1 were also applied to C_{6}D_{6} and C_{6}F_{6},
except for the CAS(6,6) calculations (since they were quite poor for C_{6}H_{6}),
and the CCSD calculations since they were too expensive for C_{6}F_{6}.
The results for these methods on two other species are shown in Table 2 below.
Table 2
Method 
Basis set 
I(n_{1})/I(n_{2}) C_{6}D_{6} 
I(n_{1})/I(n_{2}) C_{6}F_{6} 
Experimental 

0.646±0.004 
0.0178±0.0003 
RHF 
6311+g(3df,2p) 
0.662 
0.140 
RHF 
631+g(d) 
0.776 
0.134 
MP2 
augccpvdz 
0.860 
0.073 
PW91PW91 
augccpvtz 
0.908 
0.0317 
B3PW91 
augccpvtz 
0.908 
0.0604 
B3LYP 
6311+g(3df,2p) 
0.92 
0.0650 
BB1K [31] 
MG3 [26] 
0.943 
0.106 
MP2 
631+g(d) 
0.984 
0.0891 
B3LYP 
631+g(d) 
1.021 
0.0660 
MP2 
ccpvtz 
1.182 
0.0850 
The results in Table 2 show that, like the case of C_{6}H_{6}, all of the methods tested overestimate the CX stretching intensity compared to the CC stretching intensity. They also suggest that HatreeFock with a large basis set is successful for C_{6}D_{6}, as it was for C_{6}H_{6}, but that it is quite poor for C_{6}F_{6}. Table 2 also makes it clear that none of these methods is very good for C_{6}F_{6}. It seems that accurate prediction of the intensity of both the most intense, and one of the least intense bands in the Raman spectrum is quite challenging. The PW91PW91/augccpvtz calculations are much better than the others, overestimating the CF stretching intensity by only 80% compared to the CC stretching intensity; although rather poor for C_{6}H_{6}, this method did rather well on C_{6}D_{6} also.
Discussion: It is quite unexpected that the HartreeFock methods seem to be superior to methods that include the electron correlation energy for predicting Raman intensity in C_{6}H_{6}. It seems likely that this unexpected result is due to a fortuitous cancellation of errors. It is possible that this cancellation involves the anharmonic nature of the CH stretching modes. Since the automatic Raman activities in Gaussian are computed on the assumption of harmonic vibrations, and the derivation of equation 2 is also based on the assumption of both electrical and mechanical harmonicity, it may be that neglecting anharmonicity may cancel some of the error in neglecting electron correlation energy. Indeed, Maslen et al. [32] have shown that at the HartreeFock level, anharmonic corrections to Raman intensity are quite small, but they make no claim that this will be the case for correlated methods. We are investigating methods to develop an intensity formula like equation 2 that does not depend on the assumption of mechanical harmonicity.
With the methods
described here, it is possible to simulate the entire Raman spectrum from the
activities computed by Gaussian. Equation 2 may be used to convert these
activities into the corresponding intensities predicted for a given sample
temperature, laser frequency, and type of detector. A simple way to do this is
to assume that all of the Raman band will be described by Lorentzian peaks with
a common width, whose areas are proportional to the predicted intensity.
Figure 5 below shows such a simulation for the complete Raman spectrum of C_{6}D_{6},
using the results from a calculation at the RHF/6311+g(3df,2p) level of theory.
Figure 5
In this figure all of the predicted frequencies were scaled so that the
predicted and observed n_{1}
frequencies would match; the predicted intensities were scaled so that the n_{1 }intensities would match;
then a common peak width for all of the Lorentzians was chosen so that the n_{1} lineshapes would match as
closely as possible. The black curve is the experimental intensity, corrected
for spectrometer response; the blue curve is the predicted Raman intensity;
and the red curve is the predicted Raman activity, uncorrected for laser
frequency, temperature, and type of detector. The inset clearly shows that
such correction is quite important; the area under the blue curve is nearly
equal to that under the black curve, but they have different peak widths. The
red curve has significantly smaller area.
Conclusions: For predicting the Raman intensities for hydrocarbons we recommend the use of HartreeFock calculations with a large basis set including diffuse functions. For other types of molecules, or where CH stretching intensity is not of interest, perhaps a DFT method such as PW91PW91 with a basis set including diffuse functions would be better. This latter choice has the advantage that it will also produce more accurate vibrational frequencies than will the HartreeFock calculations.
Acknowledgements: A portion of the research described in this manuscript was performed at the W.R. Wiley Environmental Molecular Sciences Laboratory, a national scientific user facility sponsored by the U.S. Department of Energy’s Office of Biological and Environmental Research and located at the Pacific Northwest National Laboratory. PNNL is operated for the Department of Energy by Battelle. Computer time was supported by collaborative hpc resources funded through the University of North Carolina Office of the President.
References: