First-Principles Calculation of Terahertz Absorption with Dispersion

Article pubs.acs.org/JPCC

First-Principles Calculation of Terahertz Absorption with Dispersion Correction of 2,2′-Bithiophene as Model Compound Jongtaek Kim,†,‡ O-Pil Kwon,§ Mojca Jazbinsek,∥ Young Choon Park,† and Yoon Sup Lee*,† †

Department of Chemistry, Korea Advanced Institute of Science and Technology (KAIST), Daejeon 305-701, Korea Department of Basic Science, Korea Air Force Academy (KAFA), Cheongju 363-849, Korea § Department of Molecular Science and Technology, Ajou University, Suwon 443-749, Korea ∥ Institute of Computational Physics, Zurich University of Applied Sciences (ZHAW), 8401 Winterthur, Switzerland ‡

ABSTRACT: Terahertz absorption of organic materials is closely linked to molecular arrangements and their intermolecular interactions and is important for material identification as well as THz generation. Theoretical calculations of solid-state vibrations known as phonons help to understand intermolecular interactions responsible for THz absorption but frequently are of limited use without considering dispersion interaction. In this study, we have calculated the THz phonon modes of an organic model crystal 2,2′-bithiophene, considering dispersion intermolecular interactions assuming the fixed cell dimensions. Both energies and intensities of phonon modes at low frequencies were interpreted concentrating on the intermolecular level in conjunction with hydrogen bonds and showed an excellent agreement with the experimental results. This approach to identify the phonon modes responsible for strong THz absorptions and to interpret those modes in terms of intermolecular vibrations is also expected to be applicable to the field of THz generation using nonlinear optical organic crystals.

■

THz.17 However, since DAST exhibits an optical phonon resonance at about 1.1 THz, THz waves generated by DAST crystals show a deep gap at this frequency even at the most optimal pump wavelength.17 Therefore, THz absorption properties in organic materials are important not only for material characterization but also for THz-wave generation. THz absorption properties of organic materials can be examined by theoretical calculations by considering the corresponding solid vibrations (i.e., phonons). Since the calculation methods were recently implemented in firstprinciples programs with the advent of modern polarization and linear response theories, theoretical studies on THz absorption became possible. The phonon modes of molecular crystals in the THz frequency region strongly depend on intermolecular interactions, whereas for the conventional infrared frequency region intramolecular interactions are relevant.18,19 Therefore, the use of crystals instead of molecules20−22 or clusters23 is required to accurately reproduce THz absorption of molecular crystals, and the dispersion interactions corresponding to intermolecular interactions are very important. In many previous reports involving theoretical calculations of THz absorption properties, the dispersion interactions are

INTRODUCTION Organic materials show specific absorption characteristics in the terahertz frequency region (0.1−10 THz), which are important for spectroscopic research and applications such as characterization and identification of explosives,1−5 drugs,6−10 and biological materials.11−13 THz spectroscopy and imaging are also important for identification of defects and inclusions in organic materials. Nonlinear optical organic materials are additionally very important materials for THz wave generation using infrared and near-infrared pump laser sources.14,15 Their figures of merit for THz generation and detection are more than 1 order of magnitude higher compared to their inorganic alternatives,16 which is due to their high nonlinear optical susceptibilities combined with a low dielectric dispersion from the low-frequency to the optical frequency range. Also for these materials the absorption properties in the THz range play a very important role in efficient THz generation. First of all, efficient THz generation is possible when phase matching between pump optical and generated THz waves in the material can be achieved16 and therefore depends on the refractive indices in the THz range, which are inherently coupled to absorption. Moreover, the emitted THz wave can be absorbed by the source material itself, which will limit the efficiency even at perfect phase-matching configurations. For example, a benchmark THz source material, DAST (N,N-dimethylaminoN′-methylstilbazolium 4-methylbenzenesulfonate) crystal, can be phase matched to generate THz frequencies of about 1.1 © 2015 American Chemical Society

Received: March 19, 2015 Revised: May 15, 2015 Published: May 15, 2015 12598

DOI: 10.1021/acs.jpcc.5b02661 J. Phys. Chem. C 2015, 119, 12598−12607

Article

The Journal of Physical Chemistry C

Figure 1. Arrangements of molecules and their interactions in BT crystal: left (side view), right (top view).

The dispersion interactions originating from fluctuating charge distributions between molecules are of great importance for phonon calculations in the THz region. The dispersion forces are nonlocal electron−electron correlations omitted in the DFT with local or semilocal functionals. Out of various schemes to solve this problem, two methods including van der Waals (vdW) interactions in DFT functionals were used in this study: the DFT-D2 approach by Grimme35 and the vdW-DF approach by Dion et al.36 In the DFT-D2 method, the vdW interactions optimized for several popular DFT functionals are simply added to DFT functionals via a pairwise force field, so that it does not affect the electron density of the parent DFT functionals (eq 1). In the vdW-DF method, on the other hand, the vdW interactions are directly included to the parent DFT functional as a nonlocal correlation part, the correction of PBE correlation is replaced by the local-density approximations (LDA) correlation, and the GGA exchange functional is optimized for the vdW correction (eq 2).37,38

neglected or only partially considered, and the agreement with experimental results is limited.24−27 For example, the calculation of THz phonon modes of the oligothiophene crystal 2,2′-bithiophene, which is hereafter noted as “BT”, performed by Hermet et al.27 showed a relatively good agreement with the experimental data. However, the calculated THz phonon modes appear somewhat complicated and are still not exactly agreeing with the experimental observations. The lack of dispersion interactions in the phonon calculations in the work of Hermet et al. may be one of the reasons for the discrepancy between the calculated and the measured THz spectra.27 Herein, we calculate the THz phonon modes of an organic model compound considering dispersion interactions and investigate the effect of these interactions on the phonon modes. We used the centrosymmetric BT crystal having a relatively simple crystal structure with a small number of atoms (only two molecules) per unit cell (Z = 2) as model compound, but the same calculation principles can be straightforwardly applied to noncentrosymmetric (nonlinear optical) crystals, such as DAST, which may be used for THz-wave generation. Therefore, the THz phonon calculation methods of this study are useful not only for organic materials characterized by THz spectroscopy and imaging but also for organic nonlinear optical materials for THz generation.

