The Vibrational Spectrum of Parabanic Acid by Inelastic Neutron

Feb 18, 2010 - The incoherent inelastic neutron scattering spectrum of parabanic acid was measured and ... Neutron diffraction suggests nearly symmetr...
0 downloads 0 Views 3MB Size
Download PDF
3630

J. Phys. Chem. A 2010, 114, 3630–3641

The Vibrational Spectrum of Parabanic Acid by Inelastic Neutron Scattering Spectroscopy and Simulation by Solid-State DFT Matthew R. Hudson, Damian G. Allis, and Bruce S. Hudson* Department of Chemistry, 1-014 Center for Science and Technology, Syracuse UniVersity, Syracuse, New York 13244-4100 ReceiVed: December 1, 2009; ReVised Manuscript ReceiVed: January 29, 2010

The incoherent inelastic neutron scattering spectrum of parabanic acid was measured and simulated using solid-state density functional theory (DFT). This molecule was previously the subject of low-temperature X-ray and neutron diffraction studies. While the simulated spectra from several density functionals account for relative intensities and factor group splitting regardless of functional choice, the hydrogen-bending vibrational energies for the out-of-plane modes are poorly described by all methods. The disagreement between calculated and observed out-of-plane hydrogen bending mode energies is examined along with geometry optimization differences of bond lengths, bond angles, and hydrogen-bonding interactions for different functionals. Neutron diffraction suggests nearly symmetric hydrogen atom positions in the crystalline solid for both heavy-atom and N-H bond distances but different hydrogen-bonding angles. The spectroscopic results suggest a significant factor group splitting for the out-of-plane bending motions associated with the hydrogen atoms (N-H) for both the symmetric and asymmetric bending modes, as is also supported by DFT simulations. The differences between the quality of the crystallographic and spectroscopic simulations by isolated-molecule DFT, cluster-based DFT (that account for only the hydrogen-bonding interactions around a single molecule), and solid-state DFT are considered in detail, with parabanic acid serving as an excellent case study due to its small size and the availability of high-quality structure data. These calculations show that hydrogen bonding results in a change in the bond distances and bond angles of parabanic acid from the free molecule values. 1. Introduction Solid-state DFT methods have proven valuable as an accurate means of simulating vibrational spectra, especially the inelastic neutron scattering (INS) spectra, of molecular solids by accounting for changes to molecular geometry that cannot be properly modeled by nonperiodic approaches. As there are no selection rules governing transitions in INS spectroscopy, an INS spectrum is the complete sum of fundamental transitions of all modes in proportion to H atom motions with contributions from overtones and combinations. The result is a feature-rich spectrum that, with proper assignment, provides considerable detail about the crystal environment. The INS spectral intensity is directly proportional to the amplitude of the atomic motions and the scattering cross section of the atoms involved. For this reason, INS spectroscopy focuses primarily on hydrogen atom dynamics, as the scattering cross section of hydrogen is much larger than other nuclei, including deuterons. This makes INS both a very sensitive tool for the study of molecules that interact in the solid state via hydrogen bonds and a sensitive probe for the characterization of the relative strengths of these interactions, as molecules can be selectively deuterated to elucidate spectral information that may have been previously inaccessible for various reasons, such as vibrational mode coupling. The lack of selection rules for INS stems from the very short-range nature of the interaction and the lack of applicability of the dipole approximation common in IR and Raman spectroscopy. The scattering cross-section for a particular nucleus depends on nuclear properties and on the relative orientation of the nuclear * Corresponding author, [email protected].

and neutron spins. The random nature of this orientation results in an incoherent contribution which, for H atoms, is the dominant term. A very complete description of INS vibrational spectroscopy is available in ref 1. Introductions for those familiar with molecular vibrational spectroscopy in general have also been presented.2-4 In contrast to isolated-molecule calculations, periodic DFT geometry optimizations account for the “propagation” of molecular changes due to crystal packing and stronger intermolecular interactions. Further, solid-state normal-mode analyses also properly address the prediction of phonon modes, factor group splitting, and the changes to vibrational mode energies resulting from crystal packing and stronger intermolecular interactions. In solid-state DFT, the lattice constants determined by neutron or X-ray diffraction are used as input parameters for the geometry optimization process. This raises two important issues. The first of these is the common disparity between the temperature of a vibrational spectrum measurement and the temperature of the diffraction study used as the basis for the solid-state calculations. When these two temperatures are significantly different (as is often the case in INS spectroscopy), the unit cell parameters being used for the calculated normalmode analysis may not yield spectral simulations representative of the sample having its vibrational spectrum measured. This is of specific concern when the interaction energies between molecules in the solid (and, therefore, their vibrational motions), such as hydrogen-bonding interactions, are sensitive to the predicted separation between donor and acceptor groups. The INS spectra of molecular solids are usually taken with the samples at approximately 25 K, a temperature for which

10.1021/jp9114095  2010 American Chemical Society Published on Web 02/18/2010

Vibrational Spectrum of Parabanic Acid

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3631

SCHEME 1: Molecular Structure of Parabanic Acid

