Prediction of the Unknown Crystal Structure of Creatine Using Fully

ARTICLE pubs.acs.org/crystal

Prediction of the Unknown Crystal Structure of Creatine Using Fully Quantum Mechanical Methods Matthew D. King,† Thomas N. Blanton,‡ Scott T. Misture,§ and Timothy M. Korter*,† †

Department of Chemistry, Syracuse University, 1-014 Center for Science and Technology, 111 College Place, Syracuse, New York 13244-4100, United States ‡ Eastman Kodak Company, Research Laboratories, 66 Eastman Avenue, Rochester, New York 14650-2106, United States § Department of Materials Science & Engineering, Alfred University, New York State College of Ceramics, Alfred, New York 14802, United States

bS Supporting Information ABSTRACT: Presented in this study is the full ab initio prediction of a previously unknown organic crystal structure. For the accurate prediction of crystal structures, the weak intermolecular forces must be well-represented. Solid-state density functional theory corrected for dispersion forces (DFT-D) is perhaps the most effective method for providing accurate descriptions of such forces in periodic systems. The success of DFT-D in crystal structure prediction is demonstrated in the prediction of the crystal structure of creatine, which has never been experimentally determined. The proposed P21/c structure was unquestionably verified by powder X-ray diffraction and terahertz spectroscopy, demonstrating a high degree of accuracy for the predicted structure. This clearly demonstrates the capacity and the feasibility of using full quantum mechanical methods in crystal structure prediction.

1. INTRODUCTION Steps toward the long sought after universal solution to the challenge of crystal structure prediction have involved many approaches. Molecular mechanics (MM),1
5 ab initio evolutionary techniques6
9 and hybrid MM/density functional theory (DFT)10,11 methods are among some of the most common for organic compounds. While each of these methods has demonstrated some degree of success, there remains much ground to be covered to reach adequate predictive capabilities for applications toward truly unknown crystal structures. Many of the limitations of current crystal structure prediction methodologies are rooted in the basic foundations of the approaches set forth to address the problem.12
15 Methods typically involve the exhaustive search of possible crystalline configurations on an intricate multidimensional potential energy surface; however, in doing so, the level of computational accuracy must be diminished to even attempt this approach within a reasonable time constraint. While these surveys may be beneficial as a first approximation, a more comprehensive approach is needed to address the subtleties of the complex interactions present in crystalline structures. Common deficiencies of most computational approaches to crystal structure prediction include the neglect or marginalization of weak intermolecular forces, such as London-type dispersion forces, and the use of rigid or semirigid molecular structures.15
20 When simulating crystal structures and searching for the global energetic minimum, one can ill afford to place limitations on the r 2011 American Chemical Society

accuracy of the calculations. Rigid-body models are typically optimized in the gas phase and are insufficient for solid state calculations. Molecules with any degree of conformational flexibility will vary in structure between the gas and condensed phases, particularly zwitterionic species which are intrinsically unstable in the gas phase. Furthermore, the use of rigid-body molecules in solid-state calculations is undesirable, as there will always be a certain amount of internal strain energy contained within the evaluation of the total lattice energy. Conformational energies of solid-state molecules are typically less than 8 kJ mol
1 but have been reported in excess of 30 kJ mol
1.4
11,15
26 This range far exceeds the energy distribution of proposed structures presented in nearly all crystal prediction studies. Identical molecules in different chemical environments may contain considerable variations in internal energy and charge distribution; therefore, it is not appropriate for any comprehensive method for crystal structure prediction to apply models that either limit conformational flexibility or generalize intermolecular potentials. Presented in this study is a full ab initio prediction of the previously unknown crystal structure of the zwitterionic molecule creatine using solid-state density functional theory augmented with corrections for dispersion forces (DFT-D). This rigorous Received: October 11, 2011 Published: November 16, 2011 5733

dx.doi.org/10.1021/cg2013599 | Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design

Figure 1. Molecular structure of the creatine zwitterion.

DFT-D method has been shown to indiscriminately reproduce structures with high accuracy over diverse crystalline chemical compositions because of its ability to offer system-specific representations of conformations and interaction energies. To no surprise, this has inevitably led to the exploration into the full potential for DFT-D as a method for crystal structure prediction.10,11,25 This study is the first account, to our knowledge, of the prediction of an unknown crystal system using a complete DFT-D approach. The crystal structure of creatine has never been solved by any crystallographic techniques. The absence of structural details for a molecule of such great biological and commercial value is surprising. In the body, creatine supplies energy for cellular functions by increasing available adenosine triphosphate (ATP).27 Because of this, creatine is widely used as a dietary supplement in the monohydrate and anhydrous crystal forms. Although the crystal structure had not yet been determined, the crystallinity of anhydrous creatine has been documented by powder X-ray diffraction (PXRD).28 However, no structural information could be retrieved from these measurements. The overall scarcity of structural information on creatine can be linked to its zwitterionic molecular structure (Figure 1). Creatine is insoluble in most pure solvents other than water, making direct crystallization of the anhydrous form so far unobtainable. Crystallization of creatine from aqueous solution exclusively forms the monohydrate structure, and attempts at crystallization through sublimation or annealing lead to the internal cyclization of creatine to form creatinine through a condensation reaction. Therefore, the only pathway to reaching the anhydrous crystal form is by dehydration of the monohydrate crystal. The loss of water through this dehydration process leads to a phase change producing microporous solids not suitable for single crystal X-ray measurements.28,29 However, it is evident by PXRD data that localized microcrystalline structures of anhydrous creatine do form upon dehydration of the monohydrate. A rigorous computational approach using solid-state DFT with periodic boundary conditions was implemented to predict the crystal structure of creatine in which the periodic DFT calculations were augmented with corrections for London-type dispersion forces of the form C6/R6.30
33 The inclusion of these forces is necessary to properly model the weak intermolecular interactions that play a critical role in determining molecular arrangements in the crystalline state. Initial structures for the calculations were selectively generated in the geometries most probable for organic molecules of this type. Since creatine is structurally similar to amino acids, close attention was given to the space group symmetries and hydrogen bonding motifs found in these types of solids. The amino acids containing polar side chains are primarily found in P-centered orientations with four molecules per unit cell and one molecule in the irreducible asymmetric unit.34
37 However, the basic amino acids with charged side chains, such as arginine and lysine, have no known crystal structures, owing to the same difficulties as with creatine concerning solubility and crystallization. By logically addressing the chemical nature of the creatine molecule prior to executing

