ARTICLE pubs.acs.org/JPCA
Vibrational Spectra of r-Glucose, β-Glucose, and Sucrose: Anharmonic Calculations and Experiment Brina Brauer,† Madeleine Pincu,‡ Victoria Buch,† Ilana Bar,§ John. P. Simons,|| and R. Benny Gerber*,†,‡ †
Institute of Chemistry and The Fritz Haber Research Center, The Hebrew University of Jerusalem, Jerusalem, 91904, Israel Department of Chemistry, University of California, Irvine, California 92697, United States § Department of Physics, Ben-Gurion University of the Negev, Beer Sheva 84105, Israel Chemistry Department, Physical and Theoretical Chemistry Laboratory, Oxford University, South Parks Road, Oxford, OX1 3QZ, U.K.
)
‡
ABSTRACT: The anharmonic vibrational spectra of R-D-glucose, β-D-glucose, and sucrose are computed by the vibrational self-consistent field (VSCF) method, using potential energy surfaces from electronic structure theory, for the lowest energy conformers that correspond to the gas phase and to the crystalline phase, respectively. The results are compared with ultraviolet-infrared (UV-IR) spectra of phenyl β-Dglucopyranoside in a molecular beam, with literature results for sugars in matrices and with new experimental data for the crystalline state. Car-Parrinello dynamics simulations are also used to study temperature effects on the spectra of R-D-glucose and β-D-glucose and to predict their vibrational spectra at 50, 150, and 300 K. The effects of temperature on the spectral features are analyzed and compared with results of the VSCF calculations conducted at 0 K. The main results include: (i) new potential surfaces, constructed from Hartree-Fock, adjusted to fit harmonic frequencies from Møller-Plesset (MP2) calculations, that give very good agreement with gas phase, matrix, and solid state spectra; (ii) computed infrared spectra of the crystalline solid of R-glucose, which are substantially improved by including mimic groups that represent the effect of the solid environment on the sugar; and (iii) identification of a small number of combination-mode transitions, which are predicted to be strong enough for experimental observation. The results are used to assess the role of anharmonic effects in the spectra of the sugars in isolation and in the solid state and to discuss the spectroscopic accuracy of potentials from different electronic structure methods.
I. INTRODUCTION Vibrational spectroscopy is a very powerful tool for exploring the properties of molecular systems, including their underlying potential energy surfaces. Essential for this purpose is the combination of experiments of high resolution with accurate theoretical and computational tools of interpretation. In this paper, our aim is to present state-of-the-art calculations for three representative saccharides: R-D-glucose, β-D-glucose, and sucrose. The calculations are compared with literature data from other laboratories and also with the results of new experiments carried out in the present study. Anharmonic interactions in the potential energy surfaces of the sugars and their effects on the spectra are focal points of interest. Another aim is to explore anharmonic interactions both for the isolated sugar molecules, in a noninteracting (gas phase) or in a weakly interacting environment (solid matrices), and for saccharides in crystalline phases. The conformations adopted by the sugar molecules in the gas phase and in the crystal are different, and it is of interest to explore the differences between the anharmonic effects in the two environments. A third objective is to examine potentials derived from electronic structure calculations, with regard to their spectroscopic accuracy, and to their description of anharmonic effects. Infrared (IR) studies in the crystalline phase and in other solid phases are among the earliest approaches to the spectroscopy of saccharides.1-6 The spectra of ordered, crystalline samples are often of good resolution for several of the bands, but others can be influenced by intermolecular interactions, involving hydrogen r 2011 American Chemical Society
bonds with neighbors in the solids, which are typically quite strong, hence the lowest energy conformer in the crystal is not the same as in the gas phase. On the other hand, sugars in a raregas matrix only experience a weakly interacting environment. The vibrational frequency shifts are typically quite modest, and the sugar conformers in matrices are expected to correspond basically to those populated in the gas phase at low temperatures. Unfortunately, relatively few spectroscopic studies of sugars in matrices have been reported so far.7 The highest resolution spectra of saccharides (and of other biomolecules) are provided by UV-IR double resonance vibrational spectroscopy, typically recorded in the hydride stretch region through “ion dip” measurements in molecular beams.8-11 This enables exploration of the molecules in nearly ideal isolation at low temperatures, and it has been applied to a wide range of carbohydrates.8-11 The method requires the introduction of a UV chromophore, typically a phenyl or benzyl group, which is not present in most “natural” sugars, but in almost all cases the UV “tag” does not significantly perturb the conformational structure.11 In carrying out vibrational spectroscopy calculations, the description of anharmonic effects is a major challenge. Harmonic Special Issue: A: Victoria Buch Memorial Received: October 20, 2010 Revised: February 2, 2011 Published: March 10, 2011 5859
dx.doi.org/10.1021/jp110043k | J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A calculations are computationally feasible, even for relatively large polyatomics, but the accuracy of the harmonic approximation is insufficient for many transitions. Anharmonic effects have on the whole greater effects on vibrational spectra than on thermodynamic properties.12,13 In many applications, including the vibrational spectroscopy of sugars, empirically based scaling factors14 are applied to the calculated harmonic frequencies to improve agreement with experiment. Such an empirical approach cannot, of course, shed light on anharmonic contributions to the potential energy surface, which is a topic of much interest. In the present study, we employ the vibrational self-consistent field (VSCF) approach15-19 to compute the anharmonic vibrational spectra of glucose and sucrose, from potential surfaces of electronic structure methods. Variants of the VSCF method developed in recent years by the Gerber group19-21 are used for this purpose. These VSCF algorithms have been employed in recent years in spectroscopic calculations for a range of biological molecules, and generally yield results in very good accord with experiment.22-29 The VSCF variants are realistically applicable to molecules of the sizes targeted here and beyond, and can be used directly with potentials from electronic structure codes, generated as points on suitable multidimensional grids. The only anharmonic calculation of sugar spectra previously reported in the literature is by Gregurick et al.30 That calculation does indeed point to substantial anharmonic effects in part of the IR spectrum of glucose but since it used an Assisted Model Building with Energy Refinement (AMBER) force field, it seems important to employ more reliable potentials from ab initio theory. To assess the quality of the potential energy surfaces of the saccharides, the calculations are compared with new experimental results obtained in the present study for crystalline glucose anomers, where the environment of the sugar molecules very significantly affects the spectra and must be included in the modeling, and also with published experimental spectra recorded in weakly interacting or noninteracting environments (matrices, molecular beams).31,32 In addition to VSCF calculations for sugar molecules at 0 K, we also present calculations of spectra at finite temperatures. These are obtained in the classical approximation from molecular dynamics (MD) simulations in the Car-Parrinello (CP) approach.33-35 Classical-limit spectroscopic calculations have been carried out previously, for example, for peptides and other biomolecules by Gaigeot,36,37 and for other systems by Buch.38,39 The CP calculations provide insight into the effect of temperature on sugar spectroscopy. The structure of the paper is as follows: The theoretical and the experimental methods used here are described in section II. Results and the analysis of comparison of calculations with experiments, are presented in section III, and concluding remarks are presented in section IV.
ARTICLE
standard VSCF codes, separability is assumed for normal coordinates.17-19 This provides a substantial improvement over the harmonic approximation, but better accuracy requires the incorporation of effects beyond vibrational separability. An approach that has proved satisfactory in accuracy in a range of applications, and that remains computationally feasible for relatively large systems, first proposed by Norris et al.17 and by Jung and Gerber,18 includes the corrections for separable VSCF by second-order perturbation theory. Several different terms are used in the literature for this approach, including correlation-corrected VSCF (CC-VSCF),18,19 or vibrational MP2 theory.17 Here we will apply the term VSCF with secondorder perturbation theory corrections (VSCF-PT2). While VSCF-PT2 has proved to be of good accuracy in many applications, this algorithm is considerably more demanding computationally than separable VSCF. Pele et al.21 were able to reduce the computational effort of VSCF-PT2 by using algebraic properties to simplify the equations involved. This reduces the computational effort for large N by as much as N3, where N is the number of vibrational modes. For the systems studied in the present paper, in particular for sucrose, the simplified algorithm of Pele et al. is essential; it is implemented in the VSCF codes of the GAMESS suit of programs,42,43 and is employed in the present calculations. Several other “accelerations” of the VSCF and the VSCF-PT2 algorithms were recently proposed in the literature, and test applications using them seem encouraging.44-46 Vibrational degeneracy effects may, in principle, influence some of the transitions calculated in this study. On the basis of experience with nucleic acid bases,26 amino acids, and several other biomolecules,25,28,29 we assume here that the effects due to vibrational degeneracies are relatively small in moderate to large biomolecules. Treatments of vibrational degeneracies in the VSCF approach were introduced by Matsunaga et al.,47 by Yagi et al.,48 by Danecek and Bour49 and by Respondek and Benoit.50 However, the implementation of these algorithms for large molecules is, in the present state of the art, computationally very demanding. A calculation that incorporates vibrational degeneracy effects was carried out for glucose, in the present study, but the results provided did not improve agreement with experiment compared with the nondegenerate VSCF-PT2 treatment. For other saccharides studied here, vibrational degenerate effects were not treated. Approximate representations of the potential energy surface need to be used in VSCF calculations, in order to reduce the dimensionality of integrals in the algorithm, and bring the calculations to a feasible level of effort. Two types of potential representations proved useful: one represents the potential surface approximately, in terms of a sum of single mode potentials, and a sum of pairwise mode-mode interactions:18,19 N
V ðQ 1 , :::, Q N Þ ¼ V0 þ
II. METHODS AND MODEL SYSTEMS A. Theoretical Methodology. 1. VSCF Algorithms. The VSCF working equations and algorithms used in this work were previously described in several articles;19-21 here we briefly summarize several points of importance: The simplest level of the VSCF approximation assumes separability of the vibrational wave functions in a chosen system of coordinates40,41 and in the great majority of applications, and in the
N N
∑ VjðQ jÞ þ ∑i j∑> i Vij ðQ i , Q jÞ j¼1
ð1Þ
where Q1,...,QN are the normal modes. This representation can be conveniently generated both for empirical force fields18 and for ab initio potentials.19 While the representation is approximate and must be tested on a case by case basis, it has been applied extensively, including several VSCF calculations of biological molecules,20,22-29 and found to be very satisfactory in the great majority of cases. 5860
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Another approach is to expand the potential in a power series in the normal-mode displacements:51,52 V ðQ 1 , :::, Q N Þ ¼
∑n ∑n ::: ∑n An :::n ðQ 1 Þn :::ðQ N Þn 1
1
1
2
N
N
ð2Þ
N
Since the expansion has to be truncated after a limited number of terms, this representation is also inherently approximate. The quartic force field (QFF), due to Yagi et al.,52 is directly implemented in the VSCF code of GAMESS. In actual applications of QFF, only a limited number of quartic terms in eq 2 are employed, and it is likely that this representation is less accurate than eq 1. QFF is computationally more efficient, however, and although most of the VSCF calculations in this work were done with the pairwise coupling representation, several employed QFF. 2. Adjusted Hartree-Fock for Anharmonic Interactions. Electronic structure methods at the level of MP2 or of the hybrid density functional B3LYP have been found to give excellent agreement with experimental frequencies when used in VSCFPT2 calculations for small biological molecules22,24,53,54 such as glycine and N-methyl acetamide (NMA). The calculation of the anharmonic potential surfaces for larger biological molecules with such methods is computationally very expensive, however, and in the present study, MP2 calculations of the anharmonic potential energy surfaces were carried out for glucose, but not for sucrose. In this paper, we follow an approach proposed by Brauer et al.20 for “upgrading” potential surfaces of a relatively low level, to produce more accurate anharmonic surfaces. Most of the “upgraded” potentials constructed here are obtained from the Hartree-Fock (HF) method, which is computationally much faster than MP2. Let V HF(Q 1,...,Q N) denote the HF potential surface in the normal-mode coordinates Q 1,...,Q N. Following ref 20, the adjusted HF potential (AHF) is defined by a scaling of the coordinates in VHF: V AHF ðQ 1, :::, Q N Þ ¼ V HF ðλ1 Q 1 , :::, λN Q N Þ
ð3Þ
where the scaling coefficients λj are defined in terms of the harmonic frequencies of the “low level” potential (here HF) and of a selected “high level” potential. Here, we use MP2 as the “higher level” potential, employed in adjusting HF: =ωHF λj ¼ ωMP2 j j
ð4Þ
where ωMP2 ,ωHF are the j-mode harmonic frequencies computed j j by MP2 and HF, respectively. With this scaling, the AHF potential has the same harmonic frequencies as MP2, which are considered very accurate. The scaling of eq 3 seems reasonable provided the normal modes of the “higher level” (MP2) and of the “lower level” (HF) method are very similar. This was verified in the calculation here. We note that in previous VSCF studies of a range of other molecules, the semiempirical PM3 was “upgraded” and found to be successful in spectroscopic calculations;20,26-28 the PM3 normal modes for sugars, however, were found to differ significantly from the MP2 ones, although not for OH stretches. Thus, only for some of the systems and for several frequencies was PM3 used. The main potential surface employed was AHF. 3. “On-the-Fly” Dynamics Computations of Vibrational Spectra. Vibrational spectra can be computed from MD simulations based on Fourier transformations of system dipole moment correlation functions obtained along classical trajectories.55,56 In the calculations reported here, we used “on-the-fly” dynamics
simulations with the CP2K/QUICKSTEP package,33,57 the CP implementation of the density functional theory (DFT). Computations were carried out with ab initio potentials at the BLYP level of the Kohn-Sham formulation of DFT and with the double-ζ valence polarization (DZVP) basis set, while the pseudopotential was the Goedeker, Teter, and Hutter (GTH)type pseudo potential.58 We note that in the CP2K code as applied in this work, the one-electron orbitals were expanded in a plane-wave basis set, with a kinetic energy cutoff of 280 Ry, using the default values for convergence from the paper of Martyna and Tuckerman.59 Though expected to be less accurate at low temperatures than quantum mechanical methods such as VSCF, this quasi-classical method has several advantages. First, we note efficiency and scalability benefits that make it practical to use for relatively large systems that cannot otherwise be explored with ab initio methods such as MP2. Most importantly for the present study, it has the merit of enabling calculations of vibrational spectra at finite temperatures.36,55,60 The resulting spectral lines reflect, therefore, thermal broadening and can be directly compared to experimental data at corresponding temperatures, provided the computed trajectories explore all the conformers that contribute to the spectrum at the temperature of interest. This can be accomplished by a sufficiently long trajectory or by sampling a set of shorter trajectories from different initial conditions. Application of this method to calculations of peptide spectra at room temperature has been found to have reasonable accuracy36 despite the fact that potentials based upon the BLYP functional are not considered to be as reliable as B3LYP ones. The method has also been applied to the study of water clusters and its surface acidity,60 and several observed properties were correctly reproduced. 3.1. Computational Details. Ab initio MD trajectories of 4 ps duration were run on constant energy ensembles that consisted of R-D-glucose and β-D-glucose in the gas phase. The trajectories of nuclei moving on the Born-Oppenheimer potential energy surface were calculated at time step intervals of 1 fs. System dipole moments were computed at each time step over a cubic unit cell of side dimension 16 Å. The one-electron orbitals were expanded in a plane wave basis set with a kinetic energy cutoff of 280 Ry. Simulations were carried out at average temperatures of ∼300 and ∼150 and ∼50 K. It is to be noted that while the time step interval and the average temperature of the ensemble determine the IR spectral resolution, the duration of the simulation (td) determines the lowest wavenumber that can be accurately predicted by the MD simulation in the following way: ω > πn/td, where n is an integer (n > 2) and td is the trajectory duration. In our case, we estimate that the lowest predictable frequency is ∼315 cm-1. Assignment of spectral bands to individual atomic modes of motion was performed by matching the Fourier transform of the time dependence of either the bond length, the bond angle, or the torsion angle formed by selected groups of atoms, to corresponding IR bands. B. Experiment. 1. Measurement of IR and Raman Spectra of Solids. R-D-Glucose, β-D-glucose, and sucrose were purchased from Aldrich and were used without further purification. The IR absorption spectra of these samples were monitored, using a Nicolet Impact 410 Fourier transform IR (FTIR) spectrometer, operated under the Omnic IR data software. The spectra of transparent KBr pellets, prepared by grinding 1 mg of crystalline particles of each sample with 100 mg of dry KBr, were recorded in 5861
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A the 400-4000 cm-1 range at spectral resolution of 2 cm-1 with 2000 accumulations. The Raman spectra were collected with a Jobin-Yvon, LabRam UV HR micro-Raman spectrometer. A diode laser at 784.9 nm (∼10 mW intensity) was used for excitation and was focused onto each sample with an 50/0.75 numerical-aperture microscope objective. The scattered light was redirected from the microscope through a sharp edge long wave-pass filter that rejected the excitation laser line and the elastically scattered light, and through a confocal pinhole for increased axial resolution. The scattered light was then focused into a 0.8 m dispersive spectrometer, equipped with a 600 groove/mm grating, and the scattered Raman light was detected with an air cooled charge-coupled device (CCD), consisting of 1024 256 pixels. The spectra in the 100 - 3300 cm-1 region were monitored while scanning the spectrum across the CCD and moving the grating. Successive recordings of Raman spectra from the samples allowed monitoring of their characteristic spectral features. The definition of the measurement parameters and measuring control was done by the LabSpec 4.04 software. The Raman spectra contained a varying background overlaying the signal, and this background was removed by fitting it to a high-order polynomial and subsequently subtracting this polynomial from the spectra. 2. Spectroscopy in the Gas Phase. Detailed descriptions of the experimental strategy have been published previously.11 The structures of individual conformers, transferred through oven or laser vaporization into an expanding Ar jet and isolated in a cold molecular beam (T ∼ 10 K), were probed through a combination of ultraviolet (UV) resonant two photon ionization, and infrared “ion-dip” (IRID) vibrational spectroscopy, conducted in the near and mid-IR regions. IRID is a double resonance, IR-UV, spectroscopic technique that depletes the ground state population of the target molecule whenever the IR radiation is “in tune” with one of its fundamental vibrational modes. Monitoring the resulting dip in the resonant two photon ionization signal thereby provides a means of recording its vibrational spectrum; in the present system, the UV chromophore was provided by the phenyl tag.
III. RESULTS AND DISCUSSION A. Spectroscopy of Glucose in Ar Matrix. We begin the presentation of our results by comparing our VSCF calculations for R-D-glucose with the matrix spectroscopy experiments of Kovacs and Ivanov.7 In this case, comparison between the computational results and the experiment covers essentially the entire fragmentary range. We assume that the matrix environment does not affect the structure of the saccharide, which seems very reasonable in view of the weak guest/host interaction in this case. Further, we assume in the analysis that matrix shifts are negligible, or at least smaller than the estimated error of the calculations. This is supported by previous spectroscopic analysis of other small biomolecules in rare-gas matrices, e.g., glycine,22 NMA,53 alanine and proline,20 although the intensities may be affected.61 Figure 1 shows the structures of the three lowest energy conformers of isolated R-D-glucose, obtained in MP2/TZP calculations in the present work. The conformers are designated by symbols G and T for gauche and trans orientations of the CH2OH group as used by Talbot and Simons62 and Corchado et al.32 The plus and minus designations indicate the signs of relevant torsional angles; cc indicates a “counter-clockwise” orientation of the OH groups on the pyranose ring. The energy differences between the three structures are extremely
ARTICLE
Figure 1. Lowest energy conformers of R-D-glucose, with relative MP2/TZP energies in kcal/mol. See also ref 32.
Figure 2. β-D-Glucose conformers, with relative MP2/TZP energies in kcal/mol, as in Figure 1.
small. The energy ranking of the three conformers in our calculations differ from the B3LYP (DFT) results of Corchado et al.,32 although structures of corresponding conformers in the two calculations are very similar. MP2 is generally expected to be somewhat more accurate than DFT for organic molecules, but the energy differences are so slight that no conclusion can be made as to the true global minimum. A similar situation is found for β-D-glucose, where the calculations also predict several structures with very small differences, as shown in Figure 2. While the energy rankings are questionable for the lowest energy conformers, there are some potentially significant differences between the structures, involving for example, hydrogen bonded interactions and those associated with the orientation of the CH2OH side chain. The similarity of their energies most likely reflects a balance between hydrogen bonding strength and steric effects associated with the CH2 groups of CH2OH. It seems reasonable to assume that in the matrix experiment, all three conformers, and perhaps additional low-energy ones, are trapped with similar probability. On the basis of this assumption, VSCF-PT2 calculations were performed for each of the conformers shown in Figure 1. The potential energy surface used in all cases was obtained using the hybrid HF/MP2 method described in section II (adjusted to reproduce MP2/TZP frequencies at the harmonic level of approximation). The results were compared with the experimental data of Kovacs and Ivanov.7 The comparison of theoretical and experimental frequencies and intensities, after an assignment in terms of the three conformers of Figure 1 is shown in Appendix 1. Close agreement between the computed and experimental frequencies is seen throughout the range of the data. There are a few exceptions and transitions for which the theoretical assignment is uncertain. These could be due to the presence of additional conformers, or other factors not considered in the calculations, for example, contributions from combination mode transitions or vibrational degeneracy effects. However, for many of the transitions, agreement between the computed and experimental frequencies is closer than the estimated error bound for the computational methods employed. (We have found typical frequency errors of 1-2%, due to the 5862
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Figure 3. (Left) OH vibrational “ion dip” spectra associated with the three lowest energy conformers of phenyl β-glucose together with their structures, calculated using DFT (B3LYP/6-31þG*). (Right) corresponding mid-IR spectra; black: experimental; red: calculated (both shown as positive going traces).11.
inaccuracies of the VSCF-PT2 method, and of the potential used). The agreement seems to hold nicely down to the lowest frequency modes, which are collective vibration and torsions. The comparison for the intensities is much less satisfactory. Although there is good quantitative agreement for quite a few transitions, theory and experiment differ by a factor of 2 or so in many cases, and there are several larger discrepancies. We note, for example, that computed intensities are, in cases, quite sensitive to the above-mentioned scaling procedure. Keeping in mind the difficulties of accurate intensity calculations even for small molecules, and possible matrix effects and differing conformer populations, however, the computational results are in good accord with much of the experimental data. The main conclusion that can be drawn, especially from the close agreement of the frequencies, is that the potential surface employed is quite accurate for this system, including of course, the anharmonic interactions. That the agreement holds over the entire frequency range shows that the potential energy surface is correct along each of the normal modes. This is very encouraging for the hybrid HF/MP2 potentials introduced here, and suggests the desirability of using them also for nonspectroscopic properties of glucose, and perhaps for other saccharides. B. Spectroscopy of Phenyl β-Glucose in a Molecular Beam. Spectroscopy in molecular beams provides essentially perfect isolation for the species of interest. These experiments have therefore special conceptual importance for testing calculations and determining the limits of the computational accuracy. In this section, we compare results on the vibrational “ion dip” spectroscopy of phenyl β-D-glucose, previously published11 and reproduced in Figure 3, with scaled harmonic calculations conducted in the OH stretch region.11 The figure also shows the three lowest energy conformers.11 The main difference between the three structures is
Table 1. Comparison of VSCF-PT2 Frequencies (in cm-1) of β-Gþg Glucose with and without Phenyl: PM3 Potentials Are Used mode
frequency/cm-1 without phenyl frequency/cm-1 with phenyl
OH str 2
3661
3660
OH str 3 OH str 4
3655 3653
3682 3657
OH str 6
3663
3644
Table 2. Comparison of Computed Anharmonic OH Stretching Frequencies (in cm-1) of (a) Phenyl Substituted βD-glucose (Gþg) and (b) β-D-glucose (Gþg), Using HF/MP2 and PM3 Potentials mode
frequency/cm-1 HF/MP2
frequency/cm-1 PM3
OH str 2
3604
3660
OH str 3
3619
3682
OH str 4
3615
3657
OH str 6
3612
3644
mode
frequency/cm-1 HF/MP2
frequency/cm-1 PM3
OH str 2
3683
3661
OH str 3 OH str 4
3679 3667
3655 3653
OH str 6
3688
3663
the configuration of the CH2OH side group with respect to the glucose ring. VSCF-PT2 calculations for the OH stretching frequencies of the β-Gþg conformer (Figure 3) are shown in Table 1 where 5863
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Table 3. Comparison of Phenyl-β-D-glucose Calculated Frequencies (cm-1) with Experimental Data for the Phenyl-β-D-glucose Conformersa calculated “A” Gþg-
experiment conf. B
calculated “B” G-gþ
experiment conf. X
calculated “X” T
3640.5
3660
3634
3640
3627
3637
3645 3634
3682 3657
3645 3636
3671 3645
3642 3642
3655 3649
3632
3644
3640
3665
3604
3560
mode
experiment conf. A
OH str 2 OH str 3 OH str 4 OH str 6
VSCF-PT2 calculations (using PM3 uncorrected by MP2 potentials) for β-D-glucose conformers Gþg-, G-gþ, and T are compared with experimental data for conformers A, B, and X of phenyl β-D-glucose, respectively. a
they can be compared with the corresponding frequencies for (nonsubstituted) β-D-glucose with “OH str 2” indicating that the OH stretch is located primarily on carbon 2 according to the standard numbering system. Although the calculations used the PM3 electronic structure method, which is computationally very fast but generally of limited accuracy, it does seem to give reasonable results for the OH stretching frequencies. The computed frequencies deviate by less than 1% from the experimental value. Nevertheless, the accuracy of the calculations is insufficient for predicting the frequency ordering of the transitions. Note that the frequency ordering of the phenylated and nonphenylated systems is not the same. The close correspondence between the frequencies of substituted and nonsubstituted glucose reconfirms the previous finding11 that phenyl substitution does not significantly affect OH stretching spectra in glucose and in numerous other saccharides. To assess the accuracy of the PM3 potential surface, we compare, in Table 2, the VSCF-PT2 results for this electronic structure method with those determined for the hybrid HF/MP2 potentials described in section III. (In the HF/MP2 method, the HF potential surfaces, which are adjusted to fit MP2/TZP frequencies at the harmonic level, are expected to be much more accurate.) The agreement is very encouraging with regard to the computationally very fast PM3, although it is slightly better for β-glucose than for phenyl substituted β-glucose. It shows that this method can be used for the OH stretching spectra of saccharides. This is also in agreement with previous applications for other systems, e.g., amino-acids and nucleobases,20,25-29 which showed that PM3/MP2 or PM3/B3LYP hybrid potentials give very good results for anharmonic spectra. Finally, in Table 3, we compare the computed anharmonic frequencies with the experimental data11 for phenyl β-D-glucose, employing the computational results (VSCF-PT2) for the substituted β-D-glucose conformers using PM3 surfaces (Gþg-, Ggþ, and T), since accurate HF/MP2 hybrid potential surfaces are only available for one conformer of the phenyl-substituted β-Dglucose. Two observations can be made. The computed frequencies are in very good agreement with corresponding OH frequencies observed for the three conformers, with deviations from measurements ranging from 6 cm-1 in the best case to 44 cm-1 in the worst, and on average 19 cm-1; the majority of the measured samples deviate by 1% or less from computed frequencies. The accuracy of the frequencies is encouraging not only for the potential surface, but for the VSCF-PT2 method, which yields the anharmonic frequencies on a first-principles basis. Table 3, however also shows the limitation of the computational results at this stage. The computed frequencies cannot be used to identify the three conformers with any certainty. Despite the general agreement between experiment and theory in Table 3, there seems to be a substantial discrepancy with regard to the OH str 6 frequency, for which the computed frequency of the X
structure is too low. We attribute this to the fact that, in this case, the low-level PM3 potential was used (as opposed to the higher quality AHF potentials, in other parts of this work). The PM3 potential has well-recognized difficulties, although it was used here because of its very low computational cost. Identification of a conformer on the basis of a single frequency, or even all the OH stretching frequencies, is not feasible in the present state of the art of computational accuracy. For a more complete theoretical assignment of the experimental data, comparison for a broader range of observed frequencies is essential. Indeed, in the mid-IR data for phenyl-β-D-glucose shown in Figure 3, there does seem to be a clear distinction between the spectra, both experimental and DFT calculated, of the Gþ, G- and T conformers (A, B, and X; compare the spectra between ∼1000 and 1100 cm-1). It is not certain at this point if this will suffice, or if the entire spectrum, down to low frequencies will be required to determine the contributions of the three observed conformers. Availability of quantitative intensities should also be very helpful. Availability of data for the full vibrational spectrum, and of measured intensities make possible the complete analysis of the spectrum in the matrix, as discussed in section (II.A) above. C. Spectroscopy of Glucose and of Sucrose in the Crystalline Phase. In crystalline glucose, interactions with neighboring molecules have a major effect on the conformation of the sugar. The structures of R-D-glucose and of β-D-glucose in isolation (gas phase or matrix) involve internal hydrogen bonds. In the crystal, the OH groups instead form hydrogen bonds with neighboring molecules. For calculations of the spectra, it is essential to introduce models that represent this effect in a computationally feasible framework. Therefore, we employed a mimic system for the solid-state calculations for glucose. The geometry of a central glucose molecule with its nearest hydrogen-bonding neighbors was taken from the Cambridge Crystallographic Database.63,65 Nine neighboring glucose molecules that are hydrogen bonded by OH groups were replaced by methanol. One molecule that is hydrogen bonded through its ring oxygen was replaced by dimethyl ether. This set of molecules was structurally optimized by the semiempirical PM3 method. As it turns out, the crystal structure shows that one of the molecules is farther than the typical hydrogen bonding distance, so the methanol was excluded from the calculation. The use of CH3OH molecules as mimic groups greatly simplifies the structure optimization, and obviously also the spectroscopy calculations discussed below. The structure of the mimic used in the spectroscopic calculations is shown in Figure 4. The IR spectrum of the system was computed using PM3 potentials and the VSCF-PT2 vibrational algorithm, assuming 0 K temperature. This produces a spectrum of sharp lines that were smoothed by convolution with a Lorentzian of width of 10 cm-1, to model the broadening effects due to the temperature and the condensed phase environment in the real experimental system. 5864
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Figure 6. Structures of sucrose: (a) in crystal; (b) crystalline-like conformer, optimized in MP2; (c) in gas phase.
Figure 4. Mimic used for spectroscopy of crystalline R-D-glucose.
Figure 7. Raman spectrum of crystalline sucrose. VSCF-PT2 calculations for the hybrid HF/MP2 potential are compared with experiment.
Table 4. Comparison of OH Stretch Modes in r-D-Glucose and β-D-Glucose (Gþg), Obtained by VSCF-PT2 and DFT (BLYP) frequency/cm-1
R-D-glucose
Figure 5. IR spectra of R-D-glucose in crystal calculations, with and without mimic are compared with experiment.
Transitions due to the mimic group were not included in the spectrum. We expect the OH groups in the glucose crystal to be very strongly coupled to other modes, as the geometry suggests such strong interaction. Our vibrational method is inadequate for treating such strongly coupled vibrations, hence we do not report here the results for the OH stretching manifold. (For very strong coupling to the environment, effects such as fast vibrational energy flow may greatly affect the spectrum.) The results of the calculations are shown in Figure 5. Given the challenge of quantitative modeling of crystalline spectra, we consider the agreement to be good. The main discrepancy between theory and experiment is in the region of the CH stretching spectra. We consider it likely that this is due to structural differences between the crystal and the mimic (or the isolated R-D-glucose) that substantially influence the CH2 and CH groups of the sugar. The broadening effects in the CH stretching part of the spectrum are certainly far stronger in the
β-D-glucose
-1
mode
frequency/cm VSCF-PT2
OH1
3644
3700
OH2
3629
3630
1
OH3
3634
3670
36
OH4
3642
3670
28
OH6
3668
3670
2
OH1 OH2
3658 3683
3690 3690
32 7
OH3
3679
3680
1
OH4
3667
3660
-7
OH6
3688
3670
-18
DFT trajectory, (main peak)
frequency difference/cm-1 56
real system than our model can describe. It should be kept in mind that the treatment of broadening in our simple model is by the introduction of an ad-hoc width. The mechanisms of broadening are of course not obtainable in such a framework. There is a good degree of agreement between the computed spectra and experiment in the range of 500-1500 cm-1, this being especially the case for the results of the mimic. The main peak positions, and also to an extent the relative intensities, are well produced in this image. The calculations for the R-D-glucose without the mimic groups shown in Figure 5 are for the conformation that corresponds to the crystal, which differs considerably from the gas-phase structure, but no other effect of the environment is 5865
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Figure 8. IR spectrum obtained from system dipole by DFT/BLYP at a temperature of ∼50 K, with expanded CH and OH stretch modes. R-D-glucose - black line; β-D-glucose - red line.
introduced. It seems that the mimic groups shift peaks especially those in the low frequency range ω < 1000 cm-1, bringing them into close agreement with experiment. This region is mostly due to torsions, shearing modes, and several collective stretching modes of the ring. We conclude that modeling of environment effects such as using the conformer for the crystal, in place of the mimics, and treating collective motions are important, at least for glucose. Consider now sucrose in the crystalline phase. The initial structure for the calculations was taken from the Cambridge Crystallographic Database.63,65 Starting from this, a conformer search was carried out, using Merck molecular force field (MMFF) potentials.66,67 Energy rankings were pursued by the ab initio restricted HF (RHF) method. Finally, the lowest energy structure was optimized by MP2/6-31(d)G. The crystallographic and the optimized structures are shown in Figure 6, which also shows the lowest energy conformer of sucrose in the gas phase. The MP2 structure is quite close to the X-ray structure in the crystal. Since the influence of the crystal environment on the structure of sucrose is less strong than for glucose, the spectroscopic calculations were conducted using the “crystal-like” MP2 optimized conformer and the hybrid HF/MP2 potential surface, as described in Sec. II. The frequencies were computed from this potential by VSCF-PT2, and the intensities are from the MP2/631(d)G calculation taking the laser frequency into account.68 The calculated and experimental Raman spectra of crystalline sucrose are compared in Figure 7. The computed line shape shows several peaks corresponding to the CH bands. The effect of smoothing, due the convolution with the line shape, simplifies the structure of the spectrum. Several significant discrepancies are seen, e.g., the twin peaks in the CH stretching region in the calculations, which are not seen in the experiment. The experiment does show a structure of several fine peaks, not fully resolved, which are not seen in the computed spectrum. These could be due to contributions from combination modes or due to degeneracy effects not included in the VSCF-PT2 calculations. On the whole, however, there seems to be a good correspondence between the experimental and the computed combined peak positions and also their relative intensities. D. Sugar Spectroscopy Computed from Classical MD and Temperature Effects on the Spectra. To explore the accuracy range and the physical insights provided by classical dynamics, we carried out calculations of frequencies and relative intensities for R-D-glucose and β-D-glucose using
the “on the fly” CP2K software package.57 The calculation employed the BLYP/GTH variant of DFT.33,58 Vibrational frequencies of the lowest DFT energy conformers of R-D-glucose and β-D-glucose, obtained with MD at finite temperature ∼50 K, using DFT (BLYP), were compared with VSCF-PT2 calculations. The comparison carried out for the OH modes is shown in Table 4 and the spectrum obtained is shown in Figure 8, focusing on both the OH and CH stretches. In discussing the assignments, we note that the calculations are done in normal modes, which are collective. The OH modes are thus not purely local, but involve some small measure of mixing with other atomic displacements. The deviations between the VSCF-PT2 and MD (CP2K) frequencies for the OH stretches range up to 56 cm-1, but are mostly smaller. We attribute the differences to the more accurate HF/MP2 potential used in VSCF-PT2, rather than to the different vibrational methods used in the two cases. Our conclusions on the spectroscopic accuracy of BLYP are therefore less favorable than the findings by Gaigeot and Sprik for other biological molecules.55 Nevertheless, the comparison of Table 4 suggests that BLYP is a reasonable approximation for frequencies of saccharides.54 While relatively inaccurate for low temperature conditions, the spectroscopic calculations by CP2K/DFT have important advantages over VSCF. Above all, the method is applicable at higher temperatures, where the system probes more than one minimum of the potential. In view of this, the method can also be used to predict the temperature effects upon the spectrum. Figure 9 compares the high-range of the spectra obtained for R-D-glucose and β-D-glucose at two temperatures: ∼50 and ∼300 K. Several observations can be made. First, the spectral bands appear broader in the higher temperature spectrum (top of the figure), which is to be expected. Second, the intensities of the spectrum at the higher temperature are nearly 1 order of magnitude higher than at low temperature. We also note that the contribution of β-D-glucose to the spectrum in this frequency region is stronger than that for R-D-glucose and has different characteristics. To delve deeper into the details of these differences, we next examine individual bands in the OH range, at three temperatures, for the R-D-glucose and β-D-glucose conformers. Figure 10 shows the comparison for the OH bands at ∼50 (bottom), ∼150 (middle), and ∼300 K (top). We first note that, overall, β-D-glucose covers a tighter frequency range in the OH band than R-D-glucose, regardless of the temperature, with OH4 being the band that 5866
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Figure 9. Temperature dependence of the IR spectrum CP2K/DFT (BLYP): (a) R-D-glucose at ∼50 (black) and at ∼300 K (red); (b) β-D-glucose at ∼50 and at ∼300 K.
is most affected (shifted to the red) by the doubling of the temperature. On the other hand, the OH bands for R-D-glucose are quite spread even at 150 K (except for OH3 and OH4 bands that overlap at 50 K), but temperature doubling broadens the bands and alters their shapes. The temperature effects of Figures 9 and 10 can be interpreted as follows: At T = 50 K, the thermal energy is low compared with the internal hydrogen bonds for R-D-glucose. The motions are restricted to a small vicinity of the global minimum, and the OH spectrum is strongly affected by the H-bonds that cause a large shift to the red from free OH. At 150 K, the thermal energy makes possible larger amplitude motions, which decrease the strength of the H-bonds and the red-shift effects. At 300 K, this effect is more pronounced, and in fact the system explores different conformations. In discussing differences between R- and β-glucose, we note the following: At equilibrium, R-glucose has slightly shorter H-bonds that β-glucose, indicating stronger bonds.69 However, in the course of the dynamics, R-glucose goes through more rotameric changes, which lead to larger variations in the H-bond length for this anomer. As a result, the β-glucose peaks cover a narrower range of frequencies (Figure 10). As T is increased with more rotameric motions, the H-bonds are weakened, and the OH transitions are shifted to the blue for both anomers. While frequencies computed from classical dynamics are inevitably approximate, trajectory calculations include two anharmonic effects that are not treated by VSCF. The first of these is vibrational energy flow between different modes and into collective motions, that contribute to the width of transitions, this being typically very important in large or extended molecules. The dynamics also includes transitions of the systems between conformers. Our examination of the MD trajectory at ∼300 K indeed reveals that such transitions take place. Figure 11 shows several snapshots taken from a trajectory of R-D-glucose at ∼300 K. These snapshots demonstrate that a transition between rotamers takes place after t ∼ 3.0 ps, and it is accompanied by the formation of a new internal H-bond. The energies of the conformers are such that each of these two structures is expected to contribute to the spectrum at the temperature studied, and indeed this is the case. We conclude here that DFT makes the connection between spectra and dynamics very explicit. It is therefore an extremely
valuable theoretical tool to study dynamics and vibrational spectra at finite temperatures, as it inherently accounts for the presence of multiple conformers if the temperature makes such transitions possible. In our assessment, the effects of conformational transitions and the effect the population of different conformers have on the spectra are most important in the gas phase. In the crystalline phases studied in section III.C, the molecules are more rigid, and conformational transitions are relatively infrequent at moderate temperatures. E. Properties of Anharmonic Potentials in Glucose and Sucrose. VSCF and the classical dynamics calculations of spectroscopy both include anharmonic effects. Their role is indirectly revealed in the good agreement obtained with experiment in “first principles” vibrational calculations, that is, in methods that proceed from given potentials to calculations of spectra. The anharmonic effects come only from the potential in these calculations, and empirical corrections for harmonic frequency values, “anharmonic scaling”, which is used by many,14,70-74 are avoided. While the harmonic contribution to the fundamental transition is in general the dominant one by far, the harmonic approximation itself results in much larger discrepancies with experiment than the anharmonic calculations presented here. We have employed the representation of the potential surface as given in eq 1 for analyzing the anharmonic effects. This form describes two types of anharmonicity: one is associated with a single normal mode contribution (the terms Vi (Qi), where i is the mode label), and the other is due to anharmonic coupling Vij (Qi, Qj) between a pair of normal modes Qi, Qj. Both anharmonicities contribute significantly, in different weights, to the anharmonic frequency shifts of different modes. The physical significance is, however, different since in dynamical processes only the Vij terms can induce vibrational energy flow between different modes. It is therefore of importance to characterize the strengths of the anharmonic couplings between the different modes. We carried out such an analysis for R-D-glucose, β-D-glucose, and phenyl-β-Dglucose. The magnitude of the pairwise mode-mode couplings can be characterized by the average absolute value over {Qi, Qj} over the relevant vibrational wave functions. We refer here to the average. The largest average couplings found are between OH stretch modes and other OH stretch modes, which in some of the cases have values near 1 kcal/mol. OH stretches can also be strongly anharmonically coupled to shearing/torsional modes. The strongest such pairwise mode-mode coupling we noticed of 1.33 kcal/mol is between an 5867
dx.doi.org/10.1021/jp110043k |J. Phys. Chem. A 2011, 115, 5859–5872
The Journal of Physical Chemistry A
ARTICLE
Figure 10. Band assignments from DFT trajectories: OH modes for R-D-glucose (left) and β-D-glucose (right), at ∼50 (bottom), at ∼150 (middle), and at ∼300 K (top).
Figure 11. Snapshots from CP2K/DFT MD trajectory at ∼300 K, for R-D-glucose, showing transitions between a number of different rotamers taking place after ∼3 ps from the start of trajectory.
OH stretch and a torsional mode in the G-g conformer of R-Dglucose. OH stretching modes can also be significantly coupled to CCO bending, OCO bending, or to CO stretching modes. The magnitude of these was
