Terahertz Vibrations and Hydrogen-bonded Networks in Crystals

The development of terahertz technology in the last few decades has made it possible to obtain a clear terahertz (THz) spectrum. THz vibrations clearly show the formation of weak bonds in crystals. The simultaneous progress in the code of first-principles calculations treating noncovalent interactions has established the position of THz spectroscopy as a powerful tool for detecting the weak bonding in crystals. In this review, we are going to introduce, briefly, the contribution of weak bonds in the construction of molecular crystals first, and then, we will review THz spectroscopy as a powerful tool for detecting the formation of weak bonds and will show the significant contribution of advanced computational codes in treating noncovalent interactions. From the second section, following the Introduction, to the seventh section, before the conclusions, we describe: (1) the crystal packing forces, the hydrogen-bonded networks and their contribution to the construction of organic crystals; (2) the THz vibrations observed in hydrogen-bonded molecules; (3) the computational methods for analyzing the THz vibrations of hydrogen-bonded molecules; (4) the dispersion correction and anharmonicity incorporated into the first-principles calculations and their effect on the peak assignment of the THz spectrum (5) the temperature dependence; and (6) the polarization dependence of the THz spectrum.


Introduction
Terahertz (THz) vibrations in the frequency range of 30-200 cm −1 give direct evidence for the formation of weak bonds in crystals. The THz region of the electromagnetic spectrum lies between the microwave and infrared (Figure 1) regions. The contributions and efforts in the field of materials research and crystal engineering [1][2][3][4][5][6][7], particularly those joined with the hydrogen bonding interactions in the field of THz spectroscopy, have resulted in a wide range of implementations of this method in practice in spectroscopic imaging [8], environmental inspection and detection [9], public security [10,11], biomedical analysis [12,13] and wireless communications [14][15][16][17]. The fundamental scientific contribution and wide extent of the THz-spectroscopic method's application results from the fact that THz radiation is transmitted through most non-metallic and non-polar mediums, thus enabling one to detect packaging, corrugated cardboard, clothes, shoes and book bags, in order to probe for potentially dangerous materials, in combination with minimal human health risk. The method enables the identification of explosives, chemical and biological agents through their characteristic fingerprint in the THz spectrum That weak intermolecular interactions in crystals provide a spectroscopic profile with strongly individual characteristics in the THz region, i.e., the THz method is very sensitive to the weak intermolecular interactions in crystals, enables the unambiguous identification of materials. Its applicability, including for the purposes of pharmaceutical analysis, in identifying biologically active samples has been demonstrated recently. This review is intended to give a coherent survey of the weak bond in assembling molecular crystals, of THz spectroscopy as a powerful tool for detecting the formation of weak bonds in crystals and of the significant contribution of advanced computational codes to analyze the THz spectrum. We especially focus on the stretching vibration of the weak hydrogen bond observed in the THz spectrum. Due to the recent development of THz vibrational spectroscopy, the spectral range of the vibrational absorptions of organic molecules in the solid phase obtained with high resolution nowadays covers the crystal-lattice and intermolecular vibrations, in addition to the intramolecular fundamental vibrations ( Figure 1). The applications of THz technology in probing the low-energy vibrational modes of materials have excited much interest [18][19][20][21][22][23][24][25][26]. Several recent breakthroughs made in light generation/detection have greatly contributed to the established THz spectroscopy [27][28][29][30][31][32][33][34][35][36]. Currently, the challenge and worldwide research activities around THz technologies have grown significantly with respect to the extension of the frequency range coverage and are related to interdisciplinary areas. The main efforts are focused on inorganic materials. Very recent findings in 2014 alone evidenced the perspectives of such studies [37][38][39]. Nevertheless, along with materials research and crystal engineering, i.e., the species interactions in crystals are focused mainly on inorganic materials, organic materials research also found a place in THz generation and detection. Such materials have shown great potential and broadband THz properties (10-1000 cm −1 ), extending the possibility of THz spectroscopy for the identification of materials in the mid-IR (infrared) region. Their application enables not only a qualitative fingerprinting analysis of materials, but also a quantitative one. Instrumental developments of THz spectrometers on the basis of organic materials have already been implemented. In this respect, the developments and contributions by Günter's group on nonlinear optical crystals have had a significant impact on the field of organic materials research and crystal engineering for applications in THz spectroscopy [40][41][42][43][44][45][46].
The analysis of the congested THz spectrum with vibrational absorption peaks appearing in the narrow wavenumber region of 30-200 cm −1 , has been very challenging. The rapid development of computer power, together with the availability of efficient codes, has made the assignment of observed absorption peaks increasingly reliable [47,48]. Intermolecular interactions and global nuclear motions in the crystalline environment are now being revealed through the assigned THz vibrational mode. It has been pointed out that the stretching vibration of the weak hydrogen bond is located in the THz region [49]. Among intermolecular interactions, the weak hydrogen bond is very important for many structures and functions in chemical and biological systems [49][50][51][52][53][54][55].
Sections 2-8 consist of the following contents. Section 2 elaborates upon the packing force, which works in a crystal growth process and governs the crystal structure of organic molecules. The hydrogen bond among the packing forces is explained in detail, and polymorphs with different crystal packing patterns of the same molecules are introduced. Section 3 shows the THz vibrations of hydrogen-bonded systems observed experimentally. Section 4 states the minimum requirements of computational methods for analyzing the THz spectrum of hydrogen-bonded molecules. Section 5 describes additional corrections to the density functional theory (DFT) method, such as dispersive interaction and anharmonicity. Both are crucial in the analysis of the THz spectrum. Sections 6 and 7 introduce the temperature and polarization dependence of the THz spectrum. Conclusions are given in Section 8.