ARTICLE

the crystal structure prediction procedure, a large number of “possible” structures were immediately eliminated, leaving a smaller pool of plausible structures to be taken into consideration. Limiting the 230 possible space groups to only a handful of chemically sensible ones greatly increases the likelihood of reaching the desired outcome of predicting an unknown crystal structure with high accuracy. This is not to say it would be impossible for the creatine molecule to exist in another space group, but it would be highly improbable based on the large body of literature focused on the structural details of similar organic compounds. In this study, the crystal structure of creatine was accurately predicted using solid-state DFT-D. The validity of the predicted structure and the theoretical model were verified using two experimental techniques: PXRD and terahertz (THz) spectroscopy. The proposed structure was initially confirmed by comparing the calculated and experimental PXRD patterns. Further structural evaluation was performed using THz spectroscopy and DFT-D normal mode vibrational calculations to investigate the low-frequency lattice vibrations of the creatine solid. An anomalous temperature-dependent shifting of vibrational frequencies was observed in the experimental THz spectrum, which was well predicted by the DFT-D spectral simulations. The theoretical treatment applied to the creatine system clearly demonstrates the reliability of this technique to reproduce physically meaningful crystal structures in a manner unparalleled by other available methods, and it sets a new benchmark for the field of crystal structure prediction. We have successfully demonstrated the feasibility of using high-level computational approaches for the prediction of unknown crystal structures. The unmatched sensitivity of DFT-D to locate minimum-energy configurations eliminates the dependence on precise initial structures which are required by less robust methods. Physically meaningful models of intermolecular interactions allow for the accurate representation of an unknown chemical system without any guesswork or indefinite parametrization of computational aspects. The result is the increased probability of uncovering the true structures of molecular crystals that have eluded traditional experimental and theoretical techniques.

2. METHODS 2.1. Theoretical. All DFT calculations with periodic boundary conditions were performed using the CRYTAL09 software package.30,31 Calculations used the PBE density functional with the atom-centered 6-31G(d,p) basis set.38,39 A total energy convergence criteria of ΔE < 10
8 hartree was used for geometry optimizations, and one of ΔE < 10
11 hartree was used for normal mode calculations. A shrinking factor of 6 was used in defining the sampling rate as a function of k points for the density matrix and commensurate grid in reciprocal space in order to meet the convergence criteria.40,41 The actual number of grid points varied between calculations depending on the initial dimensions of the input unit cell structures. The angular and radial distribution of points were defined by a pruned (75,974) integration grid. Truncation tolerances were set to 10
7, 10
7, 10
7, 10
7, and 10
14 hartree for the Coulomb and HF exchange integral series.42 Frequencies of normal modes were calculated within a harmonic approximation by numerical differentiation of the analytical gradient of the potential energy with respect to the atomic position.43 The IR intensities for normal modes were calculated from the dipole moment derivatives (dμ/dQ) determined by the Berry phase approach of calculating Born charges as polarization differences between equilibrium and distorted geometries.44 5734

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design

ARTICLE

An empirical energy correction for London-type dispersion forces proposed by Grimme32,33 was included in all DFT calculations (DFT-D). The dispersion correction is given by the equation Edisp ¼
s6

ij

∑i ∑j ∑g Rij,66g fdmp ðRij, g Þ C

ð1Þ

where the energy correction is the summation over all atom pairs and g lattice vectors. Cij6 is the dispersion coefficient, and Rij is the internuclear separation of atom pair ij, s6 is a global dispersion scalar, and fdmp is a dampening function included to avoid near-singularities at small values of R. The dampening function is of the form fdmp ðRij Þ ¼ ½1 þ e
dðRij =Rvdw
1Þ 
1

ð2Þ