E DFT‐D = E KS‐DFT + Edisp

(1)

E XC = ExGGA + EcLDA + Ecnl

(2)

Calculations of dynamical matrix and normal modes were performed in the harmonic approximation using the direct method39 implemented in the PHONON 4.3 software.40 In the infrared absorption of long-wavelength light, only some modes at the Γ-point (k = 0) are infrared active. For Hellmann− Feynman force calculations in VASP code, the 48 atomic displacements of the asymmetric unit with 8 atoms of the BT crystal were generated by the PHONON program with δ = 0.05 Å displacements in each Cartesian direction. An empirical 4 cm−1 fwhm Lorentzian line shape has been applied to simulated THz spectra for comparison with the experimental spectrum.

■

COMPUTATIONAL DETAILS Periodic density functional theory (DFT) calculations of atom relaxations, Hellmann−Feynman (HF) forces, and Born effective charge tensors (BECTs) were performed in the framework of fixed unit cell dimensions using the Vienna ab initio simulation package (VASP)28−30 with a plane wave energy cutoff of 400 eV. The Perdew−Burke−Ernzerhof (PBE) exchange-correlation functional31 was used for the gradient approximation (GGA) to the exchange-correlation (XC) energy. The interactions between ions and electrons were described with the projector augmented wave (PAW) method32 in the real space representation. A fine grid 6 × 6 × 6 Monkhorst−Pack33 k-points mesh was used over the Brillouin zone for the accuracy of calculations. The BECTs were obtained by the linear response method.34

■

RESULTS AND DISCUSSION Intra- and Intermolecular Structures. For the BT crystal, X-ray diffraction parameters were determined at low temper5 ature 133 K: monoclinic space group P21/c (C2h ), two molecules per unit cell (Z = 2), and lattice parameters with a = 7.734 Å, b = 5.729 Å, c = 8.933 Å, and β = 106.72°.41 Figure 1 12599


Article

The Journal of Physical Chemistry C Table 1. Intramolecular Structural Parameters of BT exptla S1−C2 C2C3 C3−C4 C4C5 C5−S1 C2−C2′

PBEb

1.713 1.432 1.444 1.357 1.698 1.448 RMSD

PBE-D2b

Bond Lengths (Å) 1.737 1.745 1.387 1.387 1.419 1.423 1.377 1.377 1.716 1.725 1.445 1.446 0.001 0.001 Angles (deg) 110.2 110.3 113.1 113.2 112.6 112.7 111.6 111.6 92.5 92.2 120.9 120.7 128.9 128.9 6.2 6.2 Dihedral Angles (deg) 180.0 180.0 180.0 180.0 0.0 0.0

1.734 1.385 1.417 1.376 1.714 1.443 0.001

S1−C2C3 C2C3−C4 C3−C4C5 C4C5−S1 C5−S1−C2 S1−C2−C2′ C3C2−C2′

112.4 108.0 114.9 112.1 92.4 121.2 126.4 RMSD

110.2 113.1 112.6 111.6 92.5 120.8 129.0 6.2

S1−C2−C2′-S1′ C3C2−C2′C3′

180.0 180.0 RMDS

180.0 180.0 0.0

PBE-DFb

PBE-D2*c 1.734 1.386 1.419 1.378 1.715 1.444 0.001

PBEa 1.734 1.390 1.420 1.380 1.715 1.444 0.001

110.3 113.0 112.6 111.6 92.4 120.9 128.8 5.9

110.3 113.0 112.5 111.6 92.6 120.8 128.9 6.0

180.0 180.0 0.0

180.0 180.0 0.0

a c

From ref 27. bIon relaxations with the PBE functional only, with semiempirical or ab initio dispersion correction to the PBE functional, respectively. Full relaxation of ions and cell with semiempirical dispersion correction.

Table 2. Intermolecular Structural and Cell Parameters of BT hydrogen bond lengths (Å) exptl HB1 (C3···C2′) HB2 (C4···C2) HB3 (C5···S1′) HB1′ (C3−H···C2′) HB2′ (C4−H···C2) HB3′ (C5−H···S1′)

exptl (fixed cell) PBE-D2* Δ %

a

3.4628 3.6114 3.4882 RMSD 2.8894 3.0029 3.0455 RMSD

PBE

b

3.5168 3.6155 3.4583 0.0013 2.8522 2.9234 3.0141 0.0029

PBE-D2b 3.4659 3.5809 3.4637 0.0005 2.7749 2.8485 3.0163 0.0126 cell parameters

PBE-DFb

PBE-D2*c

3.5219 3.6234 3.4635 0.0014 2.8669 2.9384 3.0353 0.0016

3.3924 3.5358 3.3506 0.0099 2.7491 2.7681 2.9732 0.0267

a (Å)

b (Å)

c (Å)

β (deg)

V (Å3)

7.734 7.356 −0.378 −4.9%

5.729 5.619 −0.110 −1.9%

8.933 8.995 0.062 0.7%

106.72 105.20 −1.52 −1.4%

395.804 371.773 −24.031 −6.1%

a

Experimental structure from ref 27. bIon relaxations with the PBE functional only, with semiempirical or ab initio dispersion correction to the PBE functional, respectively. cFull relaxation of ions and cell with semiempirical dispersion correction.