Crystal Packing Forces and Hydrogen-Bonded Networks
The packing forces to assemble molecular crystals in the solid state are noncovalent, such as hydrogen bonding and van der Waals (vdW) interactions ( Figure 2). The -weak‖ intermolecular forces, like vdW or the hydrogen bond, play a significant role in a great variety of physical, chemical and biological phenomena. The effects of these forces manifest themselves, e.g., in the behavior of real gases, in the special properties of water or in the processes of protein folding. Consequently, a large number of the experimental and theoretical studies have been carried out in order to increase our knowledge about the nature of the intermolecular forces. It is interesting that different packing of the same molecules leads to the formation of a new polymorphic crystal structure. Hydrogen bonding was discovered almost 100 years ago [49][50][51][52][53][54][55], but still is a topic of vital scientific research. The reason for the long-standing interest lies in the eminent importance of hydrogen bonds for the structure, function and dynamics of a vast number of chemical systems, which range from inorganic to biological chemistry. The scientific branches involved are very diverse, and one may include mineralogy, materials science, general inorganic and organic chemistry, supermolecular chemistry, biochemistry, molecular medicine and pharmacy. In the hydrogen bond denoted by X-H···A, the group, X-H, is called proton donor or electron acceptor, and A is called the proton acceptor or electron donor. The chemical variation of the donor and/or acceptor, and possibly also of the environment, can gradually change the hydrogen bond to another interaction type. The transition to pure vdW interactions is very common. The hydrogen bond consists mainly of three interactions: electrostatic, induction and dispersion ( Figure 3). The strength of electrostatic interactions depends on the polarity of X δ− -H δ+ or A δ− . Induction is the interaction between a permanent dipole and an induced dipole. Dispersion is the interaction between induced dipoles. The polarity of X δ− -H δ+ or A δ− (or both) in the array, X-H···A, can be reduced by the suitable variation of X or A. This reduces the electrostatic part of the interaction, whereas the vdW component is much less affected. At the other end of the energy scale, there is a continuous transition to covalent bonding [56]. As a consequence, hydrogen bonds exist with a continuum of strengths. It can be useful for practical reasons to introduce a classification, such as -weak‖, -strong‖ and, possibly, also -in between‖. Jeffrey called hydrogen bonds moderate if they resemble those between water molecules or in carbohydrates (one could also call them -normal‖) and are associated with energies in the range of 4-15 kcal· mol −1 [57]. Hydrogen bonds with energies above and below this range are termed strong and weak, respectively. Weak and less common types of hydrogen bonds have been a major topic in hydrogen bond research. In the weak hydrogen bond of X-H···A, X and A are not simultaneously very electronegative atoms. The bond is very weak (ca. 1 kcal· mol −1 ), but it is important especially in the formation of molecular crystals or the higher order structure of proteins. The existence of the unconventional weak hydrogen bond was proposed early for the C-H···O contact, where a hydrogen atom is covalently bonded to a less electronegative carbon atom [58]. After years of controversy, the final crystallographic evidence for the existence of C-H···A (A = O, N, Cl) hydrogen bonds was provided by a survey of 113 neutron diffraction structures [59]. The THz spectra of crystalline solids are typically and uniquely sensitive to the molecular packing configurations, allowing for the detection of polymorphs. Polymorphism, the ability of solid material to exist in more than one form or crystal structure, is exceedingly important because of its relevance to structure-property relationships, the subtle and mysterious relationship between the structure and crystal packing of a molecule and the immense practical importance in the quality control of molecular materials industries, such as pharmaceuticals [60,61]. One is able to analyze and understand with increasing certainty a vast amount of accurate data currently available from crystallographic databases, while on the other hand, the appreciation of phenomena, like polymorphism, can lead to immediate benefits in applied areas, such as pharmaceutical development [62]. Many published works summarize the structural aspects of numerous polymorphic systems [63]. For distinguishing pharmaceuticals of a similar molecular composition, but differing crystal structures, THz spectroscopic and THz imaging studies on the polymorphism and/or cocrystal of several pharmaceuticals were intensively performed [64][65][66][67][68][69][70][71][72]. Here, spectral analysis with computational simulations is a key step for the distinction. A large number of molecules are found to exhibit polymorphism. The occurrence of polymorphic modifications in molecular compounds is manifested not only as a consequence of the minimum free energy in the crystalline phase, but also by the kinetics of the crystal nucleation and growth. Conformational polymorphs [73,74] are generated if the molecular conformation is different from one polymorph to the other; packing polymorphs [75] are those in which the molecules of a similar conformation pack into different crystal lattices and solvatomorphs (pseudopolymorphs) [76] that crystallize, including the solvents in the crystal lattice.
The polymorphs in which a weak hydrogen bond and/or other weak noncovalent interactions generate a distinct packing feature have been intensively investigated [77][78][79][80][81][82][83][84][85][86][87]. The second polymorph (Form II) of aspirin was recently found, and it is differentiated by the network of weak hydrogen bonds from the original form (Form I) [88]. The polymorphism of aspirin is, however, still an enigma. The molecular conformations are almost identical, showing only small differences of approximate two degrees in three dihedral angles, C 2 -C 1 -C 7 -O 1 , C 1 -C 2 -O 3 -C 8 and C 2 -O 3 -C 8 -C 9 ( Figure 4). The crystal structure of aspirin Form I was first determined by Wheatley in 1964 [89], but the gas-phase structure of aspirin was unveiled very recently by the rotational study of supersonically cooled aspirin, using laser ablation molecular beam Fourier transform microwave spectroscopy and its conformational analysis [90]. Aspirin Form II was predicted in 2004 [91] and observed experimentally the following year. Form II is relatively stable at 100 K, converting to the original Form I under room temperature. The crystal packing differs only slightly between the two polymorphs [92,93]. Aspirin Form I chains the carboxylic-acid and acetyl-group centrosymmetric dimers alternately in one dimension nearly along the a-axis (Figure 4a). An aspirin Form II dimer is connected via catemeric methyl C-H···O (Figure 4b). It was later found that the Form II structure is completely superimposable on the Form I structure. A very recent computational study revealed that Forms I and II are virtually degenerate [94]. That is, the crystal structure of Forms I and II is indistinguishable in structure and energy after geometry optimization. The detection of the bonding that distinguishes two polymorphs is needed to settle the dispute, and the stretching vibration of the weak hydrogen bond can be detected in the THz region. The crystallization of aspirin with weak hydrogen-bonded networks is easily affected by additives. The growth-induced twisting of aspirin crystallites grown from the melt is driven by salicylic acid, either generated from aspirin decomposition or added deliberately [95]. Banded spherulites of twisted aspirin crystals grew as either of the two known polymorphs of aspirin or as a mixture of the two. The spherulites are composed of helicoidal crystallites twisted along the <010> growth directions. The determination of the absolute sense of crystal growth in enantiopolar directions and the determination of the salicylic acid concentration variance in submicrometer-sized crystallites remain future experimental challenges.