where Rvdw is the sum of the van der Waals radii and d defines the function steepness. A cutoff distance for the dispersion energy term of 25 Å was used to truncate the summation over lattice vectors. Values for Cij6, Rvdw, and d were taken from ref 33, and the Rvdw values were scaled for the application of dispersion corrections in molecular crystals.33,45 The global dispersion scaling factor, s6, was modified to adjust for the magnitude of the dispersion forces applied to the calculations. It has been demonstrated in previous studies employing solid-state DFT-D that the proposed s6 values, originally optimized for molecules in the gas phase, overcompensate for dispersion forces in the reproduction of nonzero temperature crystal structures.46
48 Therefore, the proposed value of 0.75 for use with the PBE functional was reduced to 0.60 on the basis of the results of these previous studies as a reasonable starting point for the crystal structure predictions. Initial structures used in the survey of potential crystal structures were created by exploiting space group symmetries common to the structurally similar amino acids, which are typically found in space groups 4 (P21), 14 (P21/c, P21/a, and P21/n), and 19 (P212121). The preferential crystallization of biological molecules to limited space groups is well documented.49 The preliminary unit cell volumes were created from molecular “building blocks”. A molecule of arbitrary spatial configuration was constrained to a minimized rectangular volume. The unit cell volumes could be created from these building blocks by applying the appropriate symmetry operations at points placed at the corners of the asymmetric volumes or along a defined grid on the volume edges. Depending on the initial molecular orientation, this approach will most often produce unit cell configurations that are too large. The unit cell volume can be reduced by the symmetric reduction of the lattice vectors until reasonable intermolecular distances are achieved. The accurate representations of the intermolecular forces during the structural optimization procedures of the periodic DFT-D methods allow for such generalized initial structures. The constructed unit cells were fully optimized in which all lattice dimensions were permitted to float while maintaining space group symmetry. Structures that failed to optimize within the energy convergence criteria of ΔE < 10
8 hartree were disregarded as plausible structures. 2.2. Sample Preparation. Creatine monohydrate (Fluka 27900, g99%) was purchased from Sigma-Aldrich and used without further purification. Anhydrous creatine was obtained by dehydration of creatine monohydrate at 95 °C for several hours. Creatine samples for PXRD measurements were finely ground using a mortar and pestle. The powder was then passed through a 325 mesh nylon sieve and frontpacked into a glass holder. Samples for THz measurements were pulverized in a poly(tetrafluoroethylene) (PTFE) matrix using a stainless-steel ball mill (Dentsply Rinn 3110-3A) at a concentration of 2
6% by mass. The sample mixtures were pressed into 13  2 mm pellets under a pressure of 2000 psi using a hydraulic press (ICL EZ-Press 12) equipped with a 13 mm diameter stainless steel die. 2.3. Powder X-ray Diffraction. Powder X-ray diffraction (PXRD) measurements were obtained using a Bruker D8 Discover

diffractometer with the creatine sample mounted in a 1 mm capillary. The optics included an incident beam parabolic Goebel mirror, providing a parallel incident beam with the diffracted intensity measured using a scintillation detector. Data were collected from 2 to 70° 2θ with a step size of 0.03° 2θ for 30 s per step, using a Cu X-ray tube running at 1600 W. 2.4. Terahertz Spectroscopy. Experimental THz spectra were recorded using a time-domain pulsed THz spectrometer. A diodepumped solid-state Nd:YVO4 laser (Verdi-5, Coherent Inc.) was used to pump a Ti:sapphire oscillator (Mira-SPO, Coherent Inc.). The 800 nm output pulse seeded a Ti:sapphire regenerative amplifier (Legend-HE-USP, Positive Light), which was also pumped by a second diode-pumped solid-state Nd:YLF laser (Evolution-30, Positive Light). The 1 kHz pulsed output beam from the regenerative amplifier consisted of ∼800 nm pulses with a duration of ∼40 fs and energy of ∼2.5 mJ/ pulse. The near-infrared (NIR) pulses were propagated through a ZnTe semiconductor crystal generating ∼200 fs THz pulses with a bandwidth of 0.2
2.8 THz through the second-order nonlinear optical process of optical rectification.50 A PTFE disk located beyond the generator crystal absorbed any residual NIR radiation. The THz pulses were directed onto the sample by a pair of off-axis parabolic mirrors. The samples and blanks for measurement were held in a cryostat (VPF ST-100, Janis Research Company) equipped with 3-mm-thick polymethylpentene windows. Transmitted THz radiation was directed onto a ZnTe detector crystal by a second pair of off-axis parabolic mirrors. THz radiation was detected by free-space electro-optic sampling by measuring the polarization rotation of a NIR gate pulse overlapping spatially and temporally with the THz pulse at the detector crystal.51 The polarization of the NIR gate pulse was monitored by a balanced photoreceiver (2307, New Focus) connected to a lock-in amplifier (SR-830, Stanford Research Systems). The amplifier was locked to an optical chopper (C-980, Digirad) synchronized to a subharmonic of the regenerative amplifier repetition rate. A linear delay stage (ILS150CCHA, Newport) was used to vary the arrival time of the THz pulses relative to the NIR pulses at the detector crystal, yielding the time-domain THz waveform. Data were acquired at 293 and 78 K with samples under vacuum. The entire THz path was purged with dry air to eliminate absorption by atmospheric water vapor. Samples and blanks were scanned 32 times and data averaged for each individual data set. A 32 ps scan window consisting of 3200 data points was used to capture the THz waveform, which was then symmetrically zero-padded to a total of 6000 data points for the data transforms. The effective instrument resolution arising from the 32 ps scan length was approximately 1.0 cm
1. The ratio of the power spectra obtained from the Fourier-transformed data sets of the sample and blank yields the THz absorption spectrum. Each THz spectrum presented in this work is the average of four individual THz spectra, each representing a complete set of sample and blank measurements.