diffraction data are often not available. Neutron diffraction data providing accurate hydrogen atom positions are also not commonly available. New facilities will provide increased capability for obtaining neutron diffraction structures with small crystals at low temperature.5,6 The second issue is a consequence of the inexact prediction of molecular geometries that comes about from the approximate nature of the electronic structure treatments of the many density functionals. Specifically, the predicted size of a molecule may differ from the results of the crystallographic determination for reasons that have to do with the accuracy of the treatment of the covalent bonds. For instance, a molecule predicted to be slightly larger than that determined experimentally may shift within a fixed-lattice calculation or tilt due to steric congestion between neighboring molecules, leading to a change in position. It is expected that the largest affect on the distances and angles between atoms involved in the weakest bonds in the crystal, including the intermolecular hydrogen bonds, would suffer the greatest deformation. Parabanic acid (PBA, 1H,3H-imidazolidinetrione) is a hydantoin-derived barbiturate that produces a soporific effect7 and, despite being too readily hydrolyzed to be a viable drug candidate,8 has possible applications as an oxygen radical activity monitor in, for instance, the injured human brain.9 PBA has the same hydrogen-bonding pattern as many barbiturates and forms several of the same co-crystals.10,11 PBA belongs to a large group of molecules that contain only CdO · · · HsN hydrogen-bonding interactions in the solid state, including alloxan, urazole, and barbituric acid, all of which have been investigated in broader studies of crystal polymorphism.12 PBA does not seem to exhibit polymorphism even at elevated pressure,13 although it has been predicted that there are seven alternative structures only 2-6 kJ/mol higher in energy than the one observed.14-18 PBA is of interest from the INS point of view in that it contains only two H atoms, both involved in hydrogen bonds. The H atoms are related to each other by symmetry in the isolated molecule. In isolation, PBA is a planar C2V molecule with three carbonyl groups and two N-H groups (Scheme 1). In crystalline PBA (Figure 1), there are two infinite chains of moderately strong hydrogen bonds along the a-axis as well as parallel chains running in a diagonal “X” along the ac-plane to form hydrogenbonded two-dimensional sheets.19-21 This dense network of hydrogen bonds in a small unit cell volume, as well as the abundance of diffraction data, makes PBA an excellent system for testing the limits of DFT to accurately model extended hydrogen bonding in molecular solids. The INS spectrum of PBA at 25 K is compared with solid-state (and isolatedmolecule) DFT calculations and discussed in relation to: (1) the performance of various functionals in treating the vibrations of the solid, (2) the effect of environment on the out-of-plane vibrations in the solid, (3) optimization of the unit cell from the values provided by the diffraction experiment at 123 K, (4) a hydrated-molecule model (Scheme 2) to remove the possible influence of molecular volume changes in the solid-state calculations, and (5) possible anharmonicity in the covalent

Figure 1. A hydrogen-bonded PBA five-member cluster with atomic labeling scheme and unit cell views down each crystal axis. The two NsH · · · OdC hydrogen bonds are shown and form chains along the a-axis.

SCHEME 2: Hydrated-Molecule Structure of Parabanic Acid

N-H bonds and its influence on NsH · · · O bending mode energy predictions. 2. Experimental and Computational Methods 2.1. Sample. Parabanic acid was obtained from SigmaAldrich at 99% purity and used as received. Neutron diffraction data at 123 and 298 K reported by Craven et al.19,21 was used as the initial basis for the solid-state DFT calculations (Cambridge Structural Database ID: PARBAC). PBA crystal data are as follows for the 123 K neutron structure: monoclinic space group P21/n, Z ) 4, a ) 10.704(2) Å, b ) 8.187(2) Å, c ) 4.969(1) Å, β ) 92.32(1)°, V ) 435.1 Å3. The melting point is reported to be 516 K. 2.2. INS and Raman Spectra. The INS spectrum of the powdered material was measured using the time-focusing, crystal analyzer spectrometer TOSCA at the Rutherford Appleton Laboratory, U.K.22 at approximately 25 K over the range of 4-500 meV (32-4000 cm-1) with a resolution of 1.5-3.0% of the energy transfer (∆ω/ω). Raman spectra for the polycrystalline sample were collected at 78 K and room temperature using a Renishaw inVia Raman spectrometer with a 785 nm excitation wavelength. 2.3. Calculations. DFT geometry optimizations and vibrational normal-mode analyses were performed at the Γ-point using DMol3 (v4.4),23,24 GAUSSIAN 03 (rev. E.01),25 and CRYSTAL06 (v1.0.2).26 DMol3 calculations were performed using both the local density approximation (LDA) PWC27 and VWN28 density functionals and the generalized gradient approximation (GGA) BLYP,29,30 BOP,31 BP,29,32 HCTH,33 PBE,34,35 PW91,36 RPBE,37 and VWN-BP28,29,33 density functionals with the DNP basis set. All DMol3 calculations were performed with 3 × 3 × 4 k-point sampling, a program option “fine” grid size of 0.04 Å, and a convergence criterion ∆Eopt < 5 × 10-8 hartree. Normal mode analyses were completed within the harmonic approximation numerically for each DMol3 calculation using

3632

J. Phys. Chem. A, Vol. 114, No. 10, 2010