displays the atom labeling scheme and the packing of BT molecules. The molecular conformation of BT is planar in crystalline phase, which is quite different in gas phase having a nonplanar conformation with trans-configuration of the sulfur atoms. The intermolecular interactions of BT crystal, which are greatly important with regard to the THz absorption, are characterized by three symmetrically unique hydrogen bonds. Figure 1 shows the three weak hydrogen bonds: C3−H···π-C2′ (i.e., HB1), C4−H···π-C2′ (i.e., HB2), and C5−H···S1′ (i.e., HB3). In addition to HB1 and HB3 bonds previously pointed out by Eijck et al.,42 HB2 bond almost equivalent to HB1 bond in terms of bonding type and the distance is considered here to evaluate intermolecular interactions as well as to investigate intermolecular vibration motions. Therefore, BT molecules in

crystalline state primarily form dimers through the weak C− H···π hydrogen interactions (i.e., HB1 and HB2) like in durene or naphthalene crystals. The dimers secondarily interact with one another via the weak C−H···S hydrogen interactions (i.e., HB3), which leads to the final crystal structure with herringbone-type interactions. This secondary interaction helps the twisted BT molecule to be planar with the S−C− C−S dihedral angle of 0°. Changes of these intermolecular (i.e.HB1 and HB2) and interdimeric (i.e., HB3) interactions by temperature variation may lead to changes of molecular arrangements and planarity. Interestingly, no significant π−π stacking interactions were found, probably due to the large intermolecular distance along the π−π stacking direction. The experimental and the calculated intramolecular structural parameters of BT crystal are listed in Table 1. The results are 12600


Article


Figure 2. Experimental (9 K) and calculated THz spectra of BT crystal: PBE (ion relaxation), PBE-D2 (ion relaxation with semiempirical dispersion correction), PBE-DF (ion relaxation with ab initio dispersion correction). The experimental data (a−d) and PBE calculation result (a) obtained from ref 27 were reproduced under the permission of the author and the copyright of the publisher. The experimental data obtained at 9 K was used for comparison with calculations due to peak distinction.

correction to the DFT functional. The PBE-D2 and PBE-DF simulations are for ionic relaxation with semiempirical or ab initio dispersion correction to the DFT functional, respectively. Besides, relaxation for both ions and cell with semiempirical dispersion correction (D2) was performed with the PBE-D2* simulation, since the ab initio dispersion correction (DF) is not available for cell relaxation at present. Because X-ray diffraction does not measure the hydrogen atom positions to great accuracy, the heavy atoms were used for the evaluation of hydrogen bond lengths and RMSD. Over all other methods, the PBE-D2 provides the best performance for evaluation of hydrogen bond lengths with the smallest RMSD value of 0.0005 Å.43 The PBE and PBE-DF methods show larger similar RMSD values of 0.0013 and 0.0014 Å, respectively. The PBE-D2 method excellently reproduces the two HB1 and HB3 bonds, and later the experimental THz absorption spectrum. The PBE-D2* simulation seems to overestimate dispersion interactions,26 which causes large volume change (−6.1%) and asymmetric cell contraction and expansion resulting in large deviation of hydrogen bond lengths from the experimental ones as noted in the largest RMSD value of 0.0267 Å. The a-direction with the

similar for various DFT methods used in this study and to those of the previous work by Hermet et al.27 Bond lengths and dihedral angles are well reproduced by calculations with rootmean-square deviations (RMSDs) of 0.001 Å and 0.0°, respectively. However, the RMSDs of bond angles having similar large values around 6.0° are mostly caused by the C2 C3−C4 angle contraction and subsequently the other bondangle variations in the experimental crystal structure. Although the excellent reproduction of intramolecular parameters by theory is prerequisite for the calculation of intermolecular interactions and THz absorption spectra, to the best of our knowledge, there are no solid-state DFT methods developed to properly describe these large bond-angle variations caused by the crystal-field effect. As will be shown in the later section, the intermolecular interactions are more important and dominant factors for the THz absorption calculation than the above bond-angle variation. The calculated intermolecular structural parameters concerned with the three main hydrogen bonds are listed in Table 2 and compared with the experimental parameters. The relaxations were performed in the regime of fixed unit cell. The PBE simulation is for ionic relaxation without dispersion 12601


Article


different characteristics. The previous work relatively well describes the low-frequency region below 120 cm−1 with a small red-shift but does not properly describe the splitting of the three modes f, g, and h around 140 cm−1, which are merged at around 135 cm−1 as seen in Figure 2a and Table 3. Our PBE simulation as seen in Figure 2b shows little more red-shift in the region below 120 cm−1, which mainly contributes to larger RMSD value of 10.5 cm−1. The splitting of the three modes f, g, and h is relatively well described by the present PBE simulation (see Table 3). However, it is known that these kinds of simulations absent of dispersion interaction should not accurately reproduce the experimental THz spectra strongly depending on intermolecular interactions. Therefore, the dispersion interaction is important for accurate THz absorption simulation of the BT crystal and is introduced in this work at the semiempirical or ab initio level with the plane waves (PW) based solid-state density functional theory. The PBE-D2 simulation with semiempirical dispersion correction exhibits quite different spectral features with the smaller RMSD values of 6.4 cm−1, compared with the PBE simulations. As seen in Figure 2c, the PBE-D2 method accurately reproduces the two main peaks c (78 cm−1) and g (143 cm−1). The red-shift in the region below 120 cm−1, which is observed for the PBE simulation, is here notably reduced, and the splitting of the three modes f, g, and h around 140 cm−1 is well reproduced. In addition, the intensities of the modes f and i around the mode g are improved compared with both PBE methods, which leads to the spectral shape more similar to that of the measurement. The intensities are on the order of e ≈ i < f for both experiment and PBE-D2 and i ≪ f < e for PBE. This spectral agreement in terms of both vibration energy and intensity implies that the intermolecular interaction enhanced by dispersion force in the regime of fixed unit cell more accurately describes potential energy surfaces (PESs) responsible for lattice vibration motions. The increased intermolecular interactions induce here an upshift of frequencies and variation of intensities, leading to a more accurate reproduction of the experimental THz spectrum. It should be noted that the calculations of this work were done at 0 K with the unit cell parameters obtained at 133 K, while the experimental THz absorption spectrum measured at 9 K was used for comparison to utilize the discernible peak intensities among vibration modes. Thus, the measured spectrum shows a blue-shift because of cell contraction according to temperature decrease.27 In addition, the BT crystal with much weaker intermolecular interactions than the quaterthiophene crystal tends to be disordered at room temperature.27,44 These factors may lead to a limited reproduction of the experimental THz absorption. The calculation method PBE-D2 considering dispersion interaction at the fixed cell serves as an approximate solution for the spectral discrepancy, providing a good agreement with the experimental data. However, the calculation still does not reproduce the mode h, which exists as the second largest intensity peak in the experimental spectrum at 9 K but disappears as temperature increases.27 This peak behavior is due to the experimental causes such as thermal contraction and thermal disorder. A discussion of the mode h with much smaller intensity in calculation than in measurement was tried in detail at the end of the Vibrational Mode Analysis section. The PBE-DF method with ab initio dispersion correction gives a quite different spectral shape despite the small RMSD value of 6.9 cm−1 comparable to that of the PBE-D2 method.