3. RESULTS AND DISCUSSION 3.1. Structural Optimizations. The majority of the structural optimizations from the selection of randomly generated initial crystal geometries (approximately 40 structures) from the aforementioned pool of probable structures failed to meet the energy convergence criteria imposed in the DFT-D optimization procedures. Many of the structures failed to arrive at stable energetic minima due to unfavorable spatial configurations resulting from the initial molecular orientations and imposed space group symmetries, thus preventing the rearrangement of the molecules to achieve the most stabilizing hydrogen bonding motifs. Structures unable to converge may also include those in local minima or in metastable phases. It is important to note that enforcing 5735

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design stringent convergence criteria prevents the optimization of structures in these less stable configurations, thus reducing the final number of “possible” structures that would otherwise be indistinguishable by less adaptive computational methods. In the creatine optimizations, many potentially suitable crystal structures were obtained with a range of unit cell energies and densities (Figure 2). Stable crystalline arrangements were preferentially found in the P21/c space group. Comparison of the energies of the optimized unit cells showed a single structure that was significantly lower in energy (∼11.8 kJ mol
1) than any other successfully optimized structure. Assuming a sufficient sampling of possible structures in the investigated space groups, the much greater stability of the lowest-energy structure would make this configuration heavily favored by thermodynamics. Locating a single structure of much lower energy quite likely gives the strong hydrogen bonding exhibited by the creatine molecule. The existence of additional kinetically favored polymorphs may be possible, although it is unlikely that these would be stable given the large disparity in energies between the optimized structures. The relatively small size and zwitterionic nature of the creatine molecule might also discourage the formation of kinetically favored polymorphs.52 On the basis of the results of the geometric optimizations, the lowest energy crystal structure of space group P21/c (Z = 4) became the focal point of further structural characterization.

Figure 2. Relative energies of calculated creatine crystal structures relative to the lowest-energy structure.

ARTICLE

To ensure that this structure was indeed at a minimum energy configuration in the P21/c space group, the complete unit cell was again optimized with P1 symmetry, providing the individual molecules unrestrained conformational freedom. The P1-optimized structure was determined to have retained P21/c symmetry with lattice dimensions nearly identical to that of the original P21/c unit cell. This verified that the structure was at a true minimum energy configuration under the constraints of the P21/ c space group symmetry. To this juncture, the geometry optimizations had employed an arbitrarily scaled dispersion correction term, yielding unit cell dimensions for the minimum-energy P21/ c structure of a = 9.7666 Å, b = 5.7831 Å, c = 11.5975 Å, and β = 105.937°. The unit cell contains four molecules, with one molecule in the irreducible asymmetric unit (Figure 3A). Each creatine molecule is stabilized by eight hydrogen bonds with four neighboring molecules (Figure 3B). This three-dimensional structure exhibits a highly favorable hydrogen bonding geometry between the charged portions of adjacent molecules to an extent not observed in any of the other obtained higher energy structures. A variable factor in DFT-D methods that must be taken into consideration is the magnitude of the applied dispersion correction term for the reproduction of nonzero temperature crystal structures. Because temperature contributions are not incorporated in DFT methodologies, the thermal expansion of unit cells may be compensated for through the scaling of the dispersion correction term. While this is not an ideal solution for the limitations of the DFT model, it does allow for the accurate prediction of temperature dependencies of molecular crystals. This approach has been previously demonstrated in the reproduction of known crystal structures at finite temperatures.46
48,53 One must be wary of the physicality of modifying the dispersion forces with the intent to compensate for thermal contributions, as too small or too aggressive corrections for such forces may compromise the meaningful representation of the true chemical nature of the system under investigation. With this in mind, the influence of the dispersion term scalar (s6 in eq 1; see Supporting Information) value on the lattice dimensions of the predicted P21/c creatine crystal structure was evaluated in the range 0.40
0.70 (proposed value of 0.75 for the PBE functional33,39). This arbitrary range was selected on the basis of existing literature values between 0.52 and 0.60 used to reproduce 100 K structures

Figure 3. (A) Representation of the predicted P21/c unit cell containing four creatine molecules. (B) Hydrogen bonding configuration of the creatine molecules in the crystal structure. 5736

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design

ARTICLE

Table 1. Unit Cell Parameters of Creatine Predicted by DFT-D and the Unit Cell Parameters Optimized To Best Fit the Calculated PXRD Spectrum with the Experimental PXRD Data PXRDa

DFT-D

differenceb

space group

P21/c

P21/c

Z

4

4

a (Å)

9.76678

9.865(6)

b (Å)

5.89251

5.869(4)

0.024

c (Å)

11.63125

11.699(8)


0.068

α (deg)

90

90

β (deg) γ (deg)

105.8705 90

105.89(4) 90


0.098


0.020

volume (A3)

643.8735

651.4


7.5

density (g cm
3)

1.352

1.337

0.015

a

Figure 4. Predicted unit cell dimensions of creatine as a function of the magnitude of the applied dispersion forces: (A) calculated unit cell volumes and (B) changes in lengths of lattice vectors with respect to the magnitude of the dispersion correction term scalar.