two-point central differencing to obtain second derivatives. In order to investigate the possible anharmonicity of vibrational modes in PBA, DMol3 normal-mode analyses were completed for the BLYP functional at six different atomic displacements (from 0.05 to 0.11 Å) within the two-point approach to determine the dependence of the Γ-point frequencies on the size of the finite atomic displacements. By this approach, changes in the predicted Γ-point frequencies due to anharmonicity are a result of the different force constants predicted for each displacement (this approach does not yield anharmonic predictions, but it does provide a means to identifying modes for which the potentials are anharmonic based on successive normal-mode analyses within the harmonic approximation). This multipledisplacement approach has been previously demonstrated to be useful for such mode identifications for the case of crystalline urea.38 The CRYSTAL06 program permits unit cell parameter optimization and incorporates hybrid functionals, both features absent from DMol3. CRYSTAL06 calculations were carried out using the B3LYP and PBE0 hybrid density functionals in conjunction with two different Gaussian-type basis sets: 6-31G(d,p) and 6-311G(d,p).39-42 Calculations were performed with program option values as follows: TOLINTEG ) 8 8 8 8 16 (five parameters controlling the truncation criteria of Coulomb and HF exchange integrals, with these values corresponding to thresholds of 10-n), XLGRID (a pruned DFT grid specification with 75 radial points and 974 angular points),43 SCF convergence criterion of ∆Eopt < 10-8 hartree (between subsequent energy cycles), SHRINK ) 6 (a factor applied to set a mesh of points for which total energies are fully converged within the irreducible Brillouin zone, with “6” the optimum sampling of reciprocal space as determined by monitoring the convergence of the total energy as a function of K-point count). Geometry optimizations were performed on the neutron diffraction structure with (1) no lattice cell relaxation, (2) a constant-volume optimization constraint of the unit cell parameters and atomic positions, and (3) a full optimization with no volume constraints on the unit cell. The cell parameter optimization permits relaxation of the hydrogen-bonded PBA chains. The resulting optimized lattice parameters from (2) were also used as the basis for a geometry optimization and normalmode analysis in DMol3 to see what effect, if any, the optimized lattice parameters had on a BLYP/DNP calculation. As with DMol3, Γ-point vibrational frequencies were determined within the harmonic approximation in CRYSTAL06 by diagonalization of a mass-weighted Hessian matrix via numerical second derivatives calculated from analytic first derivatives and finite atomic displacements. A convergence criterion of ∆E < 10-11 Hartree was used for the normal mode calculations. CRYSTAL06 has a default step size of 0.001 Å and uses a centered, three-point differencing method. Isolated-molecule, hydrated-molecule, and PBA pentamer DFT calculations were performed in DMol3 at the [functional]/ DNP level (using the same optimization and normal mode criteria as the solid-state calculations). The isolated-molecule and hydrated-PBA calculations were also performed in CRYSTAL06 at the B3LYP/6-311G(d,p) and PBE0/6-311G(d,p) levels with the same optimization and normal mode criteria used in the solid-state simulations. The GAUSSIAN 03 software25 was also used for each of the aforementioned cases using program options “tight” convergence criteria (density matrix rms 3000 cm-1) and low-frequency (0-250 cm-1) regions will be considered in following sections. Comparing only the molecular vibrational mode region between the periodic and isolated-molecule cases reveals the extent to which the solid-state spectrum is in better agreement with experiment. The solid-state simulated spectrum does not agree well with the INS out-of-plane hydrogen bending mode energies (600-900 cm-1), the most intense features in this spectrum. Solid-state DFT correctly accounts for the factor group splitting of the out-of-plane hydrogen bending modes but overestimates the frequency of these modes. The relationship between the hydrogen-bonded PBA units and solid-state DFT functional(s) is presented in Figure 5 and is discussed below. The Raman spectrum of PBA is presented in Figure 6 along with the INS spectrum for comparison. The two spectroscopies can profitably be used in tandem. In Figure 6, the carbonyl

Vibrational Spectrum of Parabanic Acid

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3635

Figure 5. The incoherent INS vibrational spectrum compared to simulated solid-state spectra (red) for both LDA and GGA functionals from 400 to 1600 cm-1. All spectra simulations shown were completed using the DMol3 software program at a [functional]/DNP level of theory.

Figure 6. The INS spectrum of PBA acid (black, upper) and room temperature Raman spectrum (red, lower) between 0 and 2000 cm-1.

(CdO) stretching vibrations between 1800 and 1850 cm-1 are dominant, sharp features that are essentially non-existent in the INS spectrum, as the region is composed of overtones, combinations, and carbonyl stretching fundamentals. The CdO stretching modes, which are also weak in the INS spectrum due to the scattering cross sections of carbon and oxygen compared to hydrogen, can be assigned from the Raman spectrum. There are three sharp features in the Raman spectrum at 660, 667, and 804 cm-1. These modes are in-plane ring deformation modes which are “not visible” in the INS spectrum as the out-of-plane hydrogen bending modes dominate the spectrum. Alternatively, of the two fundamental hydrogen bending modes, only one is very weakly Raman active and is a shoulder at 689 cm-1. The Raman spectrum does not reveal the factor group splitting as clearly as the INS spectrum of PBA due to optical selection rules.

Figure 7. The INS spectrum (black) is compared with the full BLYP/ DNP and RPBE/DNP solid-state DFT simulations (red) showing the calculated fundamentals (blue) and second (green) and third (orange) overtone transitions through 3500 cm-1.

3.3. Overtones and Combinations. Vibrational motions that involve the largest hydrogen atom displacements along the normal mode coordinate generate the most intense INS features. Since there are no selection rules in INS spectroscopy, overtone and combination features can occur as clearly visible features that rival the intensities of fundamental vibrational modes at higher frequency (Figure 7). The most notable features in the PBA overtone spectrum are the 0 f 2 transitions of the outof-plane hydrogen bends. There is a grouping of peaks at slightly higher energy than the in-plane hydrogen bending fundamental transitions between 1400 and 1750 cm-1, corresponding to the first overtones of the phase-split out-of-plane hydrogen bending modes. Unsurprisingly, the 0 f 2 transitions for the out-ofplane bending modes are incorrectly simulated as the 0 f 1 fundamentals are overestimated. Several other 0 f 2 and

3636

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Hudson et al.

Figure 8. The phonon and low-frequency vibrational (3000 cm-1). A summary of the IR, Raman, and INS features above 3000 cm-1 and the DFT predictions are provided in Table 2. The high-frequency, hydrogen-stretching region of PBA (Figure 7) contains one broad feature centered at approximately 3200 cm-1 with three discernible peaks. These are attributed to one B2 mode (asymmetric stretching) and a splitting of an A1 mode (symmetric stretching) into two peaks. This assignment is based on previous IR and current Raman spectra, as the INS spectrum lacks the necessary clarity due to both the number of transitions in this region and resolution limits. There is a splitting of over 60 cm-1 for the two modes (A1, B2) in the Raman experiment. The isolated-molecule calculations predict two features in this region separated by only 1 cm-1, while the periodic simulation predicts eight modes spread over a range of 165 cm-1 that take into account the splitting of modes in the solid-state (Z ) 4). The solid-state simulation gives an A1 mode splitting of 70 cm-1, but the B2 vibration is predicted to split by only 24 cm-1. 3.5. Phonon and Low-frequency Molecular Region (0-250 cm-1). The experimental and simulated INS low-frequency region of PBA is shown in Figure 8. The isolated-molecule DFT simulation predicts two features in this region: a full molecule bend and a ring twisting (deformation) mode. As with the highfrequency region, the solid-state DFT calculation contains eight modes over approximately 40 cm-1. The simulation correctly splits the vibrations for the full molecule-bending motions centered at 212 cm-1 but over-splits the modes corresponding to ring deformation with O4/O5 out-of-plane bending. The calculation splits these four modes by 40 cm-1, (two peaks between 175 and 215 cm-1), while the broad feature in the INS