three hydrogen bonds experiences the largest contraction (−4.9%) due to the overestimation of dispersion forces, the bdirection with longer C−H···S interactions undergoes smaller contraction (−1.9%) due to slipped π−π stacking, and the cdirection with S···S and π···π interactions experiences small expansion (0.7%). Although all the RMSD magnitudes for hydrogen bond lengths are very small with similar intramolecular parameters, the differences could be significant for the THz absorption spectra which are sensitive to the delicate changes of intermolecular interactions as discussed in a later section. Taking into account semiempirical dispersion correction and heavy atom criterion into evaluating hydrogen bonds is of great importance for the exact description of intermolecular interactions and the THz absorption. It seems that the ab initio dispersion correction is not yet optimal for molecular crystals such as the BT crystal. Therefore, the PBE-D2 method using semiempirical dispersion correction in the fixed cell dimensions is recommended for THz absorption calculations, and this is demonstrated in the following THz spectra simulations using the aforementioned various DFT methods. THz Spectra Simulation. Figure 2 shows the four simulated THz spectra of BT crystal plotted together with the experimentally measured one at 9 K. The experimental THz spectrum exhibits well-resolved nine spectral features. The calculated and measured THz absorption frequencies of the BT crystal are given in Table 3 together with RMSD values for Table 3. Experimental and Calculated IR-Active THz Frequencies (cm−1) of BT Crystal mode

exptla

PBEa

PBEb

PBE-D2b

PBE-DFb

PBE-D2*c

a b c d e f g h i

47 62 78 105 119 133 143 149 161 RMSD

39 61 75 95 118 133 135 138 151 7.1

30 55 72 85 115 132 137 144 150 10.5

35 60 78 93 119 139 146 148 156 6.4

35 71 72 98 119 141 150 152 159 6.9

52 71 88 118 139 147 151 160 169 11.6

a Experimental structure. bIon relaxations with the PBE functional only, with semiempirical or ab initio dispersion correction to the PBE functional, respectively. cFull relaxation of ions and cell with semiempirical dispersion correction.

exact comparison. In view of relative intensities, all the peaks except the mode h (149 cm−1, see Table 3) are well reproduced by the four simulations with relatively different peak positions. The second strongest mode h (149 cm−1) of the experimental spectrum is calculated with similar frequencies but with negligible intensities leading to disappearance of doublet feature formed by the modes g and h. Therefore, the first and the third strongest modes g and c are used as the two main peaks to evaluate the simulation quality of the spectral shape. The results of four simulations show some differences. The PBE simulations from ref 27 (Figure 2a) and this work (Figure 2b), which do not consider dispersion interactions within fixed cell dimensions, show a moderately good agreement with the experimental values with RMSD values of 7.1 and 10.5 cm−1, respectively. The PBE calculation of this work produces the spectral shape roughly similar to that of ref 27, with somewhat 12602


Article


Figure 3. Representations of THz phonon modes of BT crystal along the b-axis. The modified Cartesian coordinates x′ and z′ are for each translational motion. The intermolecular motions are additionally shown with circles (butterfly), bold straight arrows (optical translation), and curved arrows (torsion).

unit cell dimensions induce large variations in the THz spectral region. Therefore, there is a need to adjust the semiempirical dispersion parameters with regard to cell relaxation processes to reproduce the cell volume as good as the measured one, as pointed out by Korter et al.45−50 This approach is not suitable for the present work in the absence of X-ray diffraction structure obtained at 9 K. Vibrational Mode Analysis. The vibrational motions were empirically analyzed with the eigenvectors displayed in two directions, parallel or perpendicular to the molecular bonding axis in Figures 3 and 4, respectively. Some graphical objects were added to assist the empirical analysis of intermolecular vibrations. Since the optical translation motions of the BT crystal do not correspond to the Cartesian axes of the crystal coordinate system, the modified Cartesian axes (x′, y, z′) were defined along the translation motions. The assignments of vibrational motions are listed in Table 4 in terms of both intramolecular and intermolecular components.