Figure 5. Comparison of the experimental and predicted powder X-ray diffraction patterns of creatine: (A) raw calculated PXRD pattern from predicted structure; (B) experimental PXRD pattern compared with that calculated from the crystal structrure obtained by a Rietveld refinement.

of molecular crystals using comparable DFT-D methods.46,48 It therefore seems reasonable to assume that an appropriate scalar for the creatine system lies within this range of values. The change in unit cell volume of the optimized structures varies linearly with the magnitude of the applied dispersion corrections (Figure 4A). As expected, the smallest unit cell volume was seen using the largest dispersion scalar value of 0.70 (i.e., strongest dispersion forces) and vice versa. The corresponding changes in lattice dimensions are also nearly linear in this s6 value range

Lattice parameters from powder X-ray diffraction data determined by optimization of the calculated DFT-D structure. b Deviation of calculated structure from that optimized to experimental PXRD data.

(Figure 4B). Little variation in length was observed along the a and c lattice vectors, with relatively large changes occurring along the b-axis. Between the s6 values of 0.40 and 0.50, the behavior of the b- and c-axes deviates slightly from linearity; however, this does not affect the linear trend of unit cell volume changes over the region. The differences in lattice dimensions from the lowest and highest volume structures were approximately 0.006, 0.159, and 0.042 Å along the a-, b-, and c- axes, respectively, resulting in an increase of 20 Å3 to the unit cell volume. An insignificant deviation of the β angle by 0.07° was also observed between these calculated structures (not shown in Figure 4). Using the structural data over a range of unit cell volumes, the predicted structure can be better compared to the appropriate experimental characterizations to resolve the accuracy of the crystal structure prediction. 3.2. Powder X-ray Diffraction. A distinct PXRD pattern was obtained for creatine (Figure 5). The use of capillary X-ray measurements eliminated any preferred crystallite orientation in the sample. Using the series of predicted P21/c crystal structures differing in volume, the predicted PXRD patterns were computed and compared to the experimental diffraction pattern. The raw simulated PXRD patterns all showed high agreement with the experimental PXRD pattern even before refinement of the predicted structures or crystal orientations (Figure 5A). Simple refinement of the lattice parameters yielded an extremely high correlation between the experimental and calculated PXRD patterns. The refined lattice parameters and the best-fit calculated structure agreed within 1% along all crystal axes (Table 1). These results are strong evidence that the predicted structure optimized using a dispersion term scalar of 0.40 represents the actual room temperature crystal structure of creatine. The lattice dimensions of this calculated P21/c structure were a = 9.7668 Å, b = 5.8925 Å, c = 11.6313 Å, and β = 105.871°. The predicted crystal structure was used as an initial input for the conventional Rietveld refinement. Very good agreement between the experimental PXRD and the calculated pattern from the refined structure was achieved, returning an Rwp value of 7.55% (Figure 5B). Given the large number of variables in the refinement procedure, only heavy atom positions were refined. Hydrogen positions have little influence on scattering intensities, making it difficult to determine their positions by any least-squares 5737

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design

ARTICLE

Table 2. Observed and Calculated IR-Active Vibrational Frequencies between 10 and 90 cm
1 for Creatine 293 K expa freqe

DFT-Db freq

intf

41.1

44.76

1.26

65.7

60.66 19.55 81.02

0.53

78 K exp freq

Δνd

DFT-Dc freq

int

41.8

45.51

1.84

64.4

57.07 17.97
1.3
3.59 78.43

1.09

exp 0.7

calc 0.75
2.59

a

Experimental. b Calculated using the predicted structure of Volume A. c Calculated using the predicted structure of Volume B. d Frequency differences between 293 and 78 K (experimental), and VA and VB (calculated). e Frequency (cm
1). f Intensity (km mol
1).

Figure 6. Experimental and simulated THz spectra of creatine. (A) Low-frequency THz spectra of creatine at 293 K (red) and 78 K (blue) at a concentration of 2% (w/w). Spectra in the inset taken at 6% (w/w) concentration. (B) Simulated THz spectra reproducing the anomalous behavior of the vibrational frequencies in response to changes in the unit cell volume of the predicted crystal structure. Lorentzian line shapes were convolved into the calculated spectrum using the best-fit line shapes of the experimental data.

minimization routine. However, the convergence of the Rietveld method and the high correlation in the diffraction patterns support the accuracy of the predicted creatine crystal structure. 3.3. Terahertz Spectroscopy. While the results of the PXRD measurements are convincing that the crystal structure of creatine was accurately predicted, the calculated structures were further evaluated by investigating the lattice vibrations, as these are sensitive to the intermolecular interactions of molecules in the solid state. Not only does performing normal mode calculations supply further insight into the accuracy of the predicted structure, they also help verify the stable minimum-energy confirmation of the crystalline structure. Nonstable structures, or those not at a true potential energy minimum, will normally yield imaginary frequencies in the normal mode calculations. Even with a stable structure, the accurate calculation of lattice vibrational frequencies requires highly accurate crystal structures and robust computational treatments because these low-energy vibrational modes are quite sensitive to minor deviations in atomic positions and weak intermolecular interactions. Therefore, reproduction of the low-frequency THz spectrum of creatine would be further evidence supporting the accuracy of the predicted crystal structure. The THz measurements of creatine revealed a simple spectrum with only two discernible absorptions between 10 and 90 cm
1 (Figure 6A). This is to be expected given the relatively small size of the creatine molecule and the high crystal symmetry, resulting in a low density of IR-active vibrations for such a system. The shifting of the vibrational frequencies of the two apparent