spectrum assigned to this motion spans approximately 20 cm-1 and is centered at 193 cm-1. The calculated and INS peak positions are found in Table 2. Parabanic acid has 9 translational optical modes and 12 rotational/librational modes (plus three acoustic modes) for a total of 21 features in the optical phonon region of the Γ-point simulation (plus the acoustic modes away from the Γ-point in INS). The solid-state simulation predicts these modes to all lie below 180 cm-1. As a Γ-point calculation only, the [functional]/ DNP simulated spectra do not accurately reproduce the distribution of intensity in this region, where the INS spectrum contains several broad features. No detailed assignments can be made. 4. Discussion PBA has only one known polymorphic form between 123 K and decomposition (mp ) 516 K) and from ambient pressure up to 2.1 GPa.13 Very thorough computational studies have been performed to predict the global minimum energy solid-state arrangement, yielding the neutron diffraction structure as the lowest energy form with very few other possible structures.14-16,18 Solid PBA is quite stable despite the use of only two of the three hydrogen-bond acceptors in the molecule. On the basis of previous studies and very good agreement between the Raman spectral features at 78 K and the 25 K INS spectrum, it is unlikely that the INS spectrum is that of something other than PBA in its known crystal form despite the difference in temperature between the INS and diffraction experiments. The degree of agreement between the computed structures and the diffraction experiment for PBA is of interest in its own right. Table 1 shows the degree of agreement found between the diffraction structure and that computed for the isolated molecule using four methods. In the order BLYP/DNP (DMol3), B3LYP/6-311G(2d,2p) (Gaussian03), B3LYP/cc-PVQZ (Gaussian03), and MP2/6-311G(2d,2p) (Gaussian03), the RMSDs for the bond lengths are 0.0232, 0.0221, 0.0219, and 0.0210 Å. If the two N-H bonds are deleted from this calculation, the values

Vibrational Spectrum of Parabanic Acid are 0.0217, 0.0157, 0.0154, and 0.0140 Å. The average deviation for the heavy atom bond lengths are 0.0162, 0.0015, 0.0001, and 0.0048 Å. This shows that the DNP basis set predicts the molecule to be too large overall by ca. 0.02 Å. While this is improved with the larger basis sets, this is due to underestimating the CdO bond lengths considerably; the RMSD values remain near 0.02 Å. If the PBA is surrounded by four water molecules constrained to C2V symmetry, the RMSD for BLYP/DNP (DMol3) for all atoms is 0.0152 Å compared to the 0.0232 Å value for the isolated molecule. The hydrated PBA species remains too large, with a mean heavy atom deviation of 0.016 Å compared to 0.0162 Å for the isolated molecule. However, in this case, a modest improvement in the basis set, with the B3LYP/6311G(2d,2p) (Gaussian03), gives an RMSD of 0.0056 Å (0.0031 Å with N-H bonds removed) and mean heavy atom deviation of 0.0005 Å with the CdO bond length deviations being 0.0025 Å instead of 0.0250 Å as they were for the larger basis set treatments of the isolated molecules. The basis set seems to be of greater importance rather than the functional, with some outliers, and, by comparison with the results for the isolated molecule, inclusion of the hydrogen bonding appears to be more important than the basis set. The DNP treatment of the periodic problem naturally includes the hydrogen bonding. With the BLYP functional, the RMSD is 0.0093 Å (0.0092 Å for heavy atoms) for bond lengths and the mean deviation for the heavy atoms is 0.0081 Å with all deviations being positive. For the BP functional, these values are 0.0075 Å RMSD and 0.0043 Å mean deviation. This is one of the lowest mean square deviations observed of theory and experiment for a hydrogen-bonded structure. Full tables of bond lengths and angles and corresponding deviations from the experimental structure are presented in the Supporting Information. The inaccurate reproduction of the out-of-plane bending mode energies in the solid-state DFT calculations motivated several additional modeling efforts: (1) varied step-size to sample the bend potential further from the minimum in the normal-mode analyses to consider the possible anharmonicity of the out-ofplane vibrational modes, (2) fixed lattice constant molecular geometry optimizations with multiple functionals and basis sets, (3) full optimization of the unit cell parameters and molecular geometry at constant cell volume, and (4) full optimization of the unit cell parameters and molecular geometry without the constant volume constraint (completely unrestricted optimization). 4.1. Anharmonicity. The INS features associated with the out-of-plane bending motions of the hydrogen-bonded H atoms are noticeably and almost unvaryingly overestimated in the solid-state DFT calculations, while other modes are accurately simulated at all levels of employed theory. The possibility that the bending modes are very anharmonic in nature, and thus overestimated in the harmonic approximation, was tested as done previously in ref 38. Varying the step size of the displacement of the harmonic normal-mode analysis leads to a change in the predicted energies of anharmonic normal modes because the resulting second derivative samples further from the minimum. Figure 9 shows several of the modes for the out-of-plane bending region (and one of the in-plane bending modes at approximately 1400 cm-1). Considering the out-of-plane hydrogen-bending modes, not a single mode is shifted in energy by more than a few wavenumbers. The hydrogen stretching modes (not shown in Figure 9) shift down in energy approximately 100 cm-1 upon varying the displacement of the

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3637

Figure 9. Dependence of Γ-point frequencies on the size of atomic displacements due to the anharmonicity of the predicted vibrational mode (see text).

