J o u r n a l Na me The necessity of periodic boundary conditions for the accurate calculation of crystalline terahertz spectra †
Terahertz vibrational spectroscopy has emerged as a powerful spectroscopic technique, providing valuable information regarding long-range interactions -and associated collective dynamicsoccurring in solids. However, the terahertz sciences are relatively nascent, and there have been significant advances over the last several decades that have profoundly influenced the interpretation and assignment of experimental terahertz spectra. Specifically, because there do not exist any functional group or material-specific terahertz transitions, it is not possible to interpret experimental spectra without additional analysis, specifically, computational simulations. Over the years simulations utilizing periodic boundary conditions have proven to be most successful for reproducing experimental terahertz dynamics, due to the ability of the calculations to accurately take long-range forces into account. On the other hand, there are numerous reports in the literature that utilize gas phase cluster geometries, to varying levels of apparent success. This perspective will provide a concise introduction into the terahertz sciences, specifically terahertz spectroscopy, followed by an evaluation of gas phase and periodic simulations for the assignment of crystalline terahertz spectra, highlighting potential pitfalls and good practice for future endeavors.
Introduction<!>History of Terahertz Spectroscopy as a Tool for the Chemical Sciences<!>Terahertz Vibrational Spectroscopy<!>Computational assignment of Terahertz Spectra<!>Overview<!>Predicting Accurate Mode-Types in Anhydrous Oxalic Acid<!>The Role of Subtle Interactions on Terahertz Spectra<!>Large Flexible Molecules and Polymorphism<!>Misleading Physical Interpretation with Gas-Phase Simulations<!>Conclusions and Outlook<!>Sample Preparation<!>Terahertz Time-Domain Spectroscopy<!>Periodic Simulations<!>Gas Phase Simulations
<p>Terahertz vibrational spectroscopy is a powerful technique for studying intermolecular interactions in the condensed phase, having been demonstrated to be useful for a wide-array of applications. [1][2][3][4][5][6] At terahertz frequencies (0.1 -10 THz or 3 -333 cm −1 ) large-amplitude complex motions of entire molecules are excited and sampled, and thus terahertz spectroscopy is a highly-sensitive probe of weak, often non-covalent, forces. The unique combination of sampling weak forces while probing long-range collective dynamics has spurred the growth of two, often related, lines of study. In the case of the former, weak forces in the condensed phase are inextricably dictated by the bulk packing of the individual molecules, making terahertz methods highly-sensitive to variations in three-dimensional structures. This has made terahertz spectroscopy a powerful complement to traditional structural methods (e.g., X-ray diffraction) for the identification of different packing geometries, as observed in crystalline polymorphs, for example. [7][8][9][10] On the other hand, the complex and largeamplitude motions occurring at terahertz frequencies have recently been related to a wide-variety of bulk phenomena, for example phase transformations, 11,12 mechanical responses, 5,13,14 electron-phonon coupling, 15,16 and gas-capture in porous mate-rials. 17,18 Unfortunately, the factors that make terahertz spectroscopy a valuable analytical technique simultaneously lead to one of its biggest challenges -the interpretation of experimental results. Unlike mid-infrared (mid-IR) methods, such as Fourier transform infrared spectroscopy (FTIR) and mid-IR Raman scattering, there do not exist any functional-group specific transitions at terahertz frequencies. Each individual bulk-phase geometry has a unique terahertz spectral fingerprint, and even cases where the packing between two polymorphs is very similar, there is often a dramatic reconfiguration of the terahertz vibrational density of states when transitioning between the two structures. 19 Thus, in order to effectively interpret experimental terahertz spectral data, additional analyses are required.</p><p>Over the last decade, methods for assigning and interpreting terahertz spectra have advanced significantly, in line with the growing use of terahertz spectroscopy, providing significant insight into many condensed phase phenomena. 5,[20][21][22][23][24][25][26][27][28] For example, the determination of low-frequency vibrations in organic semiconductors, in tandem with experimental terahertz spectra, has enabled a more thorough understanding of the effect of lowfrequency vibrations on charge carrier mobilities in these materials. 15,16 However, such insight is entirely dependent upon the correct prediction of both the vibrational frequencies as well as the associated mode-types, which is fundamentally related to the ability of the simulations to reconstruct the potential energy hy-persurface to a level of accuracy that enables a proper determination of the weak forces responsible for terahertz dynamics. Over the past two decades there have been a number of advances that have dramatically improved the description of terahertz dynamics, 29 including the determination of IR intensities in periodic simulations, 30 the incorporation of dispersion forces, 31,32 more advanced basis sets, 33,34 force fields, 35 and density functionals, 36 and so on. But, by-far the most critical to accurate terahertz analyses has been the development of robust periodic boundary condition simulations.</p><p>For the analysis of crystalline solids, periodic boundary conditions have been demonstrated to successfully reproduce terahertz spectra with a high degree of accuracy, as the long-range forces present in the solid are effectively captured by the simulation. 29 However, due to the complexity and (assumed) increased computational cost of periodic boundary condition simulations, gas-phase calculations involving clusters of molecules have been employed in studies describing terahertz vibrations, with varying levels of apparent success. [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] It is important to note that, in general, such assessment is based solely on the position of the calculated transitions, with little focus on the actual motions involved or their accuracy. However, there is a lack of comparison in the literature between these methods, and therefore this work aims to assess how gas-phase simulations perform in regard to reproducing not only the experimental terahertz spectra, but also the prediction of the vibrational mode-types. This is of critical importance, as the correct modeling of terahertz spectra implies that the simulation is fully capturing all of the interatomic forces with a high-degree of accuracy. When this criterion is met, experimental terahertz spectra can be used to validate the simulation, enabling access to many related properties that might be difficult to measure experimentally, for example thermodynamic parameters, 28 or interatomic forces. 14 These conditions lead to a powerful methodology for describing dynamic processes with atomic-level precision.</p><!><p>Terahertz radiation has been used in various capacities for much of the twentieth century, although technological difficulties in its generation and detection made it an elusive laboratory technique. Consequently, the terahertz region of the electromagnetic spectrum was referred to as the "terahertz gap." 56 The terahertz sciences grew out of efforts to measure ultrafast electrical signals, which in turn evolved into the generation and detection of freespace terahertz radiation. 3 While initial pioneering work utilized boutique laser sources and photoconductive antennae, 57 the advent of improved ultrafast lasers led to an increase in research related to terahertz radiation, as the femotosecond optical pulses enabled many lab-based methods to both generate and detect terahertz radiation. 56 Such a discovery lowered the barrier to terahertz experiments, and the general application of terahertz technologies began to flourish. 58 Expanding upon the initially used photoconductive antennae, new techniques based on non-linear crystals (e.g. ZnTe and GaP), such as optical rectification 59,60 and free-space electrooptic sampling, 61 led to robust tools for performing terahertz experiments in traditional academic laboratories. 62 More recently, advances in terahertz technology have resulted in countless new generation and detection schemes, with some of the more popular methods including improved photoconductive antennae, spanning a number of new materials and shapes, 63 as well as nonlinear crystals, such as 4-dimethylamino-N-methylstilbazolium tosylate (DAST), 64 (4-hydroxystyryl)-5,5dimethylcyclohex-2-enylidene)malononitrile (OH1), 65 lithium niobate (LiNbO 3 ), 66 and lithium tantalate (LiTaO 3 ), 67 to name a few, as well as air-plasma generation and detection schemes. 68,69 These sources are resulting in terahertz spectrometers that are capable of covering an ultrabroad spectral bandwidth (> 3 THz) with higher conversion efficiencies compared to traditional sources. Finally, in recent years quantum cascade lasers (QCLs) have become an area of intense study, with advances enabling high-powered narrowband -and often tunable -terahertz sources. 70 Early experiments in the terahertz sciences spanned a range of disciplines and applications. Some of the earliest experiments using terahertz radiation involved the study of gas-phase rotational dynamics -effectively an extension of microwave spectroscopy to higher frequencies. 71 For example, an early significant study produced a definitive measurement of the rotational spectrum of water vapor in the range of 0.2-1.45 THz, extending the assignment beyond the bandwidth of typical microwave spectrometers. 72 Other notable achievements in high-resolution rotational spectroscopy using terahertz technology include the measurement of rotational transitions in molecular oxygen isotopologues, and highly-accurate measurements of the rotational transitions of the NH 2 anion, critical to the subsequent investigation of interstellar spectroscopic measurements made by the Herschel Space Observatory. [73][74][75] More recently, terahertz QCLs have been used to perform ultra-high resolution gas-phase spectroscopy, with a reported resolution of < 0.5 MHz (5 × 10 −7 THz), with some studies even suggesting that QCLs can simultaneously act as both the emitter and receiver module for high-resolution spectroscopy -enabling a compact spectrometer design that is suitable for space-based applications. 76,77 A second area of terahertz spectroscopy that was quickly recognized leveraged the time-domain nature of most terahertz spectrometers. Because terahertz radiation is generated using ultrafast laser pulses, time-resolved terahertz spectroscopy quickly became a powerful tool for investigating the dynamics of materials; in particular, semiconductors. Terahertz radiation is strongly attenuated by free charge-carriers, and thus optical pump-terahertz probe spectroscopy (OPTPS) became a powerful method for characterizing the charge carrier dynamics of semiconductors with sub-picosecond temporal resolution. 3 Over the years there have been countless insight into carrier dynamics using OPTPS, ranging from the determination of carrier lifetimes -and by-extension, mobilities [78][79][80] -to the quantification of electron-phonon coupling in perovskite crystals. 81,82 This method has shown growth and utility over the years, with new developments occurring to-date. For example, a recent report highlighted that the chargecarrier dynamics in metal-organic frameworks can be obtained simply by using a powdered sample and a tape cell, making the determination of such data straightforward from a samplepreparation standpoint. 83 Outside of these two illustrative spectroscopic examples, the terahertz sciences have evolved greatly over the years, yielding a rich and diverse subset of applications, owing to the unique combination of the time-domain nature of terahertz generation, the properties of terahertz waves, and the terahertz response properties of materials. For example, the transparency of many materials to terahertz radiation, coupled with pulsed generation and detection, has provided an opportunity for terahertz pulsed imaging, which has proven to be valuable for the analysis of pharmaceutical coatings, 84,85 security and defense, 86 and the preservation of fragile artwork. 87 Recently, terahertz radiation is also being exploited for its use in telecommunications, as the higher frequencies (compared to current wireless communication standards) result in greatly increased bandwidth and speeds -making terahertz the successor to the 5G standard. 88 But by-far, one of the most widely adopted uses of terahertz radiation is vibrational spectroscopy, and in particular, condensed phase spectroscopy.</p><!><p>While terahertz spectroscopy was first applied to gases, 72 by far the most widely adopted application of terahertz spectroscopy is in the condensed phase. [13][14][15]17,28 In this regard, terahertz vibrational spectroscopy can be thought of as an extension of the more common mid-IR vibrational spectroscopy (e.g., FTIR) to lower frequencies. While far-IR spectroscopy had been performed over the course of the 20th century, the advent of terahertz timedomain spectroscopy (THz-TDS) instruments around the turn of the century reduced the barrier to accessing this region of the electromagnetic spectrum, and its popularity grew accordingly.</p><p>All vibrations in molecular systems are governed by the same set of conditions and selection rules, with the frequency of a given vibration dictated by the forces and masses involved in the motion, as shown in the solution of the harmonic oscillator Schrödinger equation,</p><p>where E v is the energy of the v th vibrational level, k is the vibrational force constant, and µ is the reduced mass. Thus, while mid-IR techniques probe (relatively strong) interactions, for example covalent bonds, terahertz spectroscopy probes weaker interactions, for example torsional motions of polymer side chains or hindered (rigid-body) translational motions of entire molecules in the condensed phase, as illustrated in Figure 1.</p><p>Unfortunately, while gas-and liquid-phase materials exhibit terahertz vibrational motions, the fact that each individual molecule exists in a different configuration and overall chemi- cal environment leads to a wide-range of terahertz transitions, which often results in a broad absorption feature commonly referred to as the vibrational density of states. 89 While there exists a large body of work in this area, the lack of discrete terahertz absorption features in disordered media limits analyzing atomiclevel details from the experimental data without significant efforts. 90 On the other hand, crystalline systems present an ideal medium for performing THz-TDS experiments, as each individual unit cell contains identical molecular conformations and chemical environments, repeated over the entire sample, leading to welldefined vibrational transitions and motions (i.e. phonons) that occur at the same (or nearly the same) frequency for each unit cell.</p><p>The dependence of terahertz vibrational transitions on not only the molecular conformation, but also the long-range structure, has made terahertz vibrational spectroscopy an ideal tool for characterizing the dynamics, and by extension, the structures, of crystalline polymorphs. 91 While terahertz spectroscopy is, in itself, not a structural probe; terahertz spectroscopy can be considered a complement to X-ray diffraction. As such, terahertz spectroscopy becomes an immensely useful tool for the characterization crystalline materials, rivaling X-ray diffraction in some cases, and surpassing X-ray diffraction in acquisition time. 23,[25][26][27]91 But while he ability of terahertz vibrational spectroscopy to characterize weak forces and complex dynamics in crystals is one of its most significant assets, it simultaneously leads to its most unfortunate aspects, namely, difficulties in interpreting and assigning the experimental spectra. As mentioned, each individual solid exhibits a unique terahertz spectrum, governed by a complex combination of forces, which ultimately does not result in any functional/structural-motif-specific vibrations or trends, unlike in the mid-IR. For example, except for the relatively-rare case of isomorphic solids (crystals with a different chemical composition but the same overall structure), 92,93 it is not possible to compare terahertz spectra from two different materials and compare them using common mid-IR methods, such as spectral peak shifting, as each solid has a completely unique spectral fingerprint. Therefore, the assignment of terahertz vibrational spectra requires additional analysis, which is most often accomplished through the use of simulation methods.</p><!><p>Unlike the mid-IR, which have well defined functional-group specific transitions, terahertz vibrations are unique to each individual solid, necessitating computational assignment. Complicating matters further is the nature of terahertz vibrations, which involve long-range weak forces and motions involving large numbers of atoms that are often of mixed mode types. Thus, where mid-IR vibrations are generally governed by the strength of a covalent chemical bond and transition frequencies can be readily predicted based on the identification of the molecule alone, for terahertz analyses, the entire solid (i.e. periodicity) must be considered in order to accurately capture all of the forces that result in the observed vibrational transitions.</p><p>Therefore, the assignment of terahertz spectra is, in itself, often more complex than the corresponding experimental measurement. Low-frequency vibrational transitions arise from extremely weak long-range forces (k), in tandem with larger bodies (i.e. entire molecules) constituting the reduced masses (µ) involved in the oscillation (see Equation 1). Because the mode-types of terahertz vibrations are so unique to the studied material and depend on some of the weakest forces in the solid, a highly-robust theoretical method must be used to capture these interactions accurately. Over the years, solid-state density functional theory (DFT), which incorporates periodic boundary conditions, has proven to be an effective methodology, as it fully takes into account the long-range forces found in crystalline solids. 94 In addition, related methods, such as periodic molecular dynamics simulations, have also been shown to be suitable for the reproduction of experimentally observed terahertz spectra. 95,96 However, many recent studies have used ab initio gas-phase simulations to assign experimental condensed phase terahertz spectra, with varying levels of apparent success. [37][38][39][40][41][42][43][44][45][46][47][48][49][50][51][52][53][54][55] It is important to note that in almost all of these cases the discussed 'agreement' with experiment is largely based upon the position of calculated vibrational transitions, with little focus placed upon the accuracy of the predicted mode-types -critical for gaining deeper insight into the way terahertz motions influence bulk material properties, for example. Therefore, this work aims to asses how gas-phase simulations reproduce terahertz spectra, as well as predicting the corresponding vibrational mode-types.</p><!><p>Due to the sensitivity of terahertz spectra to weak, long-range, non-covalent forces, a series of materials exhibiting a variety of these intermolecular forces were chosen. The vibrational spectra of these crystals were then calculated using fully-periodic DFT simulations using the CRYSTAL17 software package, 97 as well as with gas-phase DFT simulations on clusters of various sizes using GAUSSIAN09. 98 The experimental and solid-state DFT-predicted terahertz spectra of the materials chosen for this study are shown in Figure 2. Specifically, oxalic acid -one of the simplest organic solids -was chosen based on its prevalent hydrogen bonding and relatively small size. Its dihydrate was also chosen to explore how the inclu-Fig. 2 Cryogenic experimentally observed spectra (black) and theoretical spectra predicted with the B3LYP-D3/def2-SVP method (blue) for anhydrous oxalic acid, oxalic acid dihydrate, α-lactose monohydrate, and βlactose. The frequencies of the theoretically produced spectra using the def2-SVP basis set have been scaled by a factor of 0.9, while the spectra predicted using Ahlrichs VTZP basis set has been scaled by 0.94 for oxalic acid dihydrate, while the theoretical anhydrous oxalic acid spectrum is unscaled.</p><p>sion of co-crystallized water molecules influences the simulation of the spectra, particularly in the gas-phase clusters. Finally, two polymorphs of lactose monohydrate were selected, as they represent common standards for terahertz spectroscopy 99 and exhibit a range of interactions, including significant London dispersion forces and the presence of flexible torsional angles.</p><p>Solid-state simulations often employ experimental crystal structures as input geometries, which are then optimized prior to the calculation of the vibrational dynamics. 29 This represents a rather straightforward methodology, as it leave no room for user input, assuming a well-behaved material with an accurate experimental crystal structure. On the other hand, the gas-phase cluster simulations require a choice by the user, who must decide on the cluster size, configuration, number of molecules, and so on. Thus, gas-phase simulations were performed using geometries directly extracted from experimental crystal structures, as well as clusters generated in differing manners, for example single-layered sheets, multi-layered stacks, and in the case of oxalic acid dihydrate, subtle changes to only the initial water molecule configuration while keeping the organic portion fixed.</p><p>In order to effectively compare the periodic and gas-phase results, the exact same theoretical parameters were used for all simulations. Specifically, the B3LYP density functional, 100 along with Fig. 3 Predicted anhydrous oxalic acid low frequency spectra (A) gas-phase monomer and solid-state spectral predictions, with calculated vibrational frequencies (B) theoretical spectra generated from gas-phase 16-subunit cluster (red) and solid-state (black) prediction with calculated vibrational frequencies (C) sheet and stack cluster constructions of tetramer and octamers, with a narrowed FWHM for clarity (D) theoretical spectra of higher ordered clusters derived from experimental X-ray measurements, with a narrowed FWHM for clarity (E) structures of oxalic acid crystalline unit cell (i), octamer stacked cluster (ii), octamer sheet cluster (iii), and 16-subunit sheet cluster (iv).</p><p>the D3 dispersion correction, 101,102 was used with the def2-SVP basis set. 103 While the overall agreement with the experimental terahertz spectra and the periodic simulations is good after global scaling of the theoretical spectra for all systems, it is possible to improve the accuracy of the vibrational simulations by using a larger triple-zeta basis set (Ahlrich's VTZP 104 ), as shown for oxalic acid and its dihydrate in Figure 2. The larger basis set results in a better agreement with the experiment without requiring scaling, however it is important to note that the improvement is largely due to the better overall description of the intermolecular forces, as evidenced by the errors in the calculated unit cell parameters for the two methods, with the double-zeta basis set producing larger errors (∼ 2.3%) compared to the triple-zeta basis set (∼ 0.79%), which is what necessitates a larger scale factor for the smaller basis set, as discussed in previous work. 94 Because the goal of this work is to compare theoretical methods, and because the smaller basis set only results in a larger overall global shifting of the spectra, coupled with the drastically increased computational cost of the triple-zeta basis set, the def2-SVP was utilized throughout this study.</p><!><p>The ability of simulations to capture the the weak and long-range forces found in simple crystalline solids is critical to the success of the method, and the subsequent interpretation of the experimental results. Thus, anhydrous oxalic acid, one of the simplest hydrogen-bonded organic solids, was selected as the first system for study. Anhydrous oxalic acid crystallizes in the orthorhombic Pbca space group, with 4 molecules in the unit cell, and lattice parameters of a = 6.49 Å, b = 6.06 Å, c = 7.80 Å, and α = β = γ = 90 • , derived from single crystal X-ray experiments conducted at 130 K. 105 The crystalline structure was initially optimized with the fully periodic CRYSTAL17 DFT software package using the def2-SVP/B3LYP-D3 method, [100][101][102][103]106,107 which resulted in an average error of the unit cell parameters of 2.30%, compared to reported single crystal X-ray diffraction data. 105 The results of the fully periodic vibrational simulation are shown in Figure 3a. The simulation predicts 9 IR-active modes below 250 cm −1 , and only two between 0 -150 cm −1 , occurring at 76.44 and 107.85 cm −1 . These two modes correspond to antisymmetric unidirectional molecular translations of adjacent planes along the c-axis, and small translations along the a-axis with torsional motions between the carboxyl groups (Figure 4a,b), respectively, which is in agreement with the previous assignment of King and Korter. 108 To first assess the ability of gas-phase simulations to accurately predict both vibrational transition frequencies and mode-types, a single molecule of oxalic acid underwent geometric optimization and frequency analysis with GAUSSIAN09 98 using the same computational parameters as the periodic boundary condition calculations, producing the spectrum shown in Figure 3a. This simulation yields a single IR active mode in the low-frequency region at 42.83 cm −1 , and ultimately fails to reproduce the ex- perimental spectrum. Furthermore, the primary motion associated with this mode involves intramolecular torsional motion between the carboxyl groups (Figure 3b). A similar motion is observed in the mode produced by the periodic simulation at 107.85 cm −1 (+60.29% deviation), which is blue-shifted, compared to the gas-phase simulation, due to the incorporation of intermolecular forces in the periodic simulation.</p><p>However, it is important to note that the lack of intermolecular forces results in two differing, yet related, consequences. Firstly, the calculated motion demonstrates no intermolecular character (i.e. translation of adjacent molecules), which is present in the mode at 107.85 cm −1 predicted by the periodic simulation. While the two mode-types exhibit similarities in the intramolecular dynamics, the overall motion of the modes contrast. Secondly, the lack of intermolecular forces significantly softens this mode, leading to the large deviation in the calculated transition frequency between the two methods. It becomes clear that an isolated molecule fails to describe experimental terahertz spectrum, as well as the encompassing depiction of the associated motions.</p><p>In order to incorporate the weak intermolecular forces and allow for intermolecular motions in the gas phase simulations, clusters of increasing size and configurations were prepared. In doing so, the minimum number of molecules required to accurately reproduce the low-frequency spectrum in the gas-phase, as well as associated mode-types, can effectively be assessed. Clusters were constructed in three manners: hydrogenbound dimers forming stacks (Figure 3e, ii) and planar sheets (Figure 3e, iii -iv), as well as geometries generated directly from the reported single crystal X-ray structure (Figure 3e, i). Clusters ranged in their total subunit count from two to 36 molecules, with the larger clusters containing 24-and 36-subunits directly extracted from the crystalline supercell (Figure 3e). Each model underwent geometric optimization with no constraints followed by vibrational simulations, with a sample of the spectra presented in Figure 3. Predicted spectra and utilized cluster geometries are provided in the ESI.</p><p>As the number of molecules in the gas-phase simulation increases, more low-frequency vibrations are predicted -an unsurprising result given the the direct relationship between the number of atoms and vibrational degrees of freedom. However, this is not to say the increased number of predicted terahertz vibrations presents agreement with the experimental spectrum. Both the stack-and sheet-based clusters generate a large number of vibrational frequencies, offering little agreement with the experimental spectrum (Figure 3c). In the case of the 16-subunit cluster (Figure 3b), a total of 106 vibrational transition frequencies are calculated in the range of 0 − 250 cm −1 , significantly more than the 6 modes predicted from the periodic calculation. Moreover, the clusters produce IR-active modes that vary greatly in frequency, intensity, and mode-type, with no apparent convergence of such results. This finding is consistent across the clusters generated directly from the crystalline structure, with 162 vibrational modes from the 24-subunit cluster, and 246 vibrational modes from the 36-subunit cluster (Figure 3d) within the 0 -250 cm −1 region.</p><p>In some systems, there is apparent agreement between the generated vibrational spectrum predicted from gas-phase clusters and the experimental spectrum or periodic simulation. For example, the 36-subunit cluster predicts an intense transition occurring at 162.15 cm −1 , seemingly reproducing absorption feature at 160.27 cm −1 predicted by the periodic simulation (Figure 3d). Though similar in vibrational frequency and intensity, the modetypes of these vibrations differ considerably (Figure 4c). Specifically, the motions predicted by the periodic simulation exhibit antisymmetric translations primarily along the b-axis with a small contribution of translation along the a-axis. The gas-phase simulation presents no such concerted motion, with the atomic motions differing in magnitude and direction across each molecule in the cluster. Such a result indicates that while gas-phase simulations may produce vibrational frequencies falling within a close range of modes predicted from periodic simulations or experimental spectra, this does not mean the simulation is physically meaningful, and any apparent agreement exists only coinciden-Vibrational modes predicted from gas-phase clusters depict more intermolecular character as more molecules are introduced, but fail to do so on the same scale as the motions predicted from periodic boundary conditions (Figure 4a,b). Dynamics predicted from gas-phase simulations yield motions that differ in magnitude and direction throughout the cluster (Figure 4c), significantly deviating from the concerted motions expected of phonons, which are observed in the periodic simulations (Figure 4). In addition to inadequate representations of terahertz vibrations, the 36-subunit cluster was significantly more computationally costly, wherein 3N displacements are required to generate the hessian matrix from which normal modes are calculated following diagonalization, where N is the total number of atoms included in the simulation. Solid-state density functional theory exploits the symmetry of the crystalline system, wherein 3N displacements are required for the same process, but N is the number of atoms included only in the symmetry-independent unit of the primitive cell.</p><!><p>The study of anhydrous oxalic acid demonstrates the sensitivity of low-frequency vibrational frequencies and mode-types to both the correct representation of the forces, as well as the initial threedimensional structure. In order to further study the sensitivity of low-frequency spectra to these factors, oxalic acid dihydrate was investigated. To determine the effect of more subtle changes to the starting geometry upon the resultant terahertz spectra, the positioning of four water molecules was varied around a hydrogenbonded dimer of oxalic acid, while the geometry of the dimer was held constant (Figure 5a, i -iii). Even with only subtle changes in the starting positions of water molecules, the predicted theoretical terahertz spectra vary significantly, with each dimer cluster producing a unique spectrum, all of which fail to reproduce the experimental spectrum, as well as the spectrum predicted using the periodic simulation, indicating that even slight changes in the geometry are able to directly and strongly influence lowfrequency spectral predictions. This finding is of practical importance, as clearly the initial user-generated geometry significantly influences the final simulated low-frequency spectra.</p><p>With the dimer gas-phase simulations offering little agreement with the experimental spectrum, larger clusters of oxalic acid dihydrate were constructed in the same manner as the anhydrate, producing the spectra in Figure 5b. As was the case of the anhydrate, larger clusters (Figure 5A, iv -v) offer little agreement with both the experimental spectrum, as well as other cluster generated spectra (Figure 5c,d). Much like the anhydrate, some spectral features seem to resemble those of the experimental spectrum, but the corresponding mode-types fail to depict the concerted motion of molecules expected of phonons displayed in the vibrations of the solid-state simulations (ESI), and fall close to experimental frequencies only coincidentally.</p><!><p>Oxalic acid in both its anhydrous and dihydrated forms are small molecules, and the observed intramolecular contributions to lowfrequency vibrations primarily exist as torsions and rotations at higher vibrational frequencies (>100 cm −1 ). In larger molecules, such intramolecular character is significantly less hindered as conformational changes are more easily accessed, and associated motions occur at lower frequencies. 109 Such a consideration is exemplified in molecules able to adopt multiple molecular conformations, such as α-lactose monohydrate and anhydrous β -lactose (Figure 6a). It is possible that gas-phase simulations are able to more accurately predict low-frequency vibrational frequencies and motions composed of more intramolecular character, and thus the evaluation of gas-phase models was extended to both of these solids.</p><p>The periodic simulations for both α-lactose monohydrate and β -lactose are in good agreement with reported terahertz spectra. 110,111 To retain the anomeric differences and consequent long-range forces of the two solids in the gas phase simulations, clusters of lactose in both forms (Figure 6a) were constructed from a supercell constructed from single crystal X-ray experiments, from a single subunit to an octamer. 112,113 Following the trend of oxalic acid, the spectra produced from the gas-phase simulations of both lactose solids (Figure 6b,c) offer little agreement to the experimental spectra, with an excess of vibrational frequencies, absent in the experimental spectrum. Additionally, the mode types present non-uniform motions throughout the cluster, differing in both magnitude and direction (ESI), mimicking the shortcomings observed in the gas-phase simulations of oxalic acid and its dihydrate.</p><!><p>The proper assignment of terahertz spectra can have practical importance, and can strongly influence the interpretation of experimental results. A recent example highlights how improper analysis using gas-phase simulations can lead to incorrect conclusions. In a study by Ajito et al., the freeze-drying dynamics of aqueous NaCl solutions were investigated using THz-TDS, and over the course of the experiment the appearance, and subsequent disappearance, of two spectral features were observed. 39 The authors proposed a mechanism involving the formation of not only sodium chloride dihydrate (NaCl•2H 2 O) -an established reaction pathway -but also the formation of intermediate NaCl nanoparticles. The authors utilized gas phase cluster simulations based on molecular dynamics, and subsequently identified two structures that yield spectral features that are in some agreement with the experimental data. The conclusion of the observed NaCl nanopar-Fig. 7 Reported experimental (blue) terahertz spectrum of aqueous NaCl upon freezing (77 K), and the corresponding reported assignment of the spectrum to two different NaCl nanoclusters simulated using gas phase simulations (red and magenta). 39 The black curve is a periodic simulation of NaCl•2H 2 O, which is in excellent agreement with the experimental spectrum, representing a much more plausible assignment.</p><p>ticles was drawn primarily from these simulations. However, as discussed in this Perspective, it is trivial to obtain a spectrum that 'agrees' with an experimental result, by simply varying cluster geometries arbitrarily, and thus this is clearly not an appropriate methodology. In the case of this work, a fully periodic simulation of NaCl•2H 2 O results in a predicted spectrum that is in excellent agreement with the reported experimental spectrum, implying that authors instead observed the formation, and sequential dehydration, of bulk NaCl•2H 2 O crystals; a much more plausible, and fact-based conclusion than what is described in that work.</p><!><p>Across the solids surveyed in this study, the terahertz spectra generated from gas-phase clusters fail to replicate the simulations performed with periodic boundary conditions, which are in agreement with reported experimental terahertz spectra. Clusterbased simulations especially deviate from experimental spectra in the number of modes predicted in the terahertz region (0 -250 cm −1 ), in which differing motions are localized to specific molecules, rather than the concerted character detailed by vibrations occurring in the terahertz region. These discrepancies can be attributed to the the lack of incorporation of long-range structures and interactions, which are responsible for dictating the terahertz spectra of solids. Thus, while gas phase cluster simulations that seem to accurately capture the terahertz dynamics of solids can be obtained, they are not physically meaningful, as even subtle changes in the cluster configuration can drastically alter the predicted spectra. Over the last decade computational advances are enabling highly-accurate atomistic simulations of complex condensed phase materials to be performed. It is a great achievement that there exist robust tools, incorporating periodic boundary conditions, that can accurately reproduce the weak forces responsible for terahertz dynamics. While historically gas phase cluster simulations might have been the only reasonable option for simulating vibrational dynamics, this is no longer the case, and the community should strive to only recognize simulations that are of suffi-cient quality and rooted in physical phenomena for interpreting low-frequency vibrational spectroscopic data. Overall, the field of terahertz vibrational spectroscopy is more accessible than ever, and this understanding should become second-nature as this exciting field grows and matures.</p><!><p>Anhydrous oxalic acid, α-lactose monohydrate, β -lactose, were purchased and used as recieved. Oxalic acid dihydrate was obtained by recrystallization from an aqueous solution of oxalic acid, yielding large colorless needlelike crystals.</p><!><p>Samples for THz-TDS measurements were prepared by mixing with PTFE to a 10% w/w concentration, followed by homogenization via pestle and mortar. The powdered mixture was then pressed under 2 tons of pressure in a 13 mm diameter die, furnishing 3 mm thick freestanding pellets. All THz-TDS spectra were acquired using a commercial Teraflash spectrometer (Toptica Photonics AG). Cryogenic (20 K) spectra were collected using a closed-cycle liquid helium cryostat (Cryocool Industries). For sample, 30000 time-domain waveforms were collected and averaged, and the resultant waveform underwent Fourier transform to yield a terahertz power spectrum with a spectral resolution of < 0.1 cm −1 . The reported spectra represent the division of the sample power spectra by the PTFE blank spectra.</p><!><p>Solid-state density functional theory simulations, incorporating periodic boundary conditions, were performed using the CRYS-TAL17 DFT software package. 97 . The crystalline solids were modeled using the hybrid Becke 106 three-parameter, Lee-Yang-Parr 107 (B3LYP) 100 density functional, coupled with Grimme's DFT-D3 dispersion corrections, 101 including the Becke-Johnson damping function. 102 The atomic orbitals were represented using the splitvalence double-ζ def2-SVP basis set. 103 All structures were initially fully optimized (lattice parameters and atomic positions) with no constraints other than the space-group symmetry of the individual solids. In the case of anhydrous oxalic acid, geometry optimizations were also performed with Ahlrich's triple-ζ basis set (VTZP), 104 with three body dispersion contributions. This optimization offered little improvement over the def2-SVP basis set while significantly increasing computational cost, and thus def2-SVP/B3LYP-D3 was employed in all simulations. Upon complete optimization, vibrational analyses were performed, and infrared intensities were calculated using the Berry-Phase method. Energy convergence criteria were set to ∆E ≤ 10 −8 and 10 −11 hartree for the optimization and vibrational calculations, respectively.</p><!><p>Gas phase calculations were performed using the Gaussian09 DFT package 98 employing tight convergence criteria, enforcing an en-ergy criterion of ∆E ≤ 10 −13 hartree as well as an ultrafine integration grid. In order to most closely compare the gas phase and periodic simulations, all simulations performed in this manner employed the def2-SVP basis set, and the B3LYP density functional, with Grimme-D3 dispersion corrections.</p>