As seen in Figure 2d, the PBE-DF simulation is not in good agreement with the experimental spectrum especially for the two main peaks, with the mode c red-shifted (78 cm−1 → 72 cm−1) and the mode g blue-shifted (143 cm−1 → 150 cm−1). In addition, the blue-shifted g mode convolved with the weak h (152 cm−1) mode does not correspond to the second strongest h mode of the experimental THz spectrum. It seems that ab initio corrections of dispersion interactions to DFT functionals are not so properly optimized as to accurately describe intermolecular interactions in molecular crystals such as the BT crystal. For PBE-D2* simulation with cell relaxation as well as semiempirical dispersion correction, most of peaks are significantly blue-shifted with the largest RMSD value of 11.6 cm−1 as shown in Table 3, whose results hence were not plotted. As discussed in the intermolecular section above, the dispersion interactions seem to be overestimated during cell relaxation process resulting in somewhat large variations of unit cell parameters and volume. These asymmetric variations of 12603


Article


Figure 4. Representations of THz phonon modes of BT crystal along the ac-plane. The modified Cartesian coordinates x′ and z′ are for each translational motion. The intermolecular motions for bending are additionally shown with the curved arrows. The dashed lines are for the existing hydrogen bonds (i.e., HB1, HB2, and HB3), and the dotted lines are for newly forming hydrogen bonds (i.e., C−H···π-C and C−H···S) with regard to intermolecular motions. The end arrows of dashed and dotted lines stand for repulsive motions leading to breaking of hydrogen bonds. The half end arrows mean that only one end goes apart from the other.

Table 4. Experimental and PBE-D2 Calculated IR-Active Frequencies (cm−1) and Mode Descriptions of BT Crystal in Terms of Intramolecular and Intermolecular Motionsa mode

exptlb

PBE-D2c

symmetry

intramolecular motion

intermolecular motion

a b c d e f g h i

47 62 78 105 119 133 143 149 161

35 60 78 93 119 139 146 148 156

Au(I) Au(I) Bu(I) Bu(I) Au(I) Au(I) Bu(I) Bu(I) Au(I)

str. translation in ac (z′) str. torsion mdr. translation along b (y) mdr torsion str. torsion str. translation in ac (x′) mdr. torsion str. bending wk. butterfly str. buttefly str. bending mdr. bending mdr. butterfly mdr. butterfly mdr. bending

OOP translation in ac (z′) IP torsion along z′ OOP translation along b (y) OOP torsion along z′ OOP torsion along z′ OOP translation in ac (x′) IP torsion along z′ IP bending along b (y) IP butterfly along b (y) IP butterfly along x′ OOP bending along b (y) OOP bending along b (y) OOP butterfly along b (y) OOP butterfly along x′ IP bending along b (y)

a

Note: str. (strong), mdr. (moderate), wk. (weak), OOP (out of phase), IP (in phase). bExperimental structure. cIon relaxations with the PBE functional with semiempirical dispersion correction to the PBE functional.

The low-frequency region is populated by nonstretching vibrations (i.e., torsion, bending, and butterfly) more than stretching vibrations, since the former with bond angle changes has shallower potential energy surfaces than the latter with bond length changes. Usually, torsional motions are likely to

occur at low frequencies, while the bending and butterfly motions come to the upper frequency region with similar energies. Considering these molecular vibrations of BT in the crystal environment (Z = 2), each molecular vibration can be duplicated into two intermolecular dimeric vibrations depend12604


Article


those of torsional vibration modes b and d forming the third group. The mode h with the same intramolecular motions as the mode i, but with the different phase, shows negligible absorption in calculation, whereas it appears as the second largest peak in the measurement. The mode h is not observed at 290 K but shows up at 9 K with high intensity and high frequency making the largest main peak doublet.27 However, this mode h with out-of-phase intermolecular motions with no significant changes in the existing hydrogen bonds has negligible intensity in all four simulations. The reason may be due to the use of structure determination at 133 K for the calculations while the THz absorption measurement was at 9 K. The present simulations employing unit cell parameter determined at 133 K do not properly describe the cell contraction when cooling from 133 to 9 K. The cell contraction is likely to occur in asymmetric manner as found in the results of PBE-D2* simulation performed with cell and ion relaxations under semiempirical dispersion correction. Besides, it is suspected that the BT crystal thermodynamically unstable at room temperature is easily disordered27 and likely frozen into orientationally and/or spatially disordered states at 9 K. It has been shown in literature that even subtle changes in hydrogen bond52,53 and/or conformation54 induce orientational and/or spatial molecular arrangement changes,55,56 leading to significant changes of vibration energies and intensities. These disorders may change C−H···π and C−H···S hydrogen bonds,57 vibration character of mode h, and the dipole moment, which result in the enhanced absorption in experiment.