modes between the spectra collected at room temperature and 78 K is contrary to normal observations. Typically, vibrational frequencies in this spectral range increase as the measurement temperature is reduced due to contraction of unit cell volumes.46,48,54 Very few exceptions to this trend have ever been reported.54 In the creatine spectrum, the lower-frequency mode behaves as one might expect, increasing approximately 0.7 cm
1 between room temperature and 78 K (Table 2). The higherenergy vibrational mode, however, is red-shifted by about 1.3 cm
1 between these same temperatures. The origin of this anomalous shifting may be attributed to mechanical anharmonicity affects or to a phase change, although the latter is unlikely due to the preserved spectral features in the temperature data. The 293 and 78 K THz spectra were simulated by solid-state DFT-D using two separate unit cell volumes to examine the effect of unit cell size on the vibrational frequencies. The predicted room temperature structure had already been determined on the basis of the results of the PXRD data (structure now referred to as Volume A or VA). However, the experimental lattice dimensions of the 78 K structure were unknown, but for the purposes of this study, they can be approximated from the series of structures calculated using the modified dispersion corrections. For the reduced-temperature unit cell, the lattice dimensions of the structure calculated at the dispersion scalar s6 value of 0.55 were used (Volume B or VB). This s6 value is very similar to those used in the reproduction of crystal structures at 100 K in previous studies employing DFT-D methods.46,48 The lattice vectors of this unit cell were contracted from the room temperature structure along all axes, giving dimensions of a = 9.7662 Å, b = 5.8034 Å, c = 11.6033 Å, and β = 105.818°. The resulting reduction in volume of 1.73% is similar to that observed between room temperature and 100 K in strongly hydrogen bonded amino acids with known crystal structures, and this is thus a reasonable structural approximation for the difference in volume of the creatine unit cells.55,56 The spectral simulations using the VA and VB structures displayed the same mode-shifting trends as observed in the experimental THz spectra (Figure 6B). The lowest frequency mode was shifted toward higher energy upon reduction of the unit cell volume, while the high-intensity vibrational mode displayed the anomalous shift to lower energy. Since the vibrational frequencies are calculated within a harmonic approximation, this suggests that the shifting behavior is a result of the changes in unit cell dimensions and not some other phenomenon, such as anharmonicity or a solid
solid phase change. The exact reasoning for this unusual shift is difficult to uncover given the complexity of the hydrogen bonding network and the molecular 5738

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design motions associated with these lattice vibrational modes. However, the mode exhibiting the red-shifting involves the inversion of the
NH2 groups of creatine. The small differences in unit cell geometry and hydrogen bonding of the VA and VB structures may influence the inversion barrier and, therefore, the frequencies of the lattice vibrational motions in a manner contradictory to what is usually observed upon reduction of unit cell volume. The collective results of both the PXRD and THz analyses demonstrate the unequivocal accuracy and reliability of the applied DFT-D methods for the prediction of unknown crystal structures. Not only was the crystal structure of creatine predicted to a high level of precision, it was done so with no preliminary knowledge of the crystal system. In the process of crystal structure prediction, a meticulous computational treatment using DFT-D may, in some instances, be a more feasible approach in contrast to an exhaustive structural survey using generalized physical models lacking the ability to accurately model the intermolecular interactions based solely on the specific chemical nature of the system under investigation. Although the necessary time and computational resources for such a solid-state theory may far exceed other methodologies, the benefits of having well represented physical interactions outweigh these additional requirements. By applying the proper dispersion corrections to DFT calculations, the construction of precise initial crystal structures in the prediction method becomes of lesser importance to only the approximate orientations that molecules may exhibit in the true crystal structure. Relying on the meaningful physical representations of intermolecular forces and molecular structure found in DFT-D, as opposed to a comprehensive survey of possible structures, greatly increases the probability of success in crystal structure prediction.

4. CONCLUSIONS The previously unknown crystal structure of creatine was accurately predicted using solid-state DFT augmented with a correction term for London dispersion forces. The unit cell contains four molecules in the space group P21/c. The calculated lattice dimensions were a = 9.7668 Å, b = 5.8925 Å, c = 11.6313 Å, and β = 105.871°. Using the predicted dimensions, the structure was refined against the experimental PXRD pattern. The lattice dimensions of the refined P21/c crystal structure and those predicted were within 1% along all crystal axes. To further evaluate the predicted structure, the lattice vibrational modes of creatine were calculated and compared to THz measurements obtained at 293 and 78 K. The THz spectra showed an anomalous red-shifting of one vibrational mode upon cooling. This shift in vibrational frequency was reproduced by the DFT-D calculations between the room temperature structure and one of reduced volume, indicating the effect was a direct result of the contraction of the crystal lattice rather than originating in vibrational anharmonicity. The combination of the PXRD and THz data verified the accuracy of the predicted crystal structure of creatine, and it supports the practicality of using DFT-D for the purpose of crystal structure prediction. ’ ASSOCIATED CONTENT