THz Vibrations of Hydrogen-Bonded Molecules
THz spectroscopy [31,36] is synonymous with far-infrared (FIR) spectroscopy, in the sense that both cover the low-frequency region below 200 cm −1 . FIR spectroscopy first emerged in 1892 with the work of Rubens, but it took over seventy years to become widely available for scientific research [96,97]. The first commercial instruments were available in the 1960s, and extensive research in the vibrational spectroscopy of diverse organic molecules in the far-and mid-infrared frequencies proceeded over the next thirty years. Before 1990, there was a gap in the electromagnetic spectrum, 3-50 cm −1 , referred to as the -THz Gap‖, which was difficult to utilize, due to the lack of a suitable light source, and analysis using the frequencies between 50 and 500 cm −1 was usually referred to as FIR spectroscopy. The development of THz time-domain spectroscopy (THz-TDS) in the early 1990s was led by the initial discovery that THz pulses could be generated and detected [98,99], and THz spectroscopy stimulated a lot of early interest. In 1993, tunable THz laser vibration-rotation-tunneling spectroscopy emerged as a powerful tool for quantifying the structure and hydrogen bond rearrangement dynamics of small water clusters at an unprecedented level of accuracy [100]. In 2005, a THz-TDS study was reported on a thin film of liquid H 2 O, and a clear absorption mode was shown at 53 cm −1 , attributed to the bending motion of intermolecular hydrogen bond coordinated in a water cage [101]. The advantages of THz spectroscopy for analytical purposes are its nondestructive nature, the use of nonionizing radiation, the short acquisition times and the simplicity of sample preparation. The THz spectroscopy enables the characterization of solid-state materials through the excitation of soft intramolecular vibrational modes, as well as intermolecular modes and hydrogen-bonding networks inherent to the molecular assembly in the solid state [102,103].
The stretching vibration of the hydrogen bond in the THz region reflects the bonding itself [104][105][106][107]. The formation of a hydrogen bond affects the vibrational modes of the groups involved in several ways, and infrared (IR) spectroscopy is a standard method for investigating the hydrogen bond in the solid state [52]. In principle, H· · · A stretching vibration is the most direct spectroscopic indicator of hydrogen bonding. The frequency of stretching vibrations decreases as the interaction becomes weaker. The stretching vibration of a weaker kind of hydrogen bond appears in the THz region [49]. Reporting explicitly stating the stretching vibration of the weak hydrogen bond is very rare, despite the increase of publications on THz spectroscopic studies and their spectral analyses.
The O-H···O stretching vibrations of water clusters are in the THz region. The X-ray study of water indicates the tetrahedral structure of the hydrogen-bonded water pentamer [108,109]. The entire structure of five-molecule water is described exactly by the C 2 point group ( Figure 5). The water pentamer has 39 vibrational modes: 15 refer to intramolecular vibrations of H 2 O, symmetric stretching and bending and asymmetric stretching fundamentals; 15 refer to three restricted rotations or librations of individual H 2 O molecules; and the remaining nine refer to hydrogen-bond stretching and bending. Each vibrational mode may be mixed, but the description is given by the main motion. The librational modes of water were observed at 450-780 cm −1 [110][111][112] by FIR and Raman spectroscopy. The O-H···O stretching vibrations were observed at 167 (FIR) [110] and 175 (Raman) cm −1 [112] and the O-H···O bending vibrations at 53 (THz-TDS) [101] and 60 (Raman) cm −1 [112]. The IR spectrum of the water pentamer with C 2 symmetry calculated at Becke's three parameter exchange and Lee-Yang-Parr correlation functional (B3LYP) level using the correlation consistent polarized valence triple-zeta (cc-pVTZ) basis set [113,114] gives two strong peaks ( Figure 5). From the displacement vectors, the peak at the lower frequency (32 cm −1 ) among the two is assigned to the O-H···O bending vibration and the other peak at a higher frequency (174 cm −1 ) to the O-H···O stretching vibration.
The stretching vibration of intermolecular hydrogen bonds appears as intermolecular translational modes and provides evidence for the bonding, but the translational vibrations in crystals themselves are observed in the spectrum irrespective of the existence of a hydrogen bond. The translational vibration mixed with the stretching mode of intermolecular hydrogen bonds would, however, be higher in frequency than unmixed translational vibration. Both types of translational vibrations between one-dimensional chains, mixed and unmixed with the stretching mode of the weak hydrogen bond between chains, are expected in aspirin Form I in the THz spectral region [72]. Aspirin has four molecules in the unit cell, and the one-dimensional chain is constructed nearly in the direction of the a-axis (Figure 4a). In addition to some intramonomer modes, twelve intermonomer modes, six interdimer (or inter-chain) modes and three lattice modes exist in the THz region. Two translational vibrations nearly along the b-and c-axes are IR-active in the THz region ( Figure 6). The relatively weak peak at 38 cm −1 in the low-temperature THz spectrum is assigned to the translational vibration nearly along the b-axis. One of the three strong peaks between 62 and 78 cm −1 is assigned to the translational vibration nearly along the c-axis. The translational vibration along the b-axis in the present unit cell does not contain the stretching vibration of the weak hydrogen bond between chains, where the one-dimensional chain shifts mutually in the direction of the b-axis, as indicated by the arrows in Figure 6. The translational vibration along the c-axis contains the approaching/receding motion of methyl hydrogen to/from the oxygen atom of the adjacent chain. The distance between the hydrogen and oxygen atoms (2.690 Å) [92] from X-ray crystallographic data is nearly equal to a classical vdW separation (~2.7 Å) [115], supporting the existence of a weak hydrogen bond. The relatively strong peak assigned to the translational vibration in the direction of the c-axis is a candidate for the translation containing the stretching of the weak hydrogen bond. The frequency is higher than that of the other translational vibration containing no weak hydrogen bond stretching mode.
The inter-chain translational vibration mixed with the stretching of the hydrogen bond between the chains of poly-(R)-3-hydroxybutyrate (PHB) was observed near the translation of aspirin Form I containing the stretching of a weak hydrogen bond. PHB takes a lamellar helical structure in its crystal form [116,117]. The band at ~80 cm −1 in the temperature-and polarization-dependent THz spectra was perpendicularly polarized to the helical axis, suggesting its origin in intermolecular hydrogen bonds among the PHB chains [118]. On the other hand, weak peaks were observed around 40 cm −1 in the low-temperature THz spectrum of benzoic acid, and those are assigned to inter-layer translational modes unmixed with the stretching vibration of the hydrogen bond [119]. The frequency is near the unmixed inter-chain translational vibration of aspirin. Benzoic acid has a hydrogen atom instead of the acetoxy group of aspirin at the ortho position. The crystal has four molecules in the unit cell, and two of four form a centrosymmetric dimer. The dimers in the crystal are connected by weak hydrogen bonds and construct a layered structure in two dimensions. Two dimers in the unit cell are, however, not in the same layer, and no hydrogen bond is constructed between the dimers in different layers. The interdimer translational vibration, which is an inter-layer translation, contains no stretching vibration of a weak hydrogen bond. It turns out that the translational vibration mixed with the stretching vibration of the hydrogen bond, as in PHB, is higher in frequency than the unmixed translational vibrations, as in benzoic acid. Figure 6. The THz spectrum of aspirin at 4 K, and the displacement vector of two modes assigned to translational vibrations [72]. The arrows indicate the displacement of chains nearly along the b-axis (left) and the c-axis (right).
Low-frequency vibrational modes of protic ionic liquids (PIL) and ionic liquids (ILs) detect and quantify hydrogen bonds [120]. PILs are a subset of ILs and are formed by the combination of equimolar amounts of a Brønsted acid and a Brønsted base [121][122][123][124][125][126][127][128]. Hydrogen bonding in PILs is reminiscent of water [126]. The FIR spectra of dried PILs (ethylammonium nitrate, propylammonium nitrate and dimethylammonium nitrate) at 353 K give the asymmetric and symmetric stretching modes of the hydrogen bonds N-H···O around 199-224 cm −1 and around 134-159 cm −1 and the bending modes of the hydrogen bonds around 60-78 cm −1 . The measurement of the complex dielectric function of the PIL ethylammonium nitrate reveals a THz mode described as a damped harmonic oscillator with a central frequency of 43 cm −1 [129]. THz spectroscopy of ethylammonium nitrate indicates a fast process described by a slightly underdamped harmonic oscillator with a damping time constant of 2/γ = 830 fs, corresponding to the typical time scale expected for an intermolecular mode, due to hydrogen bond reformation. ILs are understood as liquids consisting entirely of ions and having melting points below 373 K. The presence of hydrogen bonding in the structure of an IL, 1-alkyl-3-methylimidazolium salt, was first reported by Seddon et al. in 1986 [130]. It is assumed that hydrogen bonding plays an important role in the properties and reaction dynamics of these Coulomb systems. The direct observation of hydrogen-bond stretching frequencies are reported from the FIR spectrum in pure imidazolium ILs containing the same anion, bis(trifluoromethylsulfonyl)imide (NTf 2 − ), but various cations: 1,2,3-trimethylimidazolium (1,2,3-trimethyl-im + ), 1,3-dimethylimidazolium (1,3-dimethyl-im + ), 1,2-dimethylimizazolium (1,2-dimethyl-im + ) and 1-methylimidazolium (1-methyl-im + ) [131]. The maximum intensities below 150 cm −1 are measured at 62.  [132]. Though it has been argued that significant vibrational contributions may occur between two and 50 cm −1 , the THz spectra of these ILs show no contribution blow 50 cm −1 [133]. The importance of hydrogen bonding in imidazolium-based ILs was the subject of highly controversial discussions. Because ILs solely exist as anions and cations, it is generally assumed that Coulomb interactions absolutely dominate the properties of these materials. Although the local and directional type of interaction contributes only one tenth of the overall interaction energy, hydrogen bonds have a significant influence on the structure of these Coulomb systems. With increasing hydrogen bond abilities and strengths, the NTf 2 − anion is found in the cis configuration.
Switching off these local interactions leads to the energetically favored trans conformation [134].