normal-mode analysis, an unsurprising result but demonstrative of the sensitivity of the predicted mode energy to step-size within the harmonic approximation. The absence of mode energy change in the out-of-plane modes with differing step-size confirms that the bending motions are no more anharmonic in nature than any of the other spectral features in regards to the calculated positions. This would suggest that the nature of the experimentally determined vibrations of the out-of-plane bending region are no more anharmonic in nature when compared to the other available vibrational modes as demonstrated by the step-size variation simulations. Interestingly, in each of the DMol3 calculations, the hydrogen stretching modes are “accurately” calculated even though these modes are known, and found, to be anharmonic. Due to resolution limitations of the INS experiment, and the abundance of overtones and combinations, not much weight can be given to this agreement. This is not to say that the experimentally determined features are not anharmonic. However, anharmonicity is not the likely cause of disagreement between the experimental and simulated spectra. 4.2. Fixed-Cell INS Simulations. It has been shown that the BLYP functional simulates experimental INS spectra more accurately than many other density functionals (specifically among those density functionals implemented in the DMol3 program).49-51 For this investigation, PBA was optimized and normal-mode analyses completed for every available density functional in the DMol3 software program, encompassing both the GGA and LDA, and are shown in Figure 5. The HCTH functional is included but will not be discussed further in the comparisons, as this functional yields very poor spectral agreement with experiment (Figure 5) and equally poor optimization of the molecular geometry. Both LDA and GGA density functionals yield simulated spectra in good agreement with experiment, including C-C and C-N stretching and C-C and C-N torsional modes. Both the LDA and GGA functionals accurately simulate the in-plane hydrogen-bending motions. The in-plane modes appear as a broad asymmetric feature at approximately 1400 cm-1. On comparison of the out-of-plane hydrogen-bending modes, the GGA functionals yield better mode energies and model the observed factor-group splitting more consistently. With the exception of the HCTH functional, the GGA functionals result in molecular geometries that are more in-line with the neutron diffraction experiment (Figures 2 and 3) than found in the two LDA calculations.

3638

J. Phys. Chem. A, Vol. 114, No. 10, 2010

Focusing on the GGA functionals, there is almost uniform overestimation of the out-of-plane hydrogen-bending modes. The RPBE functional is the only functional to reproduce even one of the eight groupings of the out-of-plane hydrogen-bending modes. If we consider the BLYP, BOP, and RPBE functionals, all three result in optimized bond lengths that are nearly identical to each other with reasonable RMSDs from the diffraction experiment (Figure 2), while the optimized bond angles differ significantly among these functionals (Figure 3). The largest deviation among the functionals in terms of molecular geometries is in the C(2/4)-N3-H3 angles, with the RPBE functional most accurately reproducing the diffraction experiment with respect to this variable, followed by the BOP functional. A comparison of the hydrogen bond lengths and angles for the DMol3 calculations is found in the Supporting Information, including the differences with the experimental diffraction results. In each case, the crystallographically determined difference in NsH · · · O angles for H1 and H3 are maintained. When the deviations from experiment are considered, several trends become clear. The first is that each functional (focusing again on the GGA functionals without HCTH) calculates the N1-H1 distance to within 0.006 Å of experiment but fails to calculate the N3-H3 distances to better than 0.01 Å. In the case of H3, whose motions are responsible for the lower-energy out-of-plane bending modes, the RPBE functional gives an NsH · · · O angle that is smaller than the experimental value and H3 · · · O5 and N3 · · · O5 distances longer than experimentally observed. Each of the other functionals (again discounting the HCTH functional) yields H3 · · · O5 and heavy-atom (N3 · · · O5) distances that range from accurate to underestimated as compared to experiment. With H1, whose motions are responsible for the higher-energy out-of-plane bending modes, the N1sH1 · · · O4 angles are in reasonable agreement, but the N1 · · · O4 and H1 · · · O4 distances are consistently too short. The result of these differences in bond angles and distances is readily apparent in the simulations (Figure 5). The slight overestimation of the H3 hydrogen-bonding distances (and underestimation of bond angle) in the RPBE function is presumed to give rise to an accurate model of the low-energy out-of-plane bending modes, while the underestimation of the H1 hydrogen-bonding distances is presumed to give rise to simulations whose energy is too high compared to experiment because a stronger hydrogen bond results in a higher frequency out-of-plane mode. This point is further illustrated in Figure 4, which compares a range of PBA simulations from isolated-molecule to the solid-state cases. The full functional comparison is given as bond lengths and angles in the Supporting Information. There are two out-of-plane bend vibrations in the isolated molecule, one corresponding to a symmetric proton motion with no dipole change and one corresponding to antisymmetric proton motion, where each H atom moves together along one axis, resulting in a net change in dipole moment in the out-of-plane direction. In the crystalline solid, due to the angle of inclination between neighboring molecules and the difference in hydrogenbonding angle, the two hydrogen atoms of PBA are different. The difference in N-H bond lengths between N1-H1 and N3-H3 is 0.004 Å, between H1 · · · O4 and H3 · · · O5 is 0.05 Å, and between N1 · · · O4 and N3 · · · O5 is 0.02 Å. The angles between the neighboring PBA molecules for H1 vs H3 via the hydrogen bonding are different by 6° and propagate along the n-glide plane. The result of this asymmetry is large site-group splitting in the out-of-plane hydrogen-bending motions, as well as a very small dipole moment in the “symmetric” modes which should be accessible in optical spectra as well as the INS and

Hudson et al.

Figure 10. The INS spectrum (black) is compared with the solidstate DFT simulations (red) computed using the CRYSTAL06 software program with B3LYP and PBE0 functionals with 6-31G(d,p) and 6-311G(d,p) basis sets. The B3LYP/6-311G(d,p) simulated spectrum marked with an asterisk is the simulated spectrum for the constant volume optimization of the unit cell parameters and molecular geometry optimization.