bS

Supporting Information. Predicted crystal structures of creatine at the two unit cell volumes used in the frequency calculations. Volume A and Volume B structures have been provided with the file names CREATINE040 and CREATINE055,

ARTICLE

respectively. Additional Supporting Information includes the lattice parameters of the optimized crystal structures included in Figure 2. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This research was funded by a grant from the National Science Foundation CAREER Program (CHE-0847405). The authors also thank Syracuse University for continued support. ’ REFERENCES (1) Van Eijck, B. P.; Kroon, J. J. Comput. Chem. 1999, 20, 799–812. (2) Arikawa, T.; Tajima, N.; Tsuzuki, S.; Tanabe, K.; Hirano, T. THEOCHEM 1995, 339, 115–24. (3) Hirano, T.; Tsuzuki, S.; Tanabe, K.; Tajima, N. Chem. Lett. 1995, 1073–4. (4) Chaka, A. M.; Zaniewski, R.; Youngs, W.; Tessier, C.; Klopman, G. Acta Crystallogr., Sect. B: Struct. Sci. 1996, B52, 165–83. (5) Price, S. L.; Wibley, K. S. J. Phys. Chem. A 1997, 101, 2198–2206. (6) Cooper, T. G.; Jones, W.; Motherwell, W. D. S.; Day, G. M. CrystEngComm 2007, 9, 595–602. (7) Deaven, D. M.; Ho, K. M. Phys. Rev. Lett. 1995, 75, 288–91. (8) Oganov, A. R.; Lyakhov, A. O.; Valle, M. Acc. Chem. Res. 2011, 44, 227–237. (9) Woodley, S. M. Struct. Bonding (Berlin, Ger.) 2004, 110, 95–132. (10) Asmadi, A.; Neumann, M. A.; Kendrick, J.; Girard, P.; Perrin, M.-A.; Leusen, F. J. J. J. Phys. Chem. B 2009, 113, 16303–16313. (11) Chan, H. C. S.; Kendrick, J.; Leusen, F. J. J. Angew. Chem., Int. Ed. 2011, 50, 2979–2981. (12) Lommerse, J. P. M.; Motherwell, W. D. S.; Ammon, H. L.; Dunitz, J. D.; Gavezzotti, A.; Hofmann, D. W. M.; Leusen, F. J. J.; Mooij, W. T. M.; Price, S. L.; Schweizer, B.; Schmidt, M. U.; Van Eijck, B. P.; Verwer, P.; Williams, D. E. Acta Crystallogr., Sect. B: Struct. Sci. 2000, B56, 697–714. (13) Motherwell, W. D. S.; Ammon, H. L.; Dunitz, J. D.; Dzyabchenko, A.; Erk, P.; Gavezzotti, A.; Hofmann, D. W. M.; Leusen, F. J. J.; Lommerse, J. P. M.; Mooij, W. T. M.; Price, S. L.; Scheraga, H.; Schweizer, B.; Schmidt, M. U.; van Eijck, B. P.; Verwer, P.; Williams, D. E. Acta Crystallogr., Sect. B: Struct. Sci. 2002, B58, 647–661. (14) Day, G. M.; Motherwell, W. D. S.; Ammon, H. L.; Boerrigter, S. X. M.; Della Valle, R. G.; Venuti, E.; Dzyabchenko, A.; Dunitz, J. D.; Schweizer, B.; van Eijck, B. P.; Erk, P.; Facelli, J. C.; Bazterra, V. E.; Ferraro, M. B.; Hofmann, D. W. M.; Leusen, F. J. J.; Liang, C.; Pantelides, C. C.; Karamertzanis, P. G.; Price, S. L.; Lewis, T. C.; Nowell, H.; Torrisi, A.; Scheraga, H. A.; Arnautova, Y. A.; Schmidt, M. U.; Verwer, P. Acta Crystallogr., Sect. B: Struct. Sci. 2005, B61, 511–527. (15) Karamertzanis, P. G.; Pantelides, C. C. Mol. Phys. 2007, 105, 273–291. (16) Copley, R. C. B.; Barnett, S. A.; Karamertzanis, P. G.; Harris, K. D. M.; Kariuki, B. M.; Xu, M.; Nickels, E. A.; Lancaster, R. W.; Price, S. L. Cryst. Growth Des. 2008, 8, 3474–3481. (17) Oganov, A. R.; Glass, C. W. J. Chem. Phys. 2006, 124, 244704/ 1–244704/15. (18) Leusen, F. J. J. Cryst. Growth Des. 2003, 3, 189–192. (19) Gavezzotti, A. CrystEngComm 2003, 5, 439–446. (20) Price, S. L. Int. Rev. Phys. Chem. 2008, 27, 541–568. (21) Buttar, D.; Charlton, M. H.; Docherty, R.; Starbuck, J. J. Chem. Soc., Perkin Trans. 2 1998, 763–772. (22) Starbuck, J.; Docherty, R.; CharltonDavidButtar, M. H.; Buttar, D. J. Chem. Soc., Perkin Trans. 2 1999, 677–692. 5739

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740

Crystal Growth & Design

ARTICLE

(23) Yu, L.; Reutzel-Edens, S. M.; Mitchell, C. A. Org. Process Res. Dev. 2000, 4, 396–402. (24) Yu, L.; Stephenson, G. A.; Mitchell, C. A.; Bunnell, C. A.; Snorek, S. V.; Bowyer, J. J.; Borchardt, T. B.; Stowell, J. G.; Byrn, S. R. J. Am. Chem. Soc. 2000, 122, 585–591. (25) Asmadi, A.; Kendrick, J.; Leusen, F. J. J. Chem.—Eur. J. 2010, 16 (12701
12709), S12701/1–S12701/29. (26) Day, G. M.; Chisholm, J.; Shan, N.; Motherwell, W. D. S.; Jones, W. Cryst. Growth Des. 2004, 4, 1327–1340. (27) Wallimann, T.; Wyss, M.; Brdiczka, D.; Nicolay, K.; Eppenberger, H. M. Biochem. J. 1992, 281, 21–40. (28) Dash, A. K.; Mo, Y.; Pyne, A. J. Pharm. Sci. 2002, 91, 708–718. (29) Sakata, Y.; Shiraishi, S.; Otsuka, M. AAPS PharmSciTech 2005, 6, E527–35. (30) Dovesi, R.; Orlando, R.; Civalleri, B.; Roetti, C.; Saunders, V. R.; Zicovich-Wilson, C. M. Z. Kristallogr. 2005, 220, 571–573. (31) 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. Crystal09 User’s Manual; University of Torino: Torino, 2009. (32) Grimme, S. J. Comput. Chem. 2004, 25, 1463–1473. (33) Grimme, S. J. Comput. Chem. 2006, 27, 1787–1799. (34) Kistenmacher, T. J.; Rand, G. A.; Marsh, R. E. Acta Crystallogr., Sect. B 1974, B30, 2573–8. (35) Shoemaker, D. P.; Donohue, J.; Schomaker, V.; Corey, R. B. J. Am. Chem. Soc. 1950, 72, 2328–49. (36) Yamada, K.; Hashizume, D.; Shimizu, T.; Yokoyama, S. Acta Crystallogr., Sect. E: Struct. Rep. Online 2007, E63, o3802–o3803. (37) Suresh, S.; Padmanabhan, S.; Vijayan, M. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1996, C52, 1313–1316. (38) Krishnan, R.; Binkley, J. S.; Seeger, R.; Pople, J. A. J. Chem. Phys. 1980, 72, 650–4. (39) Perdew, J. P.; Burke, K.; Ernzerhof, M. Phys. Rev. Lett. 1996, 77, 3865–3868. (40) Gilat, G. J. Comput. Phys. 1972, 10, 432–465. (41) Monkhorst, H. J.; Pack, J. D. Phys. Rev. 1976, 13, 5188–5192. (42) Dovesi, R.; Pisani, C.; Roetti, C.; Saunders, V. R. Phys. Rev. B: Condens. Matter 1983, 28, 5781–5792. (43) Pascale, F.; Zicovich-Wilson, C. M.; Gejo, F. L.; Civalleri, B.; Orlando, R.; Dovesi, R. J. Comput. Chem. 2004, 25, 888–897. (44) Dall’Olio, S.; Dovesi, R.; Resta, R. Phys. Rev. B: Condens. Matter 1997, 56, 10105–10114. (45) Civalleri, B.; Zicovich-Wilson, C. M.; Valenzano, L.; Ugliengo, P. CrystEngComm 2008, 10, 1693. (46) King, M. D.; Buchanan, W. D.; Korter, T. M. Phys. Chem. Chem. Phys. 2011, 13, 4250–4259. (47) King, M. D.; Korter, T. M. Cryst. Growth Des. 2011, 11, 2006–2010. (48) King, M. D.; Ouellette, W.; Korter, T. M. J. Phys. Chem. A 2011, 115, 9467–9478. (49) Mighell, A. D.; Himes, V. L.; Rodgers, J. R. Acta Crystallogr., Sect. A: Found. Crystallogr. 1983, A39, 737–40. (50) Rice, A.; Jin, Y.; Ma, X. F.; Zhang, X. C.; Bliss, D.; Larkin, J.; Alexander, M. Appl. Phys. Lett. 1994, 64, 1324–6. (51) Wu, Q.; Litz, M.; Zhang, X. C. Appl. Phys. Lett. 1996, 68, 2924–2926. (52) Csaszar, A. G.; Perczel, A. Prog. Biophys. Mol. Biol. 1999, 71, 243–309. (53) Bond, A. D.; Solanko, K. A.; van de Streek, J.; Neumann, M. A. CrystEngComm 2011, 13, 1768–1771. (54) Walther, M.; Fischer, B. M.; Uhd Jepsen, P. Chem. Phys. 2003, 288, 261–268. (55) Boldyreva, E. V.; Kolesnik, E. N.; Drebushchak, T. N.; Ahsbahs, H.; Beukes, J. A.; Weber, H.-P. Z. Kristallogr. 2005, 220, 58–65. (56) Dittrich, B.; Huebschle, C. B.; Messerschmidt, M.; Kalinowski, R.; Girnt, D.; Luger, P. Acta Crystallogr., Sect. A: Found. Crystallogr. 2005, A61, 314–320. 5740

dx.doi.org/10.1021/cg2013599 |Cryst. Growth Des. 2011, 11, 5733–5740