Computational Methods to Analyze the THz Spectra of Hydrogen-Bonded Molecules
The THz spectra of hydrogen-bonded molecules have a lot of peaks in the narrow spectral region of 30-200 cm −1 . The assignment of peaks requires a precise DFT or first-principles calculation at an accuracy of less than 10 cm −1 in frequency. Powdered samples are popularly used for THz measurements, and solid-state calculations are basically required. In some cases, the gas-phase calculations ensure the analysis. In this section, we indicate the least requirement of the calculation level for the analysis of the THz spectra of hydrogen-bonded molecules. We do not cover all the methods in detail. The calculation level corresponding to what is stated here is needed for the reliable analysis of the THz spectra. Additional crucial contributions to the molecular structure and THz spectrum of hydrogen-bonded molecules, such as dispersive interactions and anharmonic force fields, are described in the next section.
An ab initio molecular orbital method at the MP2 (second-order Møller-Plesset perturbation theory) level and DFT calculations by different functionals (B3LYP (Becke's three parameter exchange and Lee-Yang-Parr correlation functional), PBE (gradient-corrected correlation functional of Perdew, Burke and Ernzerhof) and MPWB1K (a combination of a modified Perdew-Wang exchange functional by Adamo and Barone and the B1K correlation functional formulated by Truhlar and co-workers) [135]) are very useful for assigning the vibrational absorption of gas-phase molecules. It is very well known that DFT, especially with hybrid functionals, such as B3LYP, allows for an accurate prediction of the properties of hydrogen-bonded complexes. A newly developed exchange-correlation functional, MPWB1K, is reliable for predicting intermolecular vibrations of doubly hydrogen-bonded complexes, such as DNA (deoxyribonucleic acid) base pairs [136,137]. Vibrational frequency analyses should be carried out at the fully optimized geometry and at the calculation level used in the geometry optimization. To obtain a reliable frequency value, larger basis sets, such as 6-311++G(d,p) (split-valence triple-zeta polarized basis set with diffuse functions) and cc-pVTZ, are recommended. Diffuse functions are more crucial than the use of large correlation consistent basis sets, like cc-pVTZ. Vibrational frequencies, for example, calculated at the B3LYP/6-311++G(d,p) level, produced better agreement with experimental mid-IR results than the calculated frequencies at the B3LYP/cc-pVTZ level [138,139]. The use of a correlation consistent basis set augmented with diffuse functions, such as aug (augmented)-cc-pVTZ, would give much better results, but great computer resources are required to complete the calculation. Frequency calculations at the B3LYP level using the 6-31+G(d) basis set gives sufficient results for large cluster systems, like ionic liquids [120,126,134]. Here, the anharmonicity described in the next section is corrected by the standard scaling factor, 0.96, and Grimme's dispersion-corrected method (DFT-D3), described in the next section, is recommended for the calculation of THz vibrations.
The total-energy calculation of the periodic system based on the plane-wave DFT method for the generalized gradient approximation (GGA) using several exchange-correlation functionals is used to investigate intermolecular vibrations in crystals. Extremely accurate calculations are necessary to clear up uncertainties with regard to the assignment of the complicated THz spectrum. The PBE functional is frequently selected for use, as the functional has been shown consistently to produce high-quality simulation of the solid-state structure and THz spectrum [140]. The kinetic energy cutoff, to which the valence electron wave function is expanded, is desired to be as large as possible. Electronic minimization needs to be converged to the limit of machine accuracy. Full geometry optimization for both the cell parameters and atomic coordinates is preferred. The total energy of the system needs to be converged to better than 0.1 meV/atom in terms of k-point sampling for the integration over the Brillouin zone.

Dispersive Interactions and Anharmonicity
The THz spectrum reflects a weak intermolecular interaction, and thus, the dispersive interactions and anharmonic force fields play a decisive role in the peak assignment. One of the critical issues in molecular-crystals simulation is the proper description of noncovalent interactions between molecules (i.e., hydrogen bonding, vdW interactions, dipolar effects, etc.) that dictate the crystalline structure and the thermodynamic properties controlling phase transitions. In particular, dispersion correction is important for weak hydrogen bond networks in crystals. It is generally accepted that the electrostatic and induction components are smaller in weak hydrogen bonds and that the dispersive interaction is relatively larger [49,50,141].
The high level of quantum chemical methods (e.g., a coupled cluster with single, double and triple perturbative excitations, CCSD(T)) can attain a consistently accurate description of vdW interactions. Unfortunately, CCSD(T) is limited to rather small systems, due to its high computational cost and steep O(N 7 ) scaling with system size. MP2 is the most economical wave-function-based electronic structure method beyond the Hartree-Fock (HF) approximation, which provides an approximate description of all relevant vdW interactions-electrostatics, induction and dispersion. One of the serious shortcomings of the MP2 theory is the noticeable overestimation of dispersion-interaction energy [142,143]. DFT is widely used in the computationally chemistry community, and B3LYP is the most popular density functional. However, one of the serious shortcomings of B3LYP is the inaccuracy of interactions dominated by medium-range correlation energy, such as vdW attraction. The Minnesota 05 (M05) and Minnesota 06 (M06) classes of functionals developed by the group of Truhlar at the University of Minnesota are new hybrid meta-GGA exchange-correlation functional [144][145][146] with broad applicability in chemistry. Based on the comparison of the performance among M06 class functionals and one M05 class functional (M05-2X) for diverse databases, M06-2X, M05-2X, M06-HF, M06 and M06-L are recommended for the study of noncovalent interactions [147]. M06-L is local functional with 0% HF exchange. M05, M05-2X, M06, M06-2X, M06-HF are global hybrid functionals with 26%-100% HF exchange. The structure and vibrational properties of the hydrogen-bonded molecules, complexes and clusters are studied with the use of M06 class functionals [148][149][150][151][152][153][154][155][156][157][158][159].
It is well known that DFT does not describe a noncovalent interaction, though DFT provides reasonably accurate predictions for many properties of various molecules and solids and is the most widely used method for electronic structure calculations in condensed matter physics and quantum chemistry [160][161][162]. The development of procedures for dispersion correction to DFT is an active area of research. Among the most popular methods are the DFT-D2 [163] and DFT-D3 [164] procedures by Grimme and co-workers. The dispersion-correction term is given by Equation (1): Here, R is an interatomic distance, f denotes a dumping function and C is a dispersion coefficient. In the DFT-D2 approach, dispersion corrections are introduced by adding the first damped atom-atom C 6 /R 6 term of Equation (1), and in the DFT-D3 approach, they are incorporated through the first two damped atom-atom C 6 /R 6 and C 8 /R 8 terms of Equation (1). These terms are smoothly cut off in the short range, where they are not relevant, but explicitly enforce the desired long-range asymptotic behavior. In the DFT-D2 method, the parameters are obtained by fitting to experimental C 6 coefficients and/or post-Hartree-Fock binding energy data and are independent of the electronic structures and chemical environments. This deficiency is partially remedied in DFT-D3, where the C 6 coefficients depend on the coordination number of atoms. Of the closely related approaches, the vdW-TS (Tkatchenko-Scheffler) approach [165,166] is particularly intriguing because of how the C 6 coefficients depend on the chemical environment. The vdW-TS approach is a parameter-free method for the accurate determination of long-range vdW interactions from mean-field electronic structure calculations and improves the description of weakly bonded systems in DFT for a range of exchange correlation functionals. The utility of the vdW-TS procedure has been explored for describing several crystalline materials, including noble gases, molecular crystals and layered materials, in which dispersion interactions are important [167].
The calculation of vibrational frequencies is performed at the minimum on the potential energy surface. Two alternatives in the geometry optimization process are open for including a dispersive interaction into solid-state DFT (or first-principles) calculations (Figure 7). Crystallographic data in CIF (Crystallographic Information File) format are useful as the initial structure for geometry optimization. Full optimization is basically preferred (Process A of Figure 7). Full optimization means both the cell parameters and atomic coordinates are optimized. When dispersion correction is not included in DFT or the first-principles, fixed cell parameters or fixed atomic coordinates at the experimental values are frequently used to play the role of the dispersion force (Process B of Figure 7), though this is not true dispersion correction. The reasonable performance of the latest dispersion-corrected DFT has been demonstrated for systems containing weak hydrogen bonds [168]. Potential energy curves for five complexes with weak to medium strong hydrogen bonds have been computed with dispersion-corrected DFT methods. Dispersion-corrected DFT in two totally different approaches (density based and pair-wise corrected) is able to provide a rather accurate and consistent description of weakly hydrogen-bonded systems compared to the high level of ab initio calculations. The typical deviation for intermolecular distances is 0.05 Å and the binding energies are accurate within about 5%-10%. The dispersion-corrected first-principles calculation using the Grimme method [163] and the refinement by Civalleri et al. [169] has been performed intensively by Korter and co-workers over a couple of years to analyze the THz spectra [170][171][172] (Process A of Figure 7). The effects of applying an empirical dispersion correction to solid-state DFT methods were evaluated in the simulation of the crystal structure and low-frequency THz spectrum of naproxen, the non-steroid anti-inflammatory drug. Inclusion of a modified dispersion correction enabled a high-quality simulation of the THz spectrum and the crystal structure of naproxen to be achieved without the need for artificially constraining the unit cell dimensions. The dispersion contributions calculated using existing published parameters were, however, too large for the accurate reproduction of the experimental crystal structure and THz spectrum. In comparison with the empirical dispersion correction by the Grimme approach and the refinement by Civalleri, the use of fixed lattice dimensions obtained from experimental measurements in the DFT calculations (Process B of Figure 7) is an alternative that results in high-quality simulations of the THz spectra. The dispersion-corrected first-principles calculation using the vdW-TS approach has been performed very recently to assign the THz spectrum of weakly hydrogen-bonded systems [72,173]. Full geometry optimization has been performed for both the atomic coordinates and lattice parameters (Process A of Figure 7). Unit-cell parameters are properly converged with a reduction of 1%-2% from experimental values. The results are consistent with the thermal expansion that would occur upon heating from 0 K (first-principles calculations) to a finite temperature (experiments) [174]. The calculated frequencies and the relative intensities excellently reproduce the observed THz spectrum, where no scaling factor is applied to the calculated frequency values. The deviation in frequency is a few percent normally. Some torsion modes and translational modes show relatively large deviations, but the difference from experimental values is approximately 10 cm −1 and is acceptable for the accuracy for the spectral assignment.
The frequency of harmonic vibrations overestimates the experimental fundamental frequencies at the ab initio self-consistent-field (SCF) level of theory [175][176][177][178], due to basis-set incompleteness, electron correlation effects and anharmonicity. Nowadays, the basis set converges nearly to completeness, and the electron correlation is incorporated properly by several methods. Even with the use of a basis set large enough and the adequate incorporation of electron correlations, there remains a deviation attributed to anharmonicity. Several scaling factors were proposed for harmonic frequencies obtained at different levels of ab initio, DFT and semiempirical methods [179]. The deviation of the potential energy well of a simple harmonic oscillator (Figure 8a) from the actual potential energy well of an anharmonic oscillator (Figure 8b) is compensated for by scaling factors smaller than 1.0. The simple harmonic oscillator has a potential energy well that is parabolic in shape, symmetrical about the equilibrium bond length. The potential energy function contains a single quadratic term. When treaded with an anharmonic oscillator, the potential energy well is asymmetric, and dissociation can occur. The potential energy function contains asymmetric cubic, quartic, etc., terms. If one assumes that the form of the potential energy curve is that defined by the Morse potential, this is a good, although not a perfect, approximation to the real potential. As can be seen from the shape of the potential, the deviation from the Morse potential is large at the high-frequency end, and monotonic scaling by a factor is rational. Significant and systematic overestimations of the predicted vibrational frequencies were observed in solid-state calculations of co-crystallized H 2 O molecules with a small biomolecule, L-serine [180]. Evidence provided by the comparison of the experimental and calculated vibrational frequencies for both the protonated and deuterated L-serine· H 2 O solids indicates the presence of significant anharmonicity in the observed lattice vibrations. Vibrational anharmonicity may actually play a much larger role in the interpretation of the THz spectra of hydrates in contrast to their corresponding anhydrous forms. The incorporation of water molecules into the crystal structure of pharmaceuticals must be monitored, since the hydrate form can be easily created or destroyed during the manufacturing process, and the crystalline hydrate often has drastically different physical properties in comparison to the anhydrous form [181].
The scaling factor to compensate for the overestimation of the predicted vibrational frequencies of the harmonic approximation is a few percent less than 1.0 at the latest high-level of quantum chemical calculation. By applying the scaling factor to the calculated harmonic frequencies in the THz region, the value shifts by a few cm −1 . Anharmonic vibrational spectra, including the THz region of sulfuryl halides SO 2 X 2 (X = F, Cl, Br) [182] and hydrogen-bonded systems [183], have been computed by the ab initio or first-principles calculations using the vibrational self-consistent field (VSCF) method. The anharmonic frequencies are lowered in the mid-IR region, compared with harmonic frequencies. In the THz region, the anharmonic frequencies are red-shifted in SO 2 X 2 , but those are rather blue-shifted in the hydrogen-bonded systems. By the finite displacement procedure, the anharmonic infrared spectrum of ClONO 2 and BrONO 2, including the THz region, has been computed [184]. The anharmonic frequency shows a blue shift in the THz region, although there is a red shift in the mid-IR region. Calculated anharmonic vibrational frequencies of a benzoic acid dimer linked by two hydrogen bonds also gives both red and blue shifts in the THz region compared with the harmonic frequencies [119]. Complicated potential wells seem to cause the anharmonicity in the low-frequency region, such as a double-well potential with a very low energy barrier and a very flat energy potential that causes a large amplitude oscillation.

Temperature-Dependent THz Spectrum
The temperature dependence of the THz spectrum allows access to vital information about the structure and vibrational dynamics of solids. Drastic temperature dependences are frequently observed in the THz spectral region, while in some cases, the THz spectra are insensitive to temperature.
Hydrogen-bonded networks are easily broken compared to covalent bonds and facilitate structural changes at room temperature ( Figure 9) [185]. The stretching vibration of weaker kinds of hydrogen bonds can be detected only at low temperatures. Among the spectral changes of vibrational absorptions by temperature, polarization wave, isotope, and so on, the temperature dependence is the most effective for revealing the weak intermolecular interaction and global nuclear motion appearing in the THz spectral region. There is quite a good possibility that the temperature dependence of the THz spectra is related to the formation of weak bonds. The decrease of temperature causes an increase of the number of peaks, a shift of the frequency values, a reduction of the peak width, and so on, in the THz spectra. A number of vibrational modes of hydrogen-bonded purine and adenine were found to shift to higher intensities and frequencies in the FIR vibrational spectra using THz-TDS, as the temperature was reduced [186]. Nucleic acid bases are the building blocks of DNA and ribonucleic acid (RNA), and there have been numerous experimental investigations of nucleic acid bases and nucleosides using mid-IR and Raman spectroscopy, as well as neutron inelastic scattering [187][188][189][190]. Adenine is one of two common purine bases and is found in both DNA and RNA. The intermolecular hydrogen-bonding signature of self-assembled phospholipids was studied by means of THz spectroscopy, and significant information at the molecular level has been provided [191]. Phospholipid layers allow the monitoring of the IR spectrum via well-defined assays and via the temperature-dependent phase transition coupled to specific intra-and intermolecular hydrogen bonding and to the interaction with water. The signals typical for bound water in the mid-IR region showed an analogous temperature dependence for the mixed and pure lipids, whereas the hydrogen-bonding features in the THz region showed a very different behavior for the two types of lipid films. The contribution of water molecules does not exclusively dominate the spectral range in the THz region for this type of film (or at least, it is only one of the possible contributions), and any type of hydrogen bonding may contribute.
Larger molecules sometimes show complicated THz spectra even at room temperature. Furthermore, the peak splitting due to crystal packing forces, the shift of wavenumbers and the sharpening of peaks occur at low temperature. The gas-phase calculation of such a large molecule can be an alternative to the solid-state calculation in assigning the room-temperature spectrum, if solid-state calculations are difficult to perform with reliable accuracy, due to the insufficient capacity of computer power. The measurement of related molecules that are a part of a large molecule is informative for the peak assignment by the gas-phase calculation. Riboflavin consists of a ring part and a side chain. The set of related compounds are alloxazine and lumichrome for the ring part and D-mannitol for the side chain ( Figure 10). A strong peak was observed at 99.5 cm −1 at room temperature in the THz spectrum of riboflavin, and the band has a shape and temperature-dependent behavior similar to the band of D-mannitol at 96.6 cm −1 . The strong peak of riboflavin at 99.5 cm −1 is assigned to the torsion mode of the ribityl group with the gas-phase calculations [139]. Disaccharide, α,α-trehalose dihydrate has two water molecules embedded in the molecule. Two distinct peaks were observed in the room-temperature THz spectrum; those are assigned to the torsion mode of the hydroxymethyl group linked to the hydrogen atoms of two water molecules with the gas-phase calculation [192]. Translational vibrations between molecules are clearly observed at low temperature, and the stretching-vibration mode of the weak hydrogen bond is mixed in the translational vibration. The difference in the network of weak hydrogen bonds is reflected in the low-temperature THz spectrum. Benzoic acid, salicylic acid and acetylsalicylic acid (aspirin) have different networks of weak hydrogen bonds in crystals. Those molecules have a different substituent at the ortho position. All of the three compounds take a centrosymmetric dimer in crystals, but the networks of their weak hydrogen bonds completely differ from each other. Benzoic acid has a layered structure. Salicylic acid takes on a non-layered structure, and aspirin forms a one-dimensional-chain structure sustained by alternating the carboxylic acid and acetyl groups. The stretching vibration of the weak hydrogen bond calculated with the solid-state DFT demonstrates the respective bonding. The revealed stretching vibrations of weak hydrogen bonds between molecules are absent in benzoic acid; there are three in salicylic acid and one in aspirin [72,119,173].

Polarization-Dependent THz Spectrum
The polarization-dependent THz spectrum determines the orientation of the vibrational dipole moments. The THz absorption spectrum of powdered samples reflects the average property of randomly orientated molecules. On the other hand, the THz spectrum of single crystals shows the polarization and/or angle dependence [192][193][194][195]. Another approach to detect the polarization dependence of the THz spectrum is the measurement of polycrystalline thin film, which is planar-ordered relative to the metal plate and, hence, relative to the THz polarization [196][197][198].
Several reports have been given on the polarization dependence at the low-frequency region below 100 cm −1 . The polarization dependence was observed both for the sharp and broad peaks. The FIR spectrum of crystalline α,α-trehalose dihydrate measured at 4 K with a polarized THz wave showed a clear polarization dependence, especially in the sharp peaks below 60 cm −1 [192]. The angle-dependent THz spectra of the L-cysteine and L-histidine single crystals measured with THz-TDS indicated a clear hydrogen-bond peak as a result of the hydrogen-bonded network [193]. The powder spectra measured at low temperature show more intermolecular vibrational modes than those measured at room temperature. The FIR spectra of (100), (010) and (001)-oriented 1,3,5-trinitro-S-triazine single crystals showed a number of discrete absorptions ranging from 10 to 100 cm −1 , and the high dependence of absorption on the orientation of the THz polarization with respect to the crystallographic axis [194]. The polarization-dependent spectra of a stretched polyhydroxybutyrate crystal determined the direction of the vibrational transition moments [195]. The interpretation of the observed polarization dependences is not sufficient without a spectral analysis by solid-state calculations. The spectral change of planar-ordered samples in the measurements with narrow-line waveguide THz-TDS by Grischkowsky and co-workers is well correlated with the direction of the vibrational motion with the precise first-principles calculations ( Figure 11) [192,[196][197][198]. Figure 11. Displacement vectors of the decreasing and unchanged modes with a film of aspirin [192].

Conclusions
THz vibrations are a suitable indicator for the detection of weak bonding in crystals. Though FIR spectroscopy, synonymous with THz spectroscopy, first emerged in 1892, the spectral research in this area did not progress much, due to the broadness and weakness of the peaks in this spectral region. In addition, computational codes to analyze the spectrum were incomplete until very recently. With the recent development of THz technology and an efficient computational code, THz spectroscopy has attained a position as a tool for clarifying the formation of weak bonds in crystals. The THz spectra of several materials have been accumulating day by day. From the experimental point of view, the temperature and polarization dependence provide additional information about the peak characterization. From a theoretical point of view, the incorporation of dispersive interaction and anharmonicity into the solid-state first-principles (or DFT) calculations increases the accuracy of the assignment and avoids uncertainty.
Molecular crystals are assembled with noncovalent packing forces. The packing force is very weak, but the assembled crystal shows an identical structure or, at most, the difference of polymorphs with different network patterns. The stretching vibration of weak hydrogen bonds has long been expected to be in the THz region. Explicit reference to the stretching vibration has just begun, since appropriate dispersive interactions have been incorporated into the first-principles or DFT calculations very recently. Dispersive interaction is not incorporated in pure DFT, but is crucial in weak hydrogen bonds, which are at the vdW end of the types of hydrogen bonds. Due to the development of computational codes, including dispersive interactions, lattice parameters after full optimization deviate by less than 1% from crystallographic data. The effect of dispersive interactions on the frequency calculation of THz vibrations is significant, and thus, the incorporation of dispersive interactions into the calculation of THz vibrations is unavoidable for obtaining a clear assignment. Experimental and theoretical techniques have reached the point of applying THz spectroscopy to the areas of medicine, pharmacy, agronomy, and so on, beyond basic research.