ARTICLE pubs.acs.org/crystal
Application of London-Type Dispersion Corrections in Solid-State Density Functional Theory for Predicting the Temperature-Dependence of Crystal Structures and Terahertz Spectra Matthew D. King and Timothy M. Korter* Department of Chemistry, Syracuse University, 1-014 Center for Science and Technology, Syracuse, New York, 13244-4100 ABSTRACT: Solid-state density functional theory (DFT) has been shown to be a valuable tool for the simulation of low-frequency vibrational motions in molecular crystals. While it is typically required that an experimentally known crystal structure be used as the initial input for these types of calculations, it is sometimes found that suitable crystallographic data are not easily obtainable. In this study, the low-temperature unit cell structure of the β polymorph of deuterated oxalic acid dihydrate, for which a structure has only been reported at room temperature, was predicted using solid-state DFT augmented with a modified empirical correction for weak long-range dispersive interactions. The dispersion correction parameters were optimized against the known 100 K crystal structure of the R polymorph, and then used for full-geometry optimization of the β crystal structure. Using this predicted structure for the β polymorph, the observed cryogenic THz spectrum of a mixture of deuterated oxalic acid dihydrate polymorphs was simulated.
’ INTRODUCTION Understanding the lattice vibrations that give rise to absorptions in terahertz (THz) spectra is a difficult task. The application of solid-state density functional theory (DFT) methods has proven to be a reliable approach to simulate THz absorption spectra and to reveal the underlying physical nature of these lowenergy vibrational motions.1,2 A limitation of solid-state DFT in the simulation of THz spectra is the need for preexisting crystallographic data for initial atomic positions and unit cell dimensions at the desired temperature given the strong dependence of THz absorptions on the experimental temperature.3,4 However, situations may occasionally arise where the crystal structure may not be readily obtainable for a variety of reasons. In such circumstances, it is of great interest to develop reliable computational methodologies to predict suitable crystal structures at a desired temperature for the simulation of THz spectra. In this study, an empirical dispersion correction term is applied to solid-state DFT calculations to predict the low-temperature crystal structure of the β polymorph of deuterated oxalic acid dihydrate. While there is only one known crystal structure of protonated oxalic acid dihydrate (OXDH), an additional polymorph exists exclusively in the fully deuterated form (OXDH-d6).57 The β-OXDH polymorph crystal structure has only been reported at room temperature, which is surprising because oxalic acid dihydrate has been a commonly used benchmark for X-ray and neutron diffraction studies for its simple structure and high symmetry.7,8 The lack of data can be attributed to the very similar crystal structures of the two polymorphs (Figure 1), leading to complications in selectively growing the β polymorph and also the slow conversion of the metastable β form to the Rpolymorph.6 Because of the inability to obtain a crystal structure at the required temperature, a pairwise dispersion correction r 2011 American Chemical Society
proposed by Grimme9,10 was included in the DFT calculations and modified to obtain a simulated 100 K crystal structure of the β polymorph to be used in subsequent normal mode calculations. The results of this study demonstrate that including proper dispersion interactions in DFT methods is essential for the accurate prediction of unit cell structures and THz spectra, and that using arbitrarily contracted unit cell dimensions to mimic cryogenic structures is not a reliable alternative.
’ METHODS Experimental spectra were recorded using a home-built timedomain pulsed THz spectrometer based on an amplified Ti: Sapphire femtosecond laser system. ZnTe crystals were used for generation of THz radiation by optical rectification11 and detection by free-space electro-optic sampling.12 A detailed description of the THz spectrometer has been reported elsewhere.3 Oxalic acid (g99%) and deuterium oxide (99.9%) were purchased from Sigma Aldrich and used without further purification. Oxalic acid dihydrate was crystallized by evaporation of a saturated aqueous solution at room temperature. Deuteriumsubstituted oxalic acid dihydrate was prepared by recrystallization of oxalic acid dihydrate five successive times in the same manner from D2O solutions. The descriptions of sample preparation and the THz measurements of OXDH have been previously reported.13 All DFT calculations were performed using the CRYSTAL09 software package.14,15 Calculations were performed utilizing the B3LYP density functional16,17 with the atom-centered triple-ζ 6-311G(2d,2p) basis set.18 To meet the total energy convergence Received: February 14, 2011 Published: March 15, 2011 2006
dx.doi.org/10.1021/cg200211x | Cryst. Growth Des. 2011, 11, 2006–2010
Crystal Growth & Design
ARTICLE
Figure 1. Crystal packing arrangements of the (a) R polymorph viewed down the b-axis and the (b) β polymorph viewed down the c-axis.
criteria of ΔE < 108 hartree for geometry optimizations and ΔE < 1011 hartree for normal mode calculations, shrinking factors of 10 and 7 for R-OXDH and β-OXDH, respectively, were used to specify the sampling rates as a function of k points used for calculating the density matrix and the commensurate grid in reciprocal space.19,20 A total of 175 666 grid points were used in the R-OXDH calculations, and 173 317 points were used in the β-OXDH calculations. The radial and angular distributions of grid points were defined by a pruned (75, 974) integration grid. Truncation tolerances used for Coulomb and HF exchange integral series were set to 109, 109, 109, 109, and 1018 hartree. Frequencies of normal modes were calculated within the harmonic approximation by numerical differentiation of the analytical gradient of the potential energy with respect to atomic position.21 The IR intensities for normal modes were calculated from the dipole moment derivatives (dμ/dQ) determined by the Berry phase method of calculating Born charges as polarization differences between equilibrium and distorted geometries.22 An empirical energy correction for London-type dispersion forces proposed by Grimme9,10 was used in all DFT calculations. The correction term is of the form N1
Edisp ¼ s6
N
ij
C6 ∑ ∑ ∑ ij fdmp ðRij, g Þ i ¼ 1 j ¼ iþ1 g R 6
where the energy is the sum over all atom pairs and g lattice vectors, C6 is the dispersion coefficient for atom pair ij, Rij is the interatomic distance between atoms i and j, and s6 is the functional-dependent scaling factor. A dampening function fdmp(Rij,g) = (1 þ ed(Rij,g/Rvdw1))1 is included to avoid near singularities for small R values. Rvdw is the sum of the van der Waals radii, and d is the steepness of the dampening function. Values for d, C6, and Rvdw were taken from ref 10 with the Rvdw values scaled as suggested by Civellari et al.23
’ RESULTS AND DISCUSSION The THz spectra of OXDH and OXDH-d6 collected at 78 K are shown in Figure 2. The protonated OXDH produced a THz spectrum with two absorptions located at 48.2 and 66.1 cm1. Surprisingly, the THz spectrum of OXDH-d6 contained a third absorption at 67.6 cm1 that could not be correlated to any feature seen in the OXDH spectrum as the additional mode is at higher frequency. The absorptions in the OXDH-d6 corresponding to those modes observed in the OXDH spectrum were shifted to lower frequencies of 46.7 and 63.9 cm1, respectively, as expected with the increased mass upon deuterium substitution. The relative intensities of these peaks change dramatically with the intensity of the lowest frequency mode decreasing
Figure 2. Terahertz spectra (10 to 80 cm1) of OXDH-d6 and OXDH at 78 K.
by nearly a factor of 8. It was suspected that this additional feature was a result of the presence of a small quantity of the β polymorph in the deuterated sample. Our previously reported THz spectrum and DFT simulations of protonated R-OXDH show this system to be well reproduced by solid-state DFT calculations.13 To verify the presence of the β polymorph in the sample mixture, solid-state DFT was used to simulate the THz absorption spectra of both OXDH-d6 polymorphs. For these systems, however, the conventional method of using fixed unit cell dimensions in the solid-state DFT calculations could not be employed as the necessary crystallographic data does not exist for the elusive β polymorph. The shortage of available structural information was overcome by taking a predictive computational approach in calculating a suitable crystal structure for the simulation of the THz spectrum near the experimental temperature at which the THz measurements were obtained. An empirical correction term for dispersive interactions, which are typically underestimated in DFT methods, was optimized specifically for the OXDH crystal system and included in the final DFT calculations. Using a more accurate representation of the intermolecular forces in the OXDH crystal, the dispersion parameters optimized against the known R polymorph were applied to predict the crystal structure of the β polymorph at cryogenic temperatures. It may seem unnecessary to take such a complicated approach to simulate the spectrum of the β polymorph when the room temperature crystal structure is experimentally known. However, the lattice contraction upon cooling is rarely uniform along each crystallographic axis, and the resulting affect on vibrational frequencies from any change in dimensions cannot be assumed. 2007
dx.doi.org/10.1021/cg200211x |Cryst. Growth Des. 2011, 11, 2006–2010
Crystal Growth & Design
Figure 3. (a) Background-corrected THz spectrum of OXDH-d6 obtained at 78 K. Calculated modes for β-OXDH-d6 using (b) published room temperature unit cell dimensions8 and (c) the room temperature dimensions contracted by 1.46% along each axis. Note the differences in the intensity magnitude between the experimental and theoretical data.
To demonstrate this, the vibrational frequencies of the β polymorph were calculated using the published room temperature unit cell dimensions8 and also using an assumed unit cell with these same dimensions contracted by 1.46% along each axis (Figure 3). This magnitude represents the average contraction of the experimentally known R-OXDH unit cell dimensions between 293 and 100 K. It is clearly seen in Figure 3 that the simulations based on these structures do not reasonably represent the experimental spectrum. It was observed that the frequency of the calculated vibrational modes decreased as the unit cell was contracted, which is contrary to expectations. It is also noted that the calculated intensities are extraordinarily high in comparison to the experimental absorptions. These results show that this approach to mimicking low-temperature unit cell dimensions does not yield physically meaningful structures. A more rigorous approach to the prediction of unit cell dimensions, such as the incorporation of dispersion forces in the DFT methods, is necessary. The unit cell dimensions and initial atomic positions for the ROXDH-d6 polymorph at 100 K were used as a starting point for the DFT calculations with dispersion corrections.24 The global scaling factor, s6, of the dispersion equation (see Methods section), which scales the magnitude of the overall dispersion energy correction, was adjusted to minimize the error between the fully optimized unit cell dimensions and the experimentally determined dimensions. This modification to the dispersion term is necessary due to the overcompensation of dispersion interactions when using the proposed dispersion parameters in the full-geometry optimization of unit cell dimensions.23,25 To accurately reproduce the THz spectrum, it is essential that the unit cell dimensions are comparable to that of the true crystal structure at the temperature at which the experimental THz
ARTICLE
spectrum was obtained. Small changes in unit cell volumes can have a significant impact on the calculated frequencies and intensities. By optimizing the s6 factor for the known R-OXDH-d6 system at 100 K, the same parameters for the dispersion correction term can then be used to predict the crystal structure of β-OXDH-d6 at 100 K, which has never been reported, for subsequent use in normal mode calculations. An optimum s6 value of 0.78 was obtained for the R-OXDH-d6 system, which is considerably smaller than the s6 value of 1.05 proposed for application of the dispersion corrections using the B3LYP functional.10 The fully optimized unit cell dimensions obtained using this optimized s6 value were a = 6.19950 Å, b = 3.46966 Å, c = 11.89478 Å, and β = 106.2976°. The average difference between the fully optimized and the experimental unit cell dimensions was 0.01765 Å, leading to a discrepancy in unit cell volume of 1.6407 Å3. In the full-geometry optimization, slight contractions of the b- and c- axes and expansion along the a-axis were observed. This s6 value was then applied in the fullgeometry optimization of the β-OXDH-d6 to simulate the 100 K structure. The initial unit cell parameters of a = 10.021 Å, b = 5.052 Å, c = 5.148 Å, and β = 99.27°, and space group of P21/a, were taken from the published room-temperature X-ray structure.8 In the β-OXDH-d6 full-geometry optimization, the unit cell contracted from a room-temperature volume of 257.219 Å3 to a volume of 240.632 Å3, a change of 16.587 Å3, or about 6.4%. The final fully optimized unit cell dimensions were a = 9.6309 Å, b = 5.5702 Å, c = 4.6540 Å, and β = 105.463°. Large contractions occurred along the a- and c- axes, with an expansion along the b-axis. The observed volume change is reasonable when considering the comparable decrease in volume for the R-OXDH-d6 polymorph of 10.4 Å3 or 3.9%, between these temperatures.5,24 Therefore, it was concluded that this simulated crystal structure may well represent the actual 100 K crystal structure for the β polymorph. The fully optimized ROXDH-d6 structure and that of the β-OXDH-d6 polymorph mimicking the 100 K crystal structure were then used for the simulation of the THz spectrum. The simulated THz spectra for both the R and β polymorphs are shown in Figure 4. The simulated spectra are overlaid with the baseline-corrected experimental 78 K THz spectrum of the OXDH-d6 sample. The simulated R-OXDH-d6 spectrum consists of two predicted vibrational modes at 50.3 and 69.1 cm1. The frequencies of the two vibrational modes were systematically overestimated in the calculations, consistent with that previously reported for OXDH and other hydrated molecular crystals due to the anharmonicity of the lattice vibrations.4,13 To compensate for this overestimation, a single frequency scalar of 0.925 was applied to the calculated frequencies to obtain those shown in Figure 4. This same scalar value has been used in previous reports using the B3LYP functional to calculate the vibrational frequencies of hydrated molecular crystals.4 The simulated THz spectrum for R-OXDH-d6 shows excellent agreement with the two lowerenergy features of the experimental spectrum. The normal mode calculation for β-OXDH-d6 using the predicted crystal structure produced a single vibrational mode at higher energy than those for the R-OXDH-d6 polymorph (Figure 4). This is consistent with that observed in the experimental OXDH-d6 THz spectrum. The frequency of this mode was at 76.0 cm1 after applying the same frequency scalar as determined for the R-OXDH-d6 system. There is a considerable difference in the spectrum simulated using the predicted structure from that shown in Figure 3b obtained using the room 2008
dx.doi.org/10.1021/cg200211x |Cryst. Growth Des. 2011, 11, 2006–2010
Crystal Growth & Design
ARTICLE
The ability to predict crystal structures of known and unknown polymorphs is of great interest in crystal design and engineering, and would be of great value for the analysis of THz vibrational spectra of both pure and complex sample mixtures.
Figure 4. Solid-state DFT simulations for the R and β polymorphs of OXDH-d6 overlaid with the background-corrected 78 K experimental spectrum. A scalar of 0.925 was applied to the calculated vibrational frequencies of both polymorphs. The intensity of the β-OXDH-d6 is scaled to match the observed experimental intensity. Lorentzian lineshapes were convolved into calculated modes using an empirical fullwidth-half-maximum of 1.87 cm1.
temperature β-OXDH-d6 structure. This large change in frequency and intensity of the calculated mode of the simulated low-temperature structure is a product of the variations in unit cell dimensions in the full-geometry optimizations utilizing the dispersion corrections. The vibrational character of this mode is purely intermolecular and analogous to the R-OXDH mode at 50.3 cm1.11 The vibrational frequency of the calculated mode for β-OXDH-d6 is significantly overestimated, and the IR intensity of this mode is much larger than those of the R-OXDH modes with an intensity of 4.56 km mol1 for the β-OXDH-d6, compared to the lower intensities of 0.19 and 0.18 km mol1 calculated for the R-OXDH-d6 vibrational modes. This difference in spectral intensities between the two polymorphs supports the idea that the sample mixture is composed primarily of the more stable R polymorph, with a small amount of the β polymorph. In Figure 4, the intensity of the β-OXDH-d6 vibrational mode is scaled to match the observed spectral intensity. If the calculated IR intensities are correct, this would indicate that approximately 11% of the OXDH-d6 contained in the sample is in the β form. Polymorph populations were checked by powder X-ray diffraction measurements which revealed approximately 5% of the sample to be the β polymorph. This suggests a slight underestimation of the calculated intensity of the mode. The results presented here demonstrate the capability of DFT with dispersion corrections to be used as a predictive tool in determining unknown crystal dimensions, provided that some knowledge of the crystalline system is available as a foundation. Proper optimization and application of these corrections for dispersion forces may facilitate the accurate simulations of crystalline structures and THz spectra at desired experimental temperatures. It has been shown that the lattice vibrational frequencies and intensities are sensitive to small changes in unit cell dimensions, and it is therefore imperative that crystal structures representative of those at the experimental temperature be used for the accurate simulation of THz spectra.
’ CONCLUSIONS The THz spectra of an OXDH-d6 polymorphic mixture was measured experimentally at 78 K and simulated by solid-state DFT including an empirical correction for dispersion interactions. By optimizing the magnitude of the London-type dispersion interactions in the DFT calculations using the known crystal structure for the R polymorph near the experimental temperature, a suitable low-temperature crystal structure for the unknown β-OXDH-d6 polymorph could then be predicted for subsequent use in normal mode calculations. The experimental THz spectrum of OXDH-d6 was shown to be a combination of both the R and β polymorphs present at roughly a 10:1 ratio. The results demonstrate the importance of London-type dispersion corrections in solid-state DFT calculations and its potential utilization in the predicting the temperature-dependence of crystalline structures. ’ AUTHOR INFORMATION Corresponding Author
*E-mail:
[email protected].
’ ACKNOWLEDGMENT This work was funded by a grant from the National Science Foundation CAREER Program (CHE-0847405). The authors thank Xue Han (Galson Laboratories), Ewelina M. Witko (Syracuse University) and Wayne Ouellette (Naval Surface Warfare Center, Indian Head Division) for thoughtful insights, and Syracuse University for continued support. ’ REFERENCES (1) Allis, D. G.; Zeitler, J. A.; Taday, P. F.; Korter, T. M. Chem. Phys. Lett. 2008, 463, 84–89. (2) Hakey, P. M.; Allis, D. G.; Hudson, M. R.; Korter, T. M. IEEE Sens. J. 2010, 10, 478–484. (3) Hakey, P. M.; Allis, D. G.; Ouellette, W.; Korter, T. M. J. Phys. Chem. A 2009, 113, 5119–5127. (4) King, M. D.; Buchanan, W. D.; Korter, T. M. J. Phys. Chem. A 2010, 114, 9570–9578. (5) Delaplane, R. G.; Ibers, J. A. Acta Crystallogr., Sect. B 1969, 25, 2423–2437. (6) Ebisuzaki, Y.; Angel, S. M. J. Raman Spectrosc. 1981, 11, 306–311. (7) Coppens, P.; Sabine, T. M. Acta Crystallogr., Sect. B 1969, 25, 2442–2451. (8) Iwasaki, F. F.; Saito, Y. Acta Crystallogr. 1967, 23, 64–70. (9) Grimme, S. J. Comput. Chem. 2004, 25, 1463–1473. (10) Grimme, S. J. Comput. Chem. 2006, 27, 1787–1799. (11) Rice, A.; Jin, Y.; Ma, X. F.; Zhang, X. C.; Bliss, D.; Larkin, J.; Alexander, M. Appl. Phys. Lett. 1994, 64, 1324–1326. (12) Wu, Q.; Litz, M.; Zhang, X. C. Appl. Phys. Lett. 1996, 68, 2924–2926. (13) King, M. D.; Korter, T. M. J. Phys. Chem. A 2010, 114, 7127–7138. (14) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. Z. Kristallogr. 2005, 220, 571–573. (15) Dovesi, R.; Saunders, V. R.; Roetti, C.; Orlando, R.; ZicovichWilson, C. M.; Pascale, F.; Civalleri, B.; Doll, K.; Harrison, N. M.; Bush, 2009
dx.doi.org/10.1021/cg200211x |Cryst. Growth Des. 2011, 11, 2006–2010
Crystal Growth & Design
ARTICLE
I. J.; D’Arco, P.; Llunell, M.; CRYSTAL09 User's Manual, University of Torino: Torino, 2009. (16) Becke, A. D. J. Chem. Phys. 1993, 98, 1372–1377. (17) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B: Condens. Matter 1988, 37, 785–789. (18) Frisch, M. J.; Pople, J. A.; Binkley, J. S. J. Chem. Phys. 1984, 80, 3265–3269. (19) Gilat, G. J. J. Comput. Phys. 1972, 10, 432–465. (20) Monkhorst, H. J.; Pack, J. D. Phys. Rev. 1976, 13, 5188–5192. (21) Pascale, F.; Zicovich-Wilson, C. M.; Gejo, F. L.; Civalleri, B.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 888–897. (22) Dall’Olio, S.; Dovesi, R.; Resta, R. Phys. Rev. B: Condens. Matter 1997, 56, 10105–10114. (23) Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. CrystEngComm 2008, 10, 1693. (24) Lehmann, A.; Luger, P.; Lehmann, C. W.; Ibberson, R. M. Acta Crystallogr., Sect. B: Struct. Sci. 1994, B50, 344–348. (25) King, M. D.; Buchanan, W. D.; Korter, T. M. Phys. Chem. Chem. Phys. 2011, 13, 4250–4259.
2010
dx.doi.org/10.1021/cg200211x |Cryst. Growth Des. 2011, 11, 2006–2010