ing on the interacting phases (i.e., out of phase and in phase) with different energies from the counterpart and the parent monomeric vibration. This is observed for the torsional vibration of BT corresponding to the modes b and d as seen in Figures 3 and 4 and described in Table 4. The butterfly and bending intramolecular motions with similar energies give the four intermolecular dimeric vibrations of the modes f, g, h, and i as seen in Figures 3 and 4 and described in Table 4. The torsional motion was excluded in the description of the modes f ∼ i, although observed, because it has a much lower energy than the bending and butterfly motions. This leads to weak coupling with negligible contribution. In addition to these intermolecular vibrations based on intramolecular vibrations, the lattice vibrations of rigid molecules such as translations are present in the low-frequency region of the BT crystal. This is because potential energy surfaces of weak intermolecular interactions with van der Waals forces are usually shallower than those of intramolecular interactions with strong chemical bonding forces. Exceptionally, the intramolecular torsional motion, which easily leads to the twisted structure as in gas phase without crystal field, has such a shallow potential energy surface that the two intermolecular torsional motions of b and d with different phases are found near the translational motions. Thus, the intramolecular torsional motion is easy to couple with translational motions. Since the torsional motions are perpendicular to the molecular bonding axis z′, the coupling is favored for the translational modes c and e along the x′ and y axes respectively rather than the a mode along the z′ axis, as seen in Figures 3 and 4 and described in Table 4. The mode c is the mixing of translation along y and in-phase intermolecular torsional motion (mode d), while the mode e is the mixing of translation along x′ and outof-phase intermolecular torsional motion (mode b). The energies of translational vibrations also provide information on hydrogen bonding. As the translation energies increase in the order of a (along z′) ≪ c (along y) < d (along x′) modes, BT molecules favor intermolecular interactions along x′ and y axes, which correspond to the directions of HB1, HB2, and HB3 as aforementioned in the crystal structure analysis and seen in Figure 1. Therefore, the motions and energies of crystal vibrations in the low-frequency region should be understood considering both intramolecular and intermolecular motions as supported by this work. The intensities of infrared absorptions of the BT crystal depend on the dipole moment changes by intermolecular vibrations. The lattice vibrations of BT with no dipole moment at the molecular level have usually no contributions to the dipole moment change of the whole crystal but have valid contributions when mixed with symmetry-breaking intramolecular vibrations. The lattice vibration along z′ (mode a) parallel to the molecule’s long axis has almost pure translational motion and hence the smallest intensity,51 whereas those along x′ and y (modes c and e) perpendicular to the molecule’s long axis have effective intensities due to mixing into intramolecular vibrations, with the mode c being the largest. On the other hand, the molecular vibrations of BT have dipole moment changes in the order of bending < torsion < butterfly. Hence, the intermolecular vibration with strong butterfly and bending motions (mode g) leads to the most intense absorption. These two modes c and g make the main two absorption peaks as the first group. The mode f with strong bending and moderate butterfly and the mode i with moderate butterfly and bending make the second absorption group, which is more intense than

■

CONCLUSION We have performed THz absorption calculations considering dispersion interactions in the fixed cell regime. The PBE-D2 simulation considering semiempirical dispersion interactions leads to an excellent description of intermolecular interactions as well as intramolecular structures, and hence the THz absorption spectrum of BT crystal measured at 9 K. The relevant phonon modes have been described not only at the intramolecular level but also at the intermolecular level via intermolecular motion. Our calculation with semiempirical dispersion interaction within fixed cell obtained at 133 K is in a good agreement with the experimental THz spectrum obtained at 9 K. The intensities of the THz spectrum were analyzed in terms of dipole moment change concerning hydrogen bond breaking and making caused by intermolecular vibrational motions. The two main absorptions result from symmetrybreaking phonon vibrations such as butterfly, bending, and torsion motions along the directions (x′, y) of hydrogen bonds perpendicular to the molecular long axis (z′), which leads to effective dipole moment changes and large absorption intensities.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Tel.: +82-42-350-2821. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work is based on the author’s doctoral dissertation58 and was supported by the National Research Foundation of Korea (NRF) grant funded by the Korea government (No. 201112605


Article


(19) Balbuena, P. B.; Blocker, W.; Dudek, R. M.; Cabrales-Navarro, F. A.; Hirunsit, P. Vibrational Spectra of Anhydrous and Monohydrated Caffeine and Theophylline Molecules and Crystals. J. Phys. Chem. A 2008, 112, 10210−10219. (20) Hakey, P. M.; Allis, D. G.; Hudson, M. R.; Ouellette, W.; Korter, T. M. Terahertz Spectroscopic Investigation of S-(+)-Ketamine Hydrochloride and Vibrational Assignment by Density Functional Theory. J. Phys. Chem. A 2010, 114, 4364−4374. (21) Zhang, Y.; Peng, X. H.; Chen, Y.; Chen, J.; Curioni, A.; Andreoni, W.; Nayak, S. K.; Zhang, X. C. A First Principle Study of Terahertz (THz) Spectra of Acephate. Chem. Phys. Lett. 2008, 452, 59−66. (22) Saito, S.; Inerbaev, T. M.; Mizuseki, H.; Igarashi, N.; Note, R.; Kawazoe, Y. First Principles Calculation of Terahertz Vibrational Modes of a Disaccharide Monohydrate Crystal of Lactose. Jpn. J. Appl. Phys., Part 2 2006, 45, L1156−L1158. (23) Yan, H.; Fan, W. H.; Zheng, Z. P. Investigation on Terahertz Vibrational Modes of Crystalline Benzoic Acid. Opt. Commun. 2012, 285, 1593−1598. (24) King, M. D.; Buchanan, W. D.; Korter, T. M. Understanding the Terahertz Spectra of Crystalline Pharmaceuticals: Terahertz Spectroscopy and Solid-State Density Functional Theory Study of (S)(+)-Ibuprofen and (RS)-Ibuprofen. J. Pharm. Sci. 2011, 100, 1116− 1129. (25) Motley, T. L.; Allis, D. G.; Korter, T. M. Investigation of Crystalline 2-Pyridone Using Terahertz Spectroscopy and Solid-State Density Functional Theory. Chem. Phys. Lett. 2009, 478, 166−171. (26) Ortmann, F.; Hannewald, K.; Bechstedt, F. Ab Initio Studies of Structural, Vibrational, and Electronic Properties of Durene Crystals and Molecules. Phys. Rev. B 2007, 75, 195219. (27) Hermet, P.; Bantignies, J. L.; Rahmani, A.; Sauvajol, J. L.; Johnson, M. R.; Serein, F. Far- and Mid-Infrared of Crystalline 2,2′Bithiophene: Ab Initio Analysis and Comparison with Infrared Response. J. Phys. Chem. A 2005, 109, 1684−1691. (28) Kresse, G.; Furthmüller, J. Efficiency of ab Initio Total Energy Calculations for Metals and Semiconductors Using a Plane-Wave Basis Set. Comput. Mater. Sci. 1996, 6, 15−50. (29) Kresse, G.; Furthmüller, J. Efficient Iterative Schemes for Ab Initio Total-Energy Calculations Using a Plane-Wave Basis Set. Phys. Rev. B 1996, 54, 11169−11186. (30) Kresse, G.; Hafner, J. Ab Initio Molecular Dynamics for Liquid Metals. Phys. Rev. B 1993, 47, 558−561. (31) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (32) Kresse, G.; Joubert, D. From Ultrasoft Pseudopotentials to the Projector Augmented-Wave Method. Phys. Rev. B 1999, 59, 1758− 1775. (33) Monkhorst, H. J.; Pack, J. D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188−5192. (34) Giannozzi, P.; Gironcoli, S.; Pavone, P.; Baroni, S. Ab Initio Calculation of Phonon Dispersions in Semiconductors. Phys. Rev. B 1991, 43, 7231−7242. (35) Grimme, S. Semiempirical GGA-Type Density Functional Constructed with a Long-Range Dispersion Correction. J. Comput. Chem. 2006, 27, 1787−1799. (36) Dion, M.; Rydberg, H.; Schroder, E.; Langreth, D. C.; Lundqvist, B. I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. (37) Lee, K.; Murray, E. D.; Kong, L. Z.; Lundqvist, B. I.; Langreth, D. C. Higher-Accuracy Van der Waals Density Functional. Phys. Rev. B 2010, 82, 081101. (38) Klimes, J.; Bowler, D. R.; Michaelides, A. Chemical Accuracy for the Van der Waals Density Functional. J. Phys.: Condens. Matter 2010, 22, 022201. (39) Parlinski, K. Direct Method for Phonon Calculation. Am. Inst. Phys. Conf. Proc. 1999, 479, 121. (40) Parlinski, K. Software PHONON, Cracow, Poland, 2010 . (41) Pelletier, M.; Brisse, F. Bithiophene at 133 K. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1994, 50, 1942−1945.

0027510, 2007-0056095, 2014R1A5A1009799, and 20090093826). Computational resources were in part provided by KISTI (KSC-2011-C2-18).

■

REFERENCES

(1) Witko, E. M.; Korter, T. M. Terahertz Spectroscopy of the Explosive Taggant 2,3-Dimethyl-2,3-Dinitrobutane. J. Phys. Chem. A 2012, 116, 6879−6884. (2) Burnett, A. D.; Kendrick, J.; Cunningham, J. E.; Hargreaves, M. D.; Munshi, T.; Edwards, H. G. M.; Linfield, E. H.; Davies, A. G. Calculation and Measurement of Terahertz Active Normal Modes in Crystalline PETN. ChemPhysChem 2010, 11, 368−378. (3) Hooper, J.; Mitchell, E.; Konek, C.; Wilkinson, J. Terahertz Optical Properties of the High Explosive Beta-HMX. Chem. Phys. Lett. 2009, 467, 309−312. (4) Allis, D. G.; Zeitler, J. A.; Taday, P. F.; Korter, T. M. Theoretical Analysis of the Solid-State Terahertz Spectrum of the High Explosive RDX. Chem. Phys. Lett. 2008, 463, 84−89. (5) Allis, D. G.; Korter, T. M. Theoretical Analysis of the Terahertz Spectrum of the High Explosive PETN. ChemPhysChem 2006, 7, 2398−2408. (6) Zheng, Z. P.; Fan, W. H.; Li, H.; Tang, J. Terahertz Spectral Investigation of Anhydrous and Monohydrated Glucose Using Terahertz Spectroscopy and Solid-State Theory. J. Mol. Spectrosc. 2014, 296, 9−13. (7) Pellizzeri, S.; Delaney, S. P.; Korter, T. M.; Zubieta, J. Using Solid-State Density Functional Theory and Terahertz Spectroscopy To Spectroscopically Distinguish the Various Hydrohalide Salts of 5-(4Pyridyl)tetrazole. J. Mol. Struct. 2013, 1050, 27−34. (8) King, M. D.; Buchanan, W. D.; Korter, T. M. Identification and Quantification of Polymorphism in the Pharmaceutical Compound Diclofenac Acid by Terahertz Spectroscopy and Solid-State Density Functional Theory. Anal. Chem. 2011, 83, 3786−3792. (9) Hakey, P. M.; Allis, D. G.; Hudson, M. R.; Korter, T. M. Density Functional Dependence in the Theoretical Analysis of the Terahertz Spectrum of the Illicit Drug MDMA (Ecstasy). IEEE Sensors J. 2010, 10, 478−484. (10) Hakey, P. M.; Allis, D. G.; Ouellette, W.; Korter, T. M. Cryogenic Terahertz Spectrum of (+)-Methamphetamine Hydrochloride and Assignment Using Solid-State Density Functional Theory. J. Phys. Chem. A 2009, 113, 5119−5127. (11) Zhang, F.; Kambara, O.; Tominaga, K.; Nishizawa, J.; Sasaki, T.; Wang, H. W.; Hayashi, M. Analysis of Vibrational Spectra of SolidState Adenine and Adenosine in the Terahertz Region. RSC Adv. 2014, 4, 269−278. (12) Zheng, Z. P.; Fan, W. H. First Principles Investigation of LAlanine in Terahertz Region. J. Biol. Phys. 2012, 38, 405−413. (13) Jepsen, P. U.; Clark, S. J. Precise ab Initio Prediction of Terahertz Vibrational Modes in Crystalline Systems. Chem. Phys. Lett. 2007, 442, 275−280. (14) Tonouchi, M. Cutting-Edge Terahertz Technology. Nat. Photonics 2007, 1, 97−105. (15) Kalugin, N. G.; Balandin, A. A. A Special Issue on Terahertz Techniques and Applications. J. Nanoelectron. Optoelectron. 2007, 2, i− ii. (16) Brunner, F. D. J.; Kwon, O. P.; Kwon, S.-J.; Jazbinsek, M.; Schneider, A.; Guenter, P. A Hydrogen-Bonded Organic Nonlinear Optical Crystal for High-Efficiency Terahertz Generation and Detection. Opt. Express 2008, 16, 16496−16508. (17) Schneider, A.; Neis, M.; Stillhart, M.; Ruiz, B.; Khan, R. U. A.; Guenter, P. Generation of Terahertz Pulses through Optical Rectification in Organic DAST Crystals: Theory and Experiment. J. Opt. Soc. Am. B 2006, 23, 1822−1835. (18) Oppenheim, K. C.; Korter, T. M.; Melinger, J. S.; Grischkowsky, D. Solid-State Density Functional Theory Investigation of the Terahertz Spectra of the Structural Isomers 1,2-Dicyanobenzene and 1,3-Dicyanobenzene. J. Phys. Chem. A 2010, 114, 12513−12521. 12606


Article

The Journal of Physical Chemistry C (42) van Eijck, L.; Johnson, M. R.; Kearley, G. J. Intermolecular Interactions in Bithiophene as a Model for Polythiophene. J. Phys. Chem. A 2003, 107, 8980−8984. (43) Ma, J.; Michaelides, A.; Alfe, D.; Schimka, L.; Kresse, G.; Wang, E. G. Adsorption and Diffusion of Water on Graphene from First Principles. Phys. Rev. B 2011, 84, 033402. (44) Hermet, P.; Bantignies, J. L.; Maurin, D.; Sauvajol, J. L. Terahertz Spectroscopy of the Crystalline Alpha-Quaterthiophene: A Combined Experimental and Density Functional Theory Study. Chem. Phys. Lett. 2007, 445, 47−50. (45) King, M. D.; Korter, T. M. Application of London-Type Dispersion Corrections in Solid-State Density Functional Theory for Predicting the Temperature-Dependence of Crystal Structures and Terahertz Spectra. Cryst. Growth Des. 2011, 11, 2006−2010. (46) King, M. D.; Korter, T. M. Noncovalent Interactions between Modified Cytosine and Guanine DNA Base Pair Mimics Investigated by Terahertz Spectroscopy and Solid-State Density Functional Theory. J. Phys. Chem. A 2011, 115, 14391−14396. (47) Juliano, T. R.; Korter, T. M. Terahertz Vibrations of Crystalline Acyclic and Cyclic Diglycine: Benchmarks for London Force Correction Models. J. Phys. Chem. A 2013, 117, 10504−10512. (48) King, M. D.; Korter, T. M. Modified Corrections for London Forces in Solid-State Density Functional Theory Calculations of Structure and Lattice Dynamics of Molecular Crystals. J. Phys. Chem. A 2012, 116, 6927−6934. (49) Delaney, S. P.; Smith, T. M.; Pan, D.; Yin, S. X.; Korter, T. M. Low-Temperature Phase Transition in Crystalline Aripiprazole Leads to an Eighth Polymorph. Cryst. Growth Des. 2014, 14, 5004−5010. (50) Delaney, S. P.; Smith, T. M.; Korter, T. M. Conformation Versus Cohesion in the Relative Stabilities of Gabapentin Polymorphs. RSC Adv. 2014, 4, 855−864. (51) Zhang, F.; Hayashi, M.; Wang, H. W.; Tominaga, K.; Kambara, O.; Nishizawa, J.; Sasaki, T. Terahertz Spectroscopy and Solid-State Density Functional Theory Calculation of Anthracene: Effect of Dispersion Force on the Vibrational Modes. J. Chem. Phys. 2014, 140, 174509. (52) Li, R. Y.; Zeitler, J. A.; Tomerini, D.; Parrott, E. P. J.; Gladden, L. F.; Day, G. M. A Study into the Effect of Subtle Structural Details and Disorder on the Terahertz Spectrum of Crystalline Benzoic Acid. Phys. Chem. Chem. Phys. 2010, 12, 5329−5340. (53) Takahashi, M.; Kawazoe, Y.; Ishikawa, Y.; Ito, H. Interpretation of Temperature-Dependent Low Frequency Vibrational Spectrum of Solid-State Benzoic Acid Dimer. Chem. Phys. Lett. 2009, 479, 211− 217. (54) Delaney, S. P.; Pan, D. H.; Galella, M.; Yin, S. X.; Korter, T. M. Understanding the Origins of Conformational Disorder in the Crystalline Polymorphs of Irbesartan. Cryst. Growth Des. 2012, 12, 5017−5024. (55) Nickel, D. V.; Delaney, S. P.; Bian, H. T.; Zheng, J. R.; Korter, T. M.; Mittleman, D. M. Terahertz Vibrational Modes of the Rigid Crystal Phase of Succinonitrile. J. Phys. Chem. A 2014, 118, 2442− 2446. (56) Tanno, T.; Umeno, K.; Ide, E.; Katsumata, I.; Fujiwara, K.; Ogawa, N. Terahertz Spectra of 1-Cyanoadamantane in the Orientationally Ordered and Disordered Phases. Philos. Mag. Lett. 2014, 94, 25−29. (57) Korter, T. M.; Balu, R.; Campbell, M. B.; Beard, M. C.; Gregurick, S. K.; Heilweil, E. J. Terahertz Spectroscopy of Solid Serine and Cysteine. Chem. Phys. Lett. 2006, 418, 65−70. (58) Kim, J.-T., Computational Study of Organic Nonlinear Optical (NLO) Materials for Terahertz (THz) Applications: Molecule, Cluster, and Solid-State Approaches. Ph.D Dissertation, Korea Advanced Institute of Science and Technology, February 2012.

12607


First-Principles Calculation of Terahertz Absorption with Dispersion

Recommend Documents