is seen as a very weak shoulder in the Raman spectrum at approximately 689 cm-1 (Figure 6). The next step in the investigation of the inability of the DMol3 solid-state DFT simulations to accurately model the bonds, angles, and energies of the out-of-plane hydrogen-bending modes is the exploration of hybrid density functionals or planewave (PW) methods. The use of PW methods was not explored in part due to recent findings for INS simulations for hydrogen-bonded formic and acetic acids.43,52 The CRYSTAL06 program permits the use of hybrid functionals, including B3LYP and PBE0. The results of these calculations are presented in Figure 10 and bond lengths and angles with corresponding RMSDs from experiment given in the Supporting Information. Hydrogen-bonding distances and angles with differences from experimental values are given in Table 3. The difference in the two hydrogen-bonding NsH · · · O angles of approximately 6° is maintained. The N-H distances differ between nonhybrid and hybrid calculations in that the B3LYP and PBE0 functionals result in distances that are almost identical, resulting in a PBA molecular symmetry closer to C2V (and further from that found in the crystal) than found with the GGA and LDA functionals. As with the GGA/LDA functionals, the simulation that most closely resembles the experimental spectrum is that for which one of the two heavy atom distances in the hydrogen bonds is overestimated, generated at the PBE0/6-311G(d,p) level of theory. The resulting PBE0/6-311G(d,p) simulation overestimates the H1 · · · O4 distance and underestimates the H3 · · · O5 distance (Table 3), while the RPBE/DNP simulation results in the H3 · · · O5 distance being overestimated and the H1 · · · O4 distance being underestimated. The nonhybrid RPBE/DNP INS simulation remains the one that has the best fit to experiment of the fixed unit cell parameters based on the 123 K neutron diffraction data. 4.3. Constant Volume and Full Cell Optimization INS Simulations. Geometry optimizations were performed using CRYSTAL06 with both full optimization of the lattice parameters and optimization of the lattice parameters under constant cell volume constraints. The resulting lattice constants for these calculations can be found in Table 4. The bond lengths and angles are found in the Supporting Information, and the

Vibrational Spectrum of Parabanic Acid

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3639

TABLE 3: Summary of the Hydrogen-Bonding Distances (in angstroms) and Angles (in degrees) from 123 K Neutron Diffraction21 and Solid-State DFT Calculations with Corresponding Differencesa B3LYP/ 6-31G(d,p)

B3LYP/ 6-311G(d,p)

PBE0/ 6-311G(d,p)

N3sH3, Å H3 · · · O5, Å O5sN3, Å O5sH3sN3, deg N1sH1, Å O4 · · · H1, Å O4sN1, Å N1sH1sO4, deg

1.0273 1.8206 2.8277 165.771 1.0293 1.8125 2.8337 170.962

1.0244 1.8437 2.8437 164.385 1.0274 1.8098 2.8288 170.83

1.0237 1.8472 2.8471 164.573 1.0265 1.8240 2.8424 170.964

N3sH3, Å H3 · · · O5, Å O5sN3, Å O5sH3sN3, deg N1sH1, Å O4 · · · H1, Å O4sN1, Å N1sH1sO4, deg

-0.0197 -0.0345 -0.0266 1.771 -0.0137 0.0047 0.0051 0.962

-0.0226 -0.0114 -0.0106 0.385 -0.0156 0.0020 0.0002 0.830

-0.0233 -0.0079 -0.0072 0.573 -0.0165 0.0162 0.0138 0.964

0.017 0.025 0.019 1.43

0.019 0.008 0.007 0.65

0.020 0.013 0.011 0.79

NsH, Å H · · · O, Å OsN, Å OsHsN, deg a

B3LYP/6-31G(d,p) const. vol. optimization

B3LYP/6-311G(d,p) const. vol. optimization

B3LYP/6-311G(d,p) full param. optimization

PBE0/6-311G(d,p) Full param. optimization

RBLcorrected expt21

Hydrogen Bonding 1.0283 1.0247 1.8001 1.8297 2.8138 2.8350 167.886 166.052 1.03 1.0281 1.7862 1.7859 2.808 2.8060 170.921 170.972

1.0255 1.8557 2.8739 171.509 1.0254 1.8548 2.8732 171.665

1.0256 1.8237 2.8410 170.884 1.0258 1.8472 2.8395 171.336

1.047 1.86 2.85 164 1.043 1.81 2.83 170

Differences from Experiment -0.0187 -0.0223 -0.055 -0.0254 -0.0405 -0.0193 3.886 2.052 -0.013 -0.0149 -0.0216 -0.0200 -0.0206 -0.0226 0.921 0.972

-0.0215 0.0006 0.0196 7.509 -0.0176 0.0470 0.0446 1.665

-0.0214 -0.0314 -0.0133 6.884 -0.0172 0.0372 0.0109 1.336

0.020 0.033 0.034 5.44

0.019 0.036 0.012 4.96

RMSD Values 0.016 0.042 0.032 2.82

0.019 0.024 0.021 1.61

Neutron diffraction values have been previously corrected for rigid-body librations (RBL).

TABLE 4: Unit Cell Lattice Parameters and Cell Volumes (Å3) for PBA from the 123 K Neutron Diffraction Experiment16 and Values Calculated from Optimization of These Parametersa optimized lattice parameters and volume

expt21 a, Å 10.704 b, Å 8.187 c, Å 4.969 beta, deg 92.32 3 volume, Å 435.45

difference from experimental value

% of experimental value

PBE0/ B3LYP/ B3LYP/ PBE0/ B3LYP/ B3LYP/ PBE0/ B3LYP/ B3LYP/ 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) 6-311G(d,p) const. vol. full param. full param. const. vol. full param. full param. const. vol. full param. full param. optimization optimization optimization optimization optimization optimization optimization optimization optimization 10.645 8.097 5.050 91.84 435.23

10.745 8.582 5.234 89.97 482.65

10.665 8.475 5.133 90.08 463.95

-0.059 -0.090 0.081 -0.476 -0.217

0.041 0.395 0.265 -2.350 47.194

-0.039 0.288 0.164 -2.240 28.499

0.549 1.108 1.607 0.517 0.050

0.382 4.711 5.195 2.578 10.281

0.365 3.457 3.247 2.456 6.337

a Differences between the calculated and experimental values for the lattice parameters and cell volume are given as well as the total percent difference between experiment and theory.

hydrogen-bonding information is found in Table 3. Due to the volume constraints on the system, there is no significant change in the cell parameters (the beta angle is closer to 90° by almost 1°). From Figure 10, the constant-volume optimized cell more correctly accounts for the intensities of the factor group-split out-of-plane bending modes than the hybrid functional simulations with the fixed neutron diffraction parameters, but results in higher-energy positions for the “low-energy” out-of-plane bend modes. The hydrogen-bonding H3 · · · O5 and N3 · · · O5 distances, as well as the N3sH3 · · · O5 angle, are indeed in worse agreement with the diffraction experiment than the corresponding B3LYP functional simulations with no cell parameter optimization. A full optimization of the lattice parameters and atomic positions was performed using both the B3LYP and PBE0 hybrid density functionals. A comparison of the resulting B3LYP/6-311G(d,p) and PBE0/6-311G(d,p) cell parameters to the 123 K neutron diffraction results can be found in Table 4. The corresponding bond lengths and angles are found in the Supporting Information, and the hydrogen-bonding information is in Table 3. Simulated INS plots from these calculations are not provided, as the results contain multiple negative vibrational modes (a consequence of the full cell optimization being

performed under the symmetry constraints of the PBA space group). The beta angle in each case is approximately 90° and the distances between molecules are essentially identical and symmetric. In gauging the degree of failure between the two full lattice parameter optimization simulations, the PBE0 optimization is consistently more accurate than the B3LYP optimization. The RMSD values for bond lengths and bond angles are lower for the PBE0 lattice parameter and cell volume optimizations than for B3LYP and inline with the values for both B3LYP and PBE0 geometry optimizations with experimental lattice parameters. The lattice parameters for the PBE0 full cell optimization (Table 4) are slightly better than in the B3LYP equivalent optimization, with this improvement most noticeable in the resulting cell volumes. There is almost a 50 Å3 enlargement of the B3LYP full-optimization unit cell, (11%) while the PBE0 optimization results in a 28.5 Å3 (6.6%) enlargement and better bond length and angle agreement with experiment. The constant volume cell optimization (B3LYP) leads to much better RMSDs for bond lengths and angles compared to the full lattice parameter optimizations regardless of functional. The values for the hydrogen bonding, especially the intermolecular hydrogen-bonding angles, are more accurately simulated

3640

J. Phys. Chem. A, Vol. 114, No. 10, 2010

using the constant volume optimization in place of a free lattice parameter optimization (Table 3). However, these values are not as good as those deduced from the simulations using fixed experimental lattice parameters. 5. Conclusions This study of PBA reveals quite convincingly that inclusion of intermolecular interactions greatly improves the degree of agreement between the computed and observed structure for this material. This is most clearly shown for the bond-angle variation. A reasonable explanation for this effect is that it reflects an actual structural difference between the gas-phase molecule and the molecule in the presence of intermolecular interactions. The fact that hydrated-molecule calculations provide most of this improvement in agreement between calculation and experiment supports this conclusion. It is interesting to note the parallel between this observation and previous studies of N-methylacetamide leading to the same conclusion, as both systems contain the peptide linkage.53 These observations bear, to some extent, on the issue of the development of empirical potential functions for molecular mechanics force fields as pertaining to the comparison of computations with experiment. Whereas gas-phase comparisons between theory and experiment provide a direct measure of the limits of a theoretical model for reproducing gas-phase experimental data, force field development employing quantum chemical methods but based on solid-state experimental data risks propagating structure and dynamics errors due not to the theoretical model alone but instead to the neglect of the molecular environment in the validation process. Parabanic acid is a small system in which the differences between gas-phase and solid-state geometry and spectroscopic properties are experimentally observable, making the resulting agreement between isolated-molecule calculations and crystallographic data an inappropriate metric for gauging the accuracy of quantum chemical calculations that do not include intermolecular interactions. In the specific case of peptide bonds considered here and in earlier work, the bond lengths internal to the molecule change when hydrogen bonds form. This has to be included in any force field that is to be applied to condensed phases. The high-resolution, low-temperature crystallographic data show N1-H1 and N3-H3 distances of 1.043 and 1.047 Å, respectively, with heavy atom N1-O4 and N3-O5 distances of 2.83 and 2.85 Å, respectively The O4-H1 and O5-H3 hydrogen-bonding distances are 1.81 and 1.86 Å, respectively, and the N1sH1 · · · O4 and N3sH3 · · · O5 angles are 170 and 164°, respectively. A major conclusion of this study is that it is necessary to describe this asymmetric geometry accurately in order to obtain agreement between experiment and calculation. It is clear from this work that this applies to nonhybrid and hybrid density functionals alike and that optimizing the unit cell parameters, at constant volume or free, will not achieve this description more accurately than use of the diffraction parameters. The consideration of model systems leads us to the suggestion that the problem with these calculations is that the methods used result in molecular sizes that are incompatible with the lattice dimensions. This is probably due to limitations of the basis set. Acknowledgment. The authors wish to thank Professor Timothy M. Korter of Syracuse University for computational resources and assistance. We also thank the ISIS Facility of the Rutherford Appleton Laboratory for neutron beam time on the TOSCA spectrometer as well as the Department of

Hudson et al. Chemistry at the State University of New York at Albany for the use of their Raman spectrometer. D.G.A. thanks the office of the Vice President for Research at Syracuse University for research support. This work was partially supported by NSF Grant CHE-0848790 and the National Center for Supercomputing Applications under Grant TG-DMR080000N and utilized the SGI Altix [Cobalt] system. http://www.ncsa.uiuc.edu/ UserInfo/. Supporting Information Available: Further tables of calculated bond lengths and angles, corresponding RMSD’s, and structure comparisons, the 78 K Raman spectral comparison with the 25 K INS spectrum (Figure S1), and the full 0-3500 cm-1 simulated and experimental INS spectral comparison (Figure S2). This information is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Mitchell, P. C. H.; Parker, S. F.; Ramirez-Cuesta, A. J.; Tomkinson, J. Vibrational Spectroscopy with Neutrons; World Scientific: Singapore, 2005. (2) Hudson, B. S. In Frontiers of Molecular Spectroscopy, Laane, J., Ed.; Elsevier: Amsterdam, 2009, pp 597-622. (3) Hudson, B. S. Vibrational Spectroscopy 2006, 42, 25. (4) Hudson, B. S. J. Phys. Chem. A 2001, 105, 3949. (5) http://neutrons.ornl.gov/instrument_systems/beamline_12_topaz/index. shtml. (6) Koetzle, T. F.; Piccoli, P. M. B.; Schultz, A. J. Nucl. Instrum. Methods Phys. Res., Sect. A 2009, 600, 260. (7) Lewis, H. B. J. Biol. Chem. 1915, 23, 281. (8) Andrews, J. C.; Sell, I. T. Arch. Biochem. Biophys. 1955, 56, 405. (9) Hillered, L.; Persson, L. Parabanic acid for monitoring of oxygen radical actiVity in the injured human brain; Department of Clinical Chemistry: Uppsala University Hospital, Sweden, 1995. (10) Weber, H. P.; Craven, B. M. Acta Crystallogr., Sect. B: Struct. Sci. 1987, B43, 202. (11) Weber, H. P.; Ruble, J. R.; Craven, B. M.; McMullan, R. K. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1980, B36, 1121. (12) Price, S. L. Int. ReV. Phys. Chem. 2008, 27, 541. (13) Fabbiani, F. P. A.; Allan, D. R.; Marshall, W. G.; Parsons, S.; Pulham, C. R.; Smith, R. I. J. Cryst. Growth 2005, 275, 185. (14) Day, G. M.; Chisholm, J.; Shan, N.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2004, 4, 1327. (15) Day, G. M.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2005, 5, 1023. (16) Lewis, T. C.; Tocher, D. A.; Day, G. M.; Price, S. L. CrystEngComm 2003, 5, 3. (17) Lewis, T. C.; Tocher, D. A.; Price, S. L. Cryst. Growth Des. 2005, 5, 983. (18) Price, S. L. AdV. Drug DeliVery ReV. 2004, 56, 301. (19) Craven, B. M.; McMullan, R. K. Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem. 1979, B35, 934. (20) Davies, D. R.; Blum, J. J. Acta Crystallogr. 1955, 8, 129. (21) He, X. M.; Swaminathan, S.; Craven, B. M.; McMullan, R. K. Acta Crystallogr., Sect. B: Struct. Sci. 1988, B44, 271. (22) Colognesi, D.; Celli, M.; Cilloco, F.; Newport, R. J.; Parker, S. F.; Rossi-Albertini, V.; Sacchetti, F.; Tomkinson, J.; Zoppi, M. Appl. Phys. A: Mater. Sci. Process. 2002, 74, 64. (23) Delley, B. J. Chem. Phys. 1990, 92, 508. (24) Delley, B. J. Chem. Phys. 2000, 113, 7756. (25) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al.Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03; ReVision C.02; Gaussian, Inc.: Wallingford, CT, 2004.

Vibrational Spectrum of Parabanic Acid (26) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, I. J.; D’Arco, P.; Llunell, M. CRYSTAL06 User’s Manual; University of Torino: Torino, 2006. (27) Perdew, J. P.; Wang, Y. Phys. ReV. B: Condens. Matter Mater. Phys. 1992, 45, 13244. (28) Vosko, S. H.; Wilk, L.; Nusair, M. Can. J. Phys. 1980, 58, 1200. (29) Becke, A. D. Phys. ReV. A: At., Mol., Opt. Phys. 1988, 38, 3098. (30) Lee, C.; Yang, W.; Parr, R. G. Phys. ReV. B: Condens. Matter Mater. Phys. 1988, 37, 785. (31) Tsuneda, T.; Suzumura, T.; Hirao, K. J. Chem. Phys. 1999, 110, 10664. (32) Perdew, J. P.; Yue, W. Phys. ReV. B: Condens. Matter Mater. Phys. 1986, 33, 8800. (33) Boese, A. D.; Handy, N. C. J. Chem. Phys. 2001, 114, 5497. (34) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1997, 78, 1396. (35) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (36) Perdew, J. P.; Chevary, J. A.; Vosko, S. H.; Jackson, K. A.; Pederson, M. R.; Singh, D. J.; Fiolhais, C. Phys. ReV. B: Condens. Matter Mater. Phys. 1992, 46, 6671. (37) Hammer, B.; Hansen, L. B.; Norskov, J. K. Phys. ReV. B: Condens. Matter Mater. Phys. 1999, 59, 7413.

J. Phys. Chem. A, Vol. 114, No. 10, 2010 3641 (38) Johnson, M. R.; Parlinski, K.; Natkaniec, I.; Hudson, B. S. Chem. Phys. 2003, 291, 53. (39) Adamo, C.; Barone, V. J. Chem. Phys. 1999, 110, 6158. (40) Becke, A. D. Phys. ReV. A: At., Mol., Opt. Phys. 1988, 38, 3098. (41) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (42) Hehre, W.J.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. (43) Tosoni, S.; Tuma, C.; Sauer, J.; Civalleri, B.; Ugliengoa, P. J. Chem. Phys. 2007, 127, 154102. (44) Ramirez-Cuesta, A. J. Comput. Phys. Commun. 2004, 157, 226. (45) Alemagna, A.; Lorenzelli, V. J. Chim. Phys. Phys.-Chim. Biol. 1964, 61, 884. (46) Coclers, L.; Vanclef, A.; Bouche, R. Spectrosc. Lett. 1978, 11, 799. (47) Elliott, T. H.; Natarajan, P. N. J. Pharm. Pharmacol. 1967, 19, 209. (48) Kahovec, L.; Wagner, J. Z. Phys. Chem. 1941, B49, 156. (49) Hudson, M. R.; Allis, D. G.; Ouellette, W.; Hakey, P. M.; Hudson, B. S. J. Mol. Struct. 2009, 934, 138. (50) Hudson, M. R.; Allis, D. G.; Hudson, B. S. Chem. Phys. Lett. 2009, 473, 81. (51) Verdal, N.; Rivera, S. A.; Hudson, B. S. Chem. Phys. Lett. 2007, 437, 38. (52) Johnson, M. R.; Trommsdorff, H. P. Chem. Phys. 2009, 355, 118. (53) Guo, H.; Karplus, M. J. Phys. Chem. 1992, 96, 7273.

JP9114095