Next Article in Journal
Convenient Synthesis of 6,7,12,13-Tetrahydro-5H-Cyclohepta[2,1-b:3,4-b’]diindole Derivatives Mediated by Hypervalent Iodine (III) Reagent
Next Article in Special Issue
Ethylene Measurements from Sweet Fruits Flowers Using Photoacoustic Spectroscopy
Previous Article in Journal
The Oligomeric State of the Plasma Membrane H+-ATPase from Kluyveromyces lactis
Previous Article in Special Issue
Differentiation between Enamines and Tautomerizable Imines Oxidation Reaction Mechanism using Electron-Vibration-Vibration Two Dimensional Infrared Spectroscopy
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

DFT Computed Dielectric Response and THz Spectra of Organic Co-Crystals and Their Constituent Components

by
Joseph W. Bennett
,
Michaella E. Raglione
,
Shalisa M. Oburn
,
Leonard R. MacGillivray
,
Mark A. Arnold
and
Sara E. Mason
*
Department of Chemistry, University of Iowa, Iowa City, IA 52242, USA
*
Author to whom correspondence should be addressed.
Molecules 2019, 24(5), 959; https://doi.org/10.3390/molecules24050959
Submission received: 18 January 2019 / Revised: 26 February 2019 / Accepted: 4 March 2019 / Published: 8 March 2019
(This article belongs to the Special Issue Vibrational Probes of Biomolecular Structure and Dynamics)

Abstract

:
Terahertz (THz) spectroscopy has been put forth as a non-contact, analytical probe to characterize the intermolecular interactions of biologically active molecules, specifically as a way to understand, better develop, and use active pharmaceutical ingredients. An obstacle towards fully utilizing this technique as a probe is the need to couple features in the THz regions to specific vibrational modes and interactions. One solution is to use density functional theory (DFT) methods to assign specific vibrational modes to signals in the THz region, coupling atomistic insights to spectral features. Here, we use open source planewave DFT packages that employ ultrasoft pseudopotentials to assess the infrared (IR) response of organic compounds and complex co-crystal formulations in the solid state, with and without dispersion corrections. We compare our DFT computed lattice parameters and vibrational modes to experiment and comment on how to improve the agreement between theory and modeling to allow for THz spectroscopy to be used as an analytical probe in complex biologically relevant systems.

1. Introduction

The medical industry has begun to utilize active pharmaceutical ingredients (APIs) in a co-crystalline form to enhance the functionalities of their products [1,2,3]. These solids are unique, as they contain two non-ionic components, an API and a coformer. The API and coformer can bind through van der Waals forces and non-covalent interactions such as hydrogen bonding, similar to biologics such as DNA and proteins [4,5]. Binding through non-ionic forces can enable more flexibility in crystalline engineering, removing the need to form an ionic bond. Understanding the intermolecular forces holding these co-crystals together can help to better understand, predict and control their physical properties. The creation of co-crystal APIs is a strategy that has shown to influence solubility and bioavailability [6], thermal stability [7], and has also found use in the control of product concentrations in the development of new types of therapeutics [8].
To better understand the structures of co-crystals with APIs, researchers have utilized terahertz (THz) spectroscopy and low frequency Raman spectroscopy (LFRS). Both LFRS and THz spectroscopy probe frequencies between 3–300 cm 1 , where the prominent vibrational modes are related to the packing of the constituent organics in a crystal and the intermolecular forces that hold them in place [9,10]. These spectroscopic techniques have already proven useful for identifying different polymorphs that form when mixing multiple organic compounds [10,11,12], for the in situ monitoring of crystallization and structural transformations [13,14], as well as uncovering solid state transition mechanisms [15,16]. THz spectroscopy is potentially a powerful tool in probing biological systems, because it is compatible with online measurements, remote sampling, and three-dimensional imaging [17], although little work has been published on large complex systems beyond binary and tertiary co-crystals. Most THz research has focused on smaller organic-based systems (≈50–150 atoms/unit cell) where several groups have worked towards using theory and modeling to assign the vibrational modes found in both THz absorption and low frequency Raman spectroscopies [10,18,19,20].
An open question in the theory and modeling of materials is: “How closely can calculated properties match experimentally determined information?” This is especially difficult to answer when computing the measurable properties of a complex material with multiple components and potentially nonlinear behavior, or if changes in the chemical environment, such as the formation of hydrogen bonds or solvation, must be taken into account. The atomistic detail needed to address these issues can be gleaned from using methods based on quantum mechanics, specifically density functional theory (DFT). DFT methods have been around since the mid-1960s [21,22] and allowed for calculation of properties in the solid state for simple systems, such as the ordering and electronic structure of Si [23] and ZnS semiconductors [24]. Since the 1990s improvements in pseudopotential design and construction [25,26,27,28,29,30] and exchange-correlation functionals [31,32,33], as well as the distribution of open source software packages [34,35] and the advent of highly parallelized computing have placed DFT methods at the forefront of accurate and reliable methods based in quantum mechanics.
How organic molecules pack into crystal structures is ultimately governed by intermolecular interaction patterns. While chemists can qualitatively identify systems in which, for example, π π stacking, hydrogen bonding, and long-range dispersion forces occur, it is extremely difficult to predict from such an analysis what crystal structure will be preferred for a given organic molecule, either individually or with a co-former. Here, first-principles modeling based on DFT calculations becomes attractive as a method known to be reliable for modeling the structure and properties of inorganic crystal structure materials [36,37,38,39,40,41].
DFT calculations are upheld by materials chemistry as the gold-standard in quantum mechanical modeling owing to desirable balance between computational speed and accuracy. However, directly applying these methods to organic crystalline materials is complicated by the differing intermolecular forces active in these materials, noted above. While DFT is in principle an exact method, in practice the local density approximation, LDA (local density approximation) (and extensions to it, such as the generalized-gradient approximation, GGA) to the true, universal exchange-correlation functional are employed. In effect, this means that, while all non-relativistic electronic effects are included in the formal exact DFT functional, including van der Waals interactions, the long-range, non-local correlations that contribute to van der Waals forces are not accounted for in calculations that use LDA or GGA functionals. Solutions to this problem are complicated, as local effects are, by construction, included. Therefore, extending DFT approximations to include long-range forces is not as simple as just adding additional terms. Indeed, while such efforts to extend DFT came to the forefront in the early-2000s, the methods themselves [42,43,44,45,46,47], as well as best-practices for their applications [48,49], are still active areas of research.
Of the various approaches to extend DFT to include long-range electron correlation interactions, D2 (and more recently D3) dispersion coefficient methods, [43,46,50], referred to as DFT-D, are popular based on the relatively simple mathematical form, improved predictions, and use of adjustable parameters. While it is beyond the scope of this manuscript to review the methods (those this has been done elsewhere, for example, Ref. [51]), it can be stated briefly that DFT-D empirical parameters are determined by fits done to conformational and interaction energies computed using high-level quantum mechanical calculations. While in some implementations the fit parameters are adjustable, this is not universally the case.
Here, we present calculations that use open source DFT packages (that employ pseudopotentials and a planewave basis set and include dispersion corrections) to obtain the lattice parameters and vibrational modes of organic co-crystals and their constituent components. It is impossible to know which empirically tuned variant of DFT-D will offer the best performance for the organic co-crystals under study here. Given the lack of benchmarking information in the literature, we elect to test the performance of Grimme-D2 DFT-D calculations, compared to standard DFT-GGA. Our choice of Grimme-D2 DFT-D for this study is motivated by the fact that, in some implementations, the empirical parameter controlling the dispersion corrections is tunable by the user, enabling us to explore how its variation affects predictions of the bulk lattice constants. Future studies could go on to test other DFT-D functionals, and this would contribute to the developing understanding of which correction schemes offer optimal performance for organic crystals. We compare our DFT calculated results, both with and without dispersion corrections, to experimentally determined lattice constants and absorptions over THz frequencies. We use the information obtained with DFT to analyze the dielectric response of the systems studied here in order to elucidate trends in the types of functional groups present, crystal packing, and ultimately how one could use dielectric response as a tailorable property when designing and understanding complex organic co-crystal API formulations. We focus on identifying infrared (IR) active modes, their contributions to the overall dielectric response, and how closely they match the IR active modes measured in THz absorption spectroscopy. We use the work presented here to establish a baseline for future studies where recently developed methods to model organic solids could be used to improve upon the disagreements between theory and experimental endeavors.

2. Results and Discussion

The fully relaxed monoclinic crystal structures of the components used to create the co-crystals (discussed in the next subsection) are depicted in Figure 1. All three are indexed in space group 14, where all three have monoclinic angle β > 90 . Shown are the repeat units of 1,2-bis(4-pyridyl)ethylene (BPE), 1,2-bis(4-pyridyl)ethane (BPEth), and salicylic acid (SA) where the direction containing the greatest dispersion is marked by the vertical axis. The coformers BPE and BPEth differ in the functional group connecting the planar C 5 N ring units; inspection of Figure 1a,b shows that BPE contains an alkene C 2 H 2 unit and BPEth contains an alkane C 2 H 4 unit, and that this influences the relative orientation of packing of these structures. In neither case is there a significant number of intermolecular (primary) hydrogen bonds, and a consequence of this is packing directions in which dispersion is significant. Table 1 contains the computed crystallographic parameters of the fully relaxed systems. a, b, and c are the lattice constants of the primitive unit cell, and β is the angle between a and c.
Inspection of Table 1 shows that DFT without dispersion corrections overestimated BPE lattice constants a and c by ≈ 40 and 20 %, and BPEth lattice constant b by ≈ 20%, relative to experiments. Inclusion of dispersion (DFT-D) reduces these errors to ≈−7 and −3 %, and ≈−4 %, for BPE and BPEth, respectively. Another error that is not typically reported in DFT studies is that of β , the angle between a and c. For the monoclinic BPE unit cell, DFT overestimates this angle by over 25%, and this is because the errors in a and c are so large. DFT-D decreases the error in β to ≈−2.5%. The experimentally determined lattice parameters for BPE [52], BPEth [53], and SA [54] are given in the parenthesis of Table 1, followed by the % deviation of each lattice parameter from experiments. The experimental structures in Refs. [52,53,54] were determined using X-ray diffraction at room temperature.
The SA units in Figure 1c display a packing that is influenced by hydrogen bonding between the carboxyls of each unit (mainly along y, but also present in x) and this results in a network of packed dimeric units where dispersion is greatest along z. DFT overestimates the lattice constants a and c by ≈ 6 and 32%, and this is improved using DFT-D; the difference between DFT-D computed lattice constants and experiment is reduced to −0.8 and −0.2% for a and c, respectively. Much like BPE, DFT-D also improves the agreement between computed and experimentally determined angle β ; inclusion of dispersion changes β from 75.39 to 93.54 , more in line with experiment.
Beyond agreement of lattice parameters, inclusion of dispersion also affects the vibrational modes of a material. We now focus on comparing DFT and DFT-D computed vibrational modes, as well as their relative contributions to the dielectric constant, ϵ , for each of the three co-crystal components, BPE, BPEth, and SA. Table 2 reports the isotropically averaged values of ϵ for DFT and DFT-D, compared to the experimentally determined values. We find that the DFT-D values are closer to experiment than DFT values, where for BPE and BPEth DFT underestimates ϵ by −22 and −14%, respectively, while for DFT-D the values are overestimated by +14 and 4%. The opposite trend is observed for SA, where DFT overestimates ϵ by −20%, and DFT-D underestimates it by 10%. Both experiment and DFT-D show that ϵ increases in the order BPE > BPEth > SA.
We now turn to an analysis of the frequencies and oscillator strength of the DFT and DFT-D computed vibrational modes that govern ϵ . Details of how to compute ϵ from first-principles using this information are given in the Materials and Methods subsection called Computational Details.
The experimentally determined THz spectra of BPE contains two peaks centered at 59.4 and 92.73 cm 1 , as depicted by the blue curves in Figure 2a,b. Figure 2 is a comparison of the vibrational modes determined by an analysis of the dielectric response, as discussed in the Methodology section, between DFT (a) and DFT-D (b). We find that, in the range of 20–110 cm 1 , DFT predicts nine total IR active modes with a contribution to ϵ of more than 0.01, and that the peak heights (scaled by contribution to ϵ in this region) are very much not in alignment. Table 3 shows that the largest contribution to ϵ in DFT is a mode at 31.1 cm 1 , which accounts for more than 50% of the total response along y. This is not the case using DFT-D, which predicts four total IR active modes with a contribution larger than 0.01 to ϵ in the range of 20–110 cm 1 . These four peaks occur as two sets of peaks at 67.6/69.5 cm 1 and 105.6/106.4 cm 1 , where they are within 1–2 cm 1 of each other and most likely not resolvable in experiments. Figure 2b shows that DFT-D with an s 6 = 0.75 obtains two sets of modes that are shifted by +10 and +13 cm 1 from the experimentally determined spectra.
We are able to differentiate the peaks computed by DFT-D into contributions to the x (a), y (b), and z (c) directions of ϵ , as shown in Table 3 and focus on the modes at 105.6 and 106.4 cm 1 , which contribute ≈ 33 and 50% of the overall response in x- and y-directions, respectively. These vibrational modes are depicted in Figure 3 and can be described as asymmetric stretches where the largest displacements are localized on nitrogen and the carbon in the C 2 H 2 unit, where these atoms displace in opposite directions.
While it is helpful to separate the components of ϵ by direction for an analysis of the vibrational modes, most experiments will report an isotropically averaged value, as discussed in the Methodology section. The DFT and DFT-D computed isotropically averaged ϵ μ (ionic response) of BPE are 0.27 and 0.29, while the directionally averaged ϵ α (electronic response) are 2.41 and 3.63. Inclusion of dispersion in DFT-D has increased ϵ μ by 7.4% and ϵ α by 50.6%. Summed together, the averaged components of ϵ for BPE using DFT and DFT-D are 2.68 and 3.92, respectively.
The experimentally determined THz spectra of BPEth contains three peaks centered at 58.99, 66.26, and 91.92 cm 1 , as depicted by the blue curves in Figure 4a,b. Figure 4 is a comparison of the vibrational modes determined by an analysis of the dielectric response, as discussed in the Methodology section, between DFT (a) and DFT-D (b). We find that in the range of 20–110 cm 1 , DFT predicts seven total IR active modes with a contribution to ϵ of more than 0.01, and that the peak heights (scaled by contribution to ϵ in this region) are not in alignment, much like for BPE.
Table 4 shows that the largest contributions to ϵ in DFT are three modes below 42 cm 1 (at 31.8, 32.5, and 41.6 cm 1 ) and two modes at 68.8 and 70.8 cm 1 , where all five vibrational modes have a comparable contribution to ϵ . Inclusion of dispersion in DFT-D decreases the number of IR active vibrational modes to six and changes the distribution to higher frequencies, as well as their respective contribution to ϵ . This case is depicted in Figure 4b. For example, the lowest frequency peak is at 45.4 cm 1 , and the modes with the largest contribution to ϵ occur at 61.2 and 70.9 cm 1 , in line with the experimentally determined spectra. Also of note in Figure 4b is the appearance of two peaks at 93.8 and 108.9 cm 1 in the DFT-D computed spectra, which were absent in the DFT calculation without dispersion, and coincidental with peaks in the experimentally determined THz spectra. While the DFT-D spectra of BPEth has more inconsistencies when compared to the experimental THz spectra than BPE, the frequency distribution and contribution is better than when using only DFT.
The vibrational modes of BPEth at 61.2 and 70.9 cm 1 , which are the two largest contributors to ϵ μ in the y- and x-directions, respectively, are shown in Figure 5. The vibrational mode at 61.2 cm 1 can be described as an asymmetric stretch in which the largest displacements are localized on the nitrogen in the aromatic ring and the carbon in the C 2 H 4 unit connected the two aromatic rings. Unlike the mode identified for BPE, the motion of these atoms are in the same direction, and not opposed. The vibrational mode at 70.9 cm 1 has displacements in the x z -direction and can be described as an asymmetric stretching mode mixed with some aromatic ring rotation. In this mode, some of the aromatic carbon displace just as far as the nitrogen and the carbon in the C 2 H 4 unit.
The DFT and DFT-D computed isotropically averaged ϵ μ (ionic response) of BPEth are 0.32 and 0.27, while the directionally averaged ϵ α (electronic response) are 2.36 and 2.98. Inclusion of dispersion in DFT-D has decreased ϵ μ by 15.6% and increased ϵ α by 26.3%. Summed together, the averaged components of ϵ for BPEth using DFT and DFT-D are 2.68 and 3.25, respectively, comparable to the values of 2.68 and 3.92 for BPE.
Turning from the nitrogen-containing heterocycles BPE and BPEth to the carboxylic acid component of the co-crystals, we investigate the dielectric response of salicylic acid. As depicted in Figure 1c, the unit cell for SA has hydrogen bonding that takes place primarily along the y-axis, and this manifests as an increased contribution to ϵ μ along y relative to the other two axes for both DFT and DFT-D. Table 5 shows that for DFT the y component of ϵ μ is 1.01, while for x and z the sums are 0.79, and 0.44, respectively. Adding dispersion changes, the magnitudes of these three values, but not the trend; the trend in ϵ μ is still y > x > z.
Figure 6 shows that neither DFT or DFT-D computed spectra are a good match to the experimentally determined THz spectra of SA. The seven experimentally determined THz peaks of SA are centered at 36.8, 46.1, 53.5, 60.0, 68.9, 77.0, and 89.3 cm 1 in the range of 20–110 cm 1 . DFT computes too many vibrational modes at lower frequencies (30–50 cm 1 ), resulting in too many total modes, and no modes above 80 cm 1 . DFT-D has relative peaks whose contributions to ϵ are not aligned to the experimental structure, but computed the correct number of absorption features. If we compare them to the seven DFT-D computed peaks at 37.9, 49.7, 66.5, 76.1, 78.4, 92.3, and 100.2 cm 1 , then the differences in frequency are 1.1, 3.6, 13.0, 16.1, 9.5, 15.3, 10.9 cm 1 going from 20 to 110 cm 1 .
We identify the three vibrational modes of SA that contribute a majority of the dielectric response, which occur at 37.9, 76.1, and 92.3 cm 1 and are shown in Figure 7. All three modes have significant contributions from the oxygen atoms of both the phenolic OH and COOH acid functional groups. At 37.9 cm 1 , the OH of the COOH acid group displaces along y, opposite to OH on the aromatic ring, while the carbon in the aromatic ring breathes in and out. Higher in frequency, at 76.1 cm 1 , the oxygen of the COOH carbonyl and the OH on the aromatic ring displace in the same direction along x z while carbon in the aromatic ring undergo asymmetric stretches. The highest frequency in Figure 7, 92.3 1 , has the OH of the COOH acid group moving along y, opposed to the OH on the aromatic ring. In this vibrational mode, the carbon of CO functional group also moves, and overall the mode represents a partial rotation of SA.
For SA, we find that the difference between DFT and DFT-D computed contributions to ϵ are opposite that of BPE and BPEth; both ϵ μ and ϵ decrease significantly with the addition of dispersion. ϵ μ decreases from 0.75 to 0.51, and ϵ decreases from 2.95 to 2.28. This means that the overall averaged ϵ goes from 3.70 to 2.79, a decrease of 24.6% when including dispersion corrections. The directionally averaged value of ϵ (0.51) is ≈20% of the DFT-D calculated dielectric response.
In general, all three comparisons between DFT and DFT-D in Figure 2, Figure 4, and Figure 6 show that DFT will yield high-intensity, low frequency modes in the range of ≈25–40 cm 1 whose intensity and frequency change with the application of dispersion corrections in DFT-D. While the agreement between DFT-D and experimentally determined THz spectra is better for BPE than for BPEth and SA, DFT-D does yield results in which the three components are differentiable; BPE has IR active modes at 106 cm 1 , BPEth has modes in the range of 60–70 cm 1 , and SA has modes that span the range of 40–90 cm 1 , where the modes at 38, 76, and 92 cm 1 all involve displacements of the oxygen in both the OH and COOH functional groups. The numerical values of dielectric response of BPE and BPEth are more similar to each other than to SA, and this may be because the only difference between them is the alkene (C 2 H 2 ) vs alkane (C 2 H 4 ) connection between heterocyclic rings. The contributions to ϵ μ are greatest for SA, which may be explained by the increased hydrogen bonding of SA relative to BPE and BPEth.

2.1. Co-Crystals

Here, we perform the same types of analyses as in the previous section but for the larger co-crystal systems 2(SA)·BPE Forms I and II, and 2(SA)·BPEth. As will be discussed elsewhere, all three form monoclinic crystal structures (space group 14), but 2(SA)·BPE-II and 2(SA)·BPEth display packing similar to each other, while 2(SA)·BPE-I is different. The synthesis and characterization of these three co-crystals systems will be described elsewhere. While all three co-crystal structures display hydrogen bonding between the carboxyl unit of SA and the nitrogen of BPE and BPEth, 2(SA)·BPE-II and 2(SA)·BPEth also have the H of the aromatic OH of SA in close proximity to the carbonyl O of the COO unit creating a stable six member ring via resonance. This additional hybridization is absent in 2(SA)·BPE-I. The DFT-D relaxed structures are shown in Figure 8 for (a) 2(SA)·BPE-I, (b) 2(SA)·BPE-II and (c) 2(SA)·BPEth, looking down the y-axis at the x z plane.
We relax the co-crystal compositions using both DFT and DFT-D, and the lattice parameters of these relaxations are tabulated in Table 6. Much like the components BPE, BPEth, and SA discussed in the previous section of the Results, relaxations using DFT yield deviations from experimentally determined structures on the order of 5–25% overestimation, and the largest errors coincide with the directions where dispersion will be the greatest. This is observed for lattice constant a of 2(SA)·BPE-I and lattice constant c of 2(SA)·BPE-II and 2(SA)·BPEth. Inclusion of dispersion corrections in DFT-D decreases the error of the lattice parameters, in most cases to within the 1–2% acceptable for DFT-GGA calculations in the solid state, but DFT-D underestimates lattice constant c of 2(SA)·BPE-II and 2(SA)·BPEth by ≈ 8%.
As the disagreement between experimentally determined and DFT-computed lattice parameters with no dispersion correction is significantly large, and we showed earlier that DFT-D computed vibrational modes of the components were needed for a better description of the THz spectral features, we compute the vibrational modes of the co-crystals using DFT-D. Figure 9 depicts the DFT-D computed vibrational modes of the co-crystals, compared to experimentally determined THz spectra. The comparison of vibrational modes in Figure 9a shows that DFT-D predicts low frequency modes for 2(SA)·BPE-I (below 40 cm 1 ) close to experiment. Experiment shows ten THz peaks of 2(SA)·BPE-I centered at 29.5, 33.5, 38.8, 44.9, 54.6, 68.1, 73.7, 92.7, 99.6, and 102.4 cm 1 , compared to DFT-D which predicts IR active modes at 28.9, 33.9, 34.9, 48.5, 50.9/51.1, 61.5, 69.8/70.3, 87.8, 109.5 cm 1 . The smallest difference in frequencies between THz spectroscopy and DFT-D is on the order of 1–2 cm 1 and the largest difference in frequencies is on the order of 10 cm 1 . The minor deviations in vibrational mode frequencies, from the experimentally determined values, may be related to the close agreement between experimentally determined and DFT-D computed lattice parameters, which are on the order of ± 1.5%. The experimentally determined THz spectra of 2(SA)·BPE-II (Figure 9b) and 2(SA)·BPEth (Figure 9c) are similar with frequencies centered at 35.2, 43.6, 61.6, 83.4, and 110.7 cm 1 for 2(SA)·BPEth and at 24.9, 35.4, 42.8, 60.8, 84.0, and 108.7 cm 1 for 2(SA)·BPE-II. The experimentally determined frequencies are listed in Table 7.
We compare both the experimentally determined and DFT-D computed spectra of the three co-crystals to each other in Figure 10. When comparing the experimental spectra (top of Figure 10), it seems that orientation and packing effects dominate the THz moreso than the difference in ligand identity as BPE and BPEth. The comparison of 2(SA)·BPE-I and 2(SA)·BPE-II in Figure 10 (left) highlights the difference between co-crystals with the same components, but different crystal packing and, thus, different lattice parameters. The experimental spectra show coincidental peaks only at ≈35 and 43 cm 1 , while the DFT-D do not coincide at low frequency (30–60 cm 1 ), coincide in the range of 70–90 cm 1 , and begin to differ again above 90 cm 1 .
For 2(SA)·BPE-II and 2(SA)·BPEth, the DFT-D calculated IR active modes are consistently shifted higher, at least 5–8 cm 1 than experiments, if they match up at all. This mismatch between the THz and DFT-D for 2(SA)·BPEth and 2(SA)·BPE-II co-crystals may be caused by the ≈ 8% underestimation of lattice constant c. Figure 10 (right) shows that the experimental THz spectra of 2(SA)·BPE-II and 2(SA)·BPEth coincide at all frequencies above 30 cm 1 and that for the DFT-D computed spectra they (roughly) coincide at low frequencies (30–60 cm 1 ), do not coincide in the range of 70–90 cm 1 , and begin to coincide again above 90 cm 1 .
We compare the three ranges of frequencies computed via DFT-D to determine if our methods could yield spectral resolution that goes beyond the resolution capabilities of experimental THz data collection at room temperature, and our analysis points towards characteristic frequency ranges for these co-crystal systems that may indicate where differences in packing vs. differences in bipyridyl coformer (as ligand identity) will dictate vibrational modes. We find that packing arrangement may dominate between 30–60 cm 1 , ligand substitutions may dominate in the range of 70–90 cm 1 , and then packing may dominate again above 90 cm 1 . Example vibrational modes computed using DFT-D in these three regions are shown in Figure 11, Figure 12 and Figure 13 for 2(SA)·BPEth, 2(SA)·BPE-I, and 2(SA)·BPE-II, respectively. The differences between the the experimentally determined and DFT-D computed vibrational modes motivate the need for further investigations that could include low temperature THz spectra collection and methodology that goes beyond the Grimme D2 formalism used here, which will be discussed in the next section.
The DFT-D calculated contributions to ϵ μ for 2(SA)·BPEth are tabulated in Table 8 and Table 9 for 2(SA)·BPE polymorphs I and II. The directionally averaged values of ϵ μ and ϵ for 2(SA)·BPEth are 1.14 and 2.97, 1.48 and 3.05 for 2(SA)·BPE-I, and 0.81 and 3.33 for 2(SA)·BPE-II, and when added together yield ϵ that are 4.11, 4.53, and 4.14, respectively. ϵ μ for the co-crystals is 20–30% of the total dielectric response, an increase when compared to the behaviors of their components, such as SA (18%) and BPEth (8%). The isotropically averaged values of ϵ for the three co-crystals are tabulated in Table 10, compared to the experimentally determined values. We find that all three DFT-D computed ϵ are underestimated compared to experiments, and that, unlike the ϵ of the components, the trends in magnitude of ϵ do not agree between experiments and DFT-D. Moreover, the DFT-D calculated values are closer to each other than those measured experimentally.
The overall dielectric responses of the co-crystals and their components also demonstrate that intermolecular interactions, highly dependent on packing orientation, result in a dielectric response that is not merely an average of the components. For example, the directionally averaged values of ϵ μ for SA and BPE are 0.51 and 0.29, respectively, and an average of those two numbers is 0.40. The values for 2(SA)·BPE-I and 2(SA)·BPE-II are 1.48 and 0.81, which is more than even when the ϵ μ are added together (0.79). This implies a tunable, nonlinear dielectric response can be achieved, and potentially optimized, in co-crystal formulations.
Analysis of the vibrational modes of the co-crystals show that all modes in the range of 20–110 cm 1 have significant contributions from both the nitrogen containing heterocycles BPE and BPEth, and SA. A majority of the largest displacements in each mode above 40 cm 1 is localized on both of the O atoms of the COOH acid, and the phenolic OH of SA, and in many cases the nitrogen of the heterocyclic ring as well. Below 40 cm 1 , such as the modes of 2(SA)·BPE-I in Figure 12, O motion contributes less than the motions present in the heterocyclic ring. Vibrational modes with frequencies higher than 90 cm 1 , such as the modes of 2(SA)·BPE-II in Figure 13 also tend to have larger carbon displacements than modes at lower frequencies.

2.2. Improving Agreement with Experiment

In the two previous sections of the Results and Discussion, we showed that using DFT-D to compute the lattice parameters and vibrational spectra of the 2(SA)·BPE and 2(SA)·BPEth co-crystals, as well as their components, results in a marked improvement over the same parameters computed using DFT. There still remains some inconsistencies between DFT-D computed spectra and experimentally determined lattice parameters and THz spectra, such as 8% lattice constant underestimation for 2(SA)·BPE-II and 2(SA)·BPEth and the upshift in frequency for the IR active modes of BPE.
One method to improve agreement between experimental THz spectra and DFT-D modeling is to adjust the global dispersion scaling factor, s 6 , from the default value of 0.75. Systematically changing s 6 results in changes in lattice parameters and frequencies of IR-active modes, where better agreement can be reached between modeling efforts and THz spectra [4,18,19,55]. Here, we turn to another open source planewave DFT code in which varying s 6 is allowed in the input file. We use Quantum Espresso [35] and the GBRV (Garrity-Bennett-Rabe-Vanderbilt) [29] set of ultrasoft pseudopotentials [27] to map out how changing the global dispersion correction will affect lattice parameters for BPE, BPEth, SA, and 2(SA)·BPEth. We keep the k-grid sampling and energy convergence criteria the same as for the calculations employing ABINIT [34] and the ONCV-type (optimized norm-conserving Vanderbilt) of pseudopotential [28].
Figure 14 shows that s 6 values of 0.4, 0.5, 0.6, and 0.35 for BPE (a), BPEth (b), SA (c), and 2(SA)·BPEth (d) independently minimize the percent error when lattice parameters are compared to experiment. If we increase the maximum % error to ± 4%, then we obtain ranges of s 6 values of 0.25–0.54 for BPE, 0.33–0.55 for BPEth, 0.42–0.95 for SA, and 0.22–0.54 for 2(SA)·BPEth. In all cases, but SA, an s 6 = 0.75 is not present in any of these ranges. Using these four test cases as an example, we obtain a range of 0.42 < s 6 < 0.54, which is in line with the values used in previous studies of organic co-crystals [4,18,19,55].

3. Materials and Methods

3.1. Pellet Fabrication

Pellets were prepared by the methods described before [56] for the single components as well as the co-crystals, of which the details of syntheses and complete structural characterization will be detailed elsewhere. Briefly, a commercial PTFE (Polytetrafluoroethylene) powder, with particle diameters ranging from 9–13 μ m, was purchased from Micro Powders (FLUO 625 CTX2, Micro Powders, INC., New York, NY, USA). The powder was dried at 60 C and stored in a desiccator before use. For each analyte, a mixture was prepared by co-grinding 40 mg of analyte with 1700 mg of dried PTFE for five minutes in an agate mortar and pestle. Three sample pellets were prepared from each mixture by placing 400–450 mg of the mixture into a 13 mm diameter stainless steel die and using a Specac hydraulic press (model number 15011, Kent, UK) to apply a 5-ton load (0.34 GPa) for 5 min. Freshly formed sample pellets were removed from the die and placed in a desiccator until spectra were collected.

3.2. Dielectric Measurements

Dielectric constants were extracted as described previously [57]. Briefly, values of refractive index ( η ) are extracted from phase information embedded within the time-domain spectra.
The refractive index of the pellet is related to the difference in phases for the sample pellet ( ϕ s ) and reference air ( ϕ r ), according to Equation (1), where c is the speed of light, b is the sample thickness, or path length, and ω is the frequency of the electromagnetic radiation in Hz. The dielectric constant for the sample can be obtained according to Equation (2), where ϵ ( ν ) represents the dielectric constant as a function of frequency in wavenumber and k( ν ) is the frequency dependent extinction coefficient for the pellet material. Generally, k( ν ) is negligible [58] and the dielectric constant can be taken as the square of the refractive index:
η ( ω ) = 1 c 2 π ω b ( ϕ s ( ω ) ϕ r ( ω ) ) ,
ϵ ( ν ) = η ( ν ) 2 + k ( ν ) 2 .
As the pellets are mostly composed of PTFE, the analyte dielectric must be extracted using the Landau, Lifshitz and Looyenga (LLL) model [59]. This model is defined in Equation (3) using the Looyenga power law, specifically for our three-component system. In Equation (3), the dielectric constant of the pellet is composed of a volume fraction (v of each component; PTFE, sample, and air). In this work, the volume of air is obtained by subtracting the total volume of the pellet from the sum of volumes of the crystals and PTFE. The volumes of both the PTFE and crystals were determined from the mass of each and their known densities. For the crystalline samples examined here, the crystal structure density was used and, for PTFE, the value of 2.26 g/cm 3 was used for its density. The density of PTFE was measured using the volume and mass of a pure pellet:
ϵ 1 / 3 = ϵ air 1 / 3 v air + ϵ PTFE 1 / 3 v PTFE + ϵ analyte 1 / 3 v analyte .

3.3. Terahertz Spectroscopy

THz transmission spectra were collected and analyzed by the methods described before [56] using with a Teraview TPS 1000D time-domain terahertz spectrometer (TeraView Limited, Cambridge, UK). For each sample, three pellets were tested by rotating the samples into the sample holder three times and taking three triplicate measurements. This resulted in a total of 27 time-domain spectra for each analyte. Between each consecutive measurement, an air background was collected. Each spectrum was collected as 1800 co-added scans attained over one minute. Purging the sample compartment with dried air avoided the presence of confounding water vapor lines.
Time-domain spectra were truncated just prior to the etalon feature and zero-filled to 8192 (2 13 ) points. The truncated data were treated with a boxcar apodization function followed by the Fourier transformation to yield the corresponding frequency-domain electric field spectrum. Absorbance spectra were then calculated as twice the negative base ten logarithm of the ratio of the sample to air electric field spectra [60]. Twice the negative base ten logarithm is required to square the electric field values in order to realize intensities. The resolution of the resulting spectra was 1.2 cm 1 over a spectral range of 10–110 cm 1 .

3.4. Computational Details

Density functional theory (DFT) calculations with periodic repeat boundary conditions were carried out using the ABINIT open source software package [34], unless otherwise noted. The generalized gradient approximation (GGA) of Perdew, Burke and Ernzerhof (PBE) [32] was used as the exchange-correlation functional for all calculations, and DFT structural optimizations that included dispersion corrections used the Grimme-D2 implementation [43]. All structural optimizations employed a variable cell relaxation where all lattice constants, lattice angles, and Wyckoff positions were optimized concurrently starting from the experimentally determined monoclinic crystals, and monoclinic symmetry was maintained throughout all structural relaxations. We differentiate between the two types of calculations as DFT (with no dispersion correction) and DFT-D (which includes dispersion correction). For all DFT-D calculations, the functional dependent scaling factor, s 6 , was 0.75, unless otherwise noted.
All atoms were represented as optimized norm-conserving Vanderbilt (ONCV) pseudopotentials [28], that were chosen because of their accuracy and computational efficiency [30]. All calculations used a plane wave cutoff of 40 Ry, and all atoms were allowed to fully relax during structural optimizations. The convergence criteria for structural optimizations was a maximum residual force of 5 meV/Angstrom per atom. The input structures for the dispersion corrected DFT calculations will be detailed elsewhere. A 4 ×4 × 4 k-point mesh [61] was used for all structural optimizations because the difference in total energy between a 4 × 4 × 4 and 6 × 6 × 6 k-point mesh was below 3 meV per co-former unit.
Here, we follow the methodology outlined for first-principles calculations of static dielectric properties [62,63] as described in previous work on solid state bulk materials [64,65,66]. Briefly, once a structure is fully relaxed using either DFT or DFT-D, response function calculations [67,68] are carried out to generate D α β ( i , j ) , the mass weighted dynamical matrix:
D α β ( i , j ) = 2 E m i m j τ i α τ j β .
In Equation (4), E is the DFT (or DFT-D) computed total energy, m i ( m j ) is the mass of atom i (j), and τ i α is the displacement of atom i (j) in direction α ( β ). This creates a square matrix of dimension 3N×3N, where N is the total number of atoms, and each eigenvalue of matrix D is ω μ 2 , the frequency squared of a normalized eigenvector a μ . The means that the square root of the eigenvalues are the frequencies that we report in cm 1 , and its eigenvector is a vibrational mode. For each atom, the Born effective charge tensor, Z i α β * , is also calculated and used to compute the effective charge for each mode μ . The mode effective charge, Z μ α * , is defined as:
Z μ α * = i β Z i α β * ( a μ ) i β m i .
Contributions to the static dielectric response come only from vibrational modes that are IR active, and thus have a non-zero Z μ α * . It is these modes that are used to calculate ϵ μ , the ionic contribution to the dielectric constant from mode μ using the following equation:
ϵ μ α β = Z μ α * Z μ β * 4 π 2 ϵ 0 V ω μ 2 ,
where ϵ 0 is the permittivity of free space and V is the volume of the bulk solid. The total (isotropically averaged) ionic contribution to the dielectric constant is therefore given by:
ϵ μ = 1 3 α ϵ μ α α .
In addition to computing ω μ , a μ , and Z i α β * , the response function capabilities of ABINIT can also produce ϵ , which is the directionally averaged electronic contribution to the dielectric tensor. This is sometimes referred to as the high-field limit response, where the application of a high frequency alternating electric field prevents ionic motion and any response is purely electronic in nature. This separation in behavior, ionic vs. electric, means that the total static dielectric response can therefore be written as:
ϵ = ϵ + μ ϵ μ .
Here, we focus on characterizing the IR-active modes in the region of 20–110 cm 1 , in line with the collected THz data, with meaningful contributions to the dielectric response. For all comparative IR active mode plots, the DFT-computed spectral features are generated using a Gaussian distribution with a full width at half max of 1.5 cm 1 . The relative DFT-computed peak intensities, plotted as a phonon density of states (DOS), are determined by their specific contributions to ϵ in the range of 20–110 cm 1 , which is specified for the spectroscopic analysis presented here.

4. Conclusions

The comparison of DFT to DFT-D for the organic compounds BPE, BPEth, and SA shows that the inclusion of dispersion using the semiempirical Grimme-D2 dispersion correction will yield lattice parameters and vibrational modes whose values are closer to those measured in experiment. We find that the agreement between DFT-D computed vibrational modes and THz spectra was best for BPE and worst for SA. As pointed out in the final subsection of the Results and Discussion, this may be caused by using an s 6 global scaling factor that was too large, 0.75 vs. ≈0.50, or could also be because of the increased amount of hydrogen bonding in the SA crystal structure. With regards to the co-crystals 2(SA)·BPEth and 2(SA)·BPE forms I and II, DFT-D did pick out differentiable trends between the systems, but the lattice constant disagreement (for 2(SA)·BPEth and 2(SA)·BPE-II) between theory and experiment was still too large to be able to assign specific vibrational features to the THz spectra. Moreover, the experimental spectra showed that atomic packing and overall crystal lattice structure dominated the THz region moreso than changes in ligand identity as BPE vs. BPEth, while DFT-D showed that both atomic packing and ligand identity led to differentiable regions where each tended to dominate. One way to obtain more information about the DFT-D computed vibrational modes, and to match more closely to experiments, would be to use what is presented here as input for molecular dynamics calculations that take into account effects such as temperature or by quasi-harmonic DFT approximations [69] that map out the effects of temperature and volume changes on solid-state properties related to low frequency vibrations [70,71,72].
Our study determines the optimal ranges of s 6 for the set of materials under investigation and we are able to discern differences in s 6 that may be correlated to crystal structure/property relations and be used as a guide for further investigations. Organic heterocycles like BPE and BPEth with π π stacking interactions can be optimized with s 6 ≈ 0.3–0.5, and organics that demonstrate significant H-bonding, like SA, can be optimized with s 6 ≈ 0.4–0.9. Comparing the s 6 of the combination, these two different types of bonding schemes (to create co-crystals) also yield insights; the lattice parameters of 2(SA)·BPEth co-crystal are closest to experiments when s 6 approaches the values used for organic materials with π π stacking and not those with significant hydrogen-bonding. This suggests that, in complex materials, containing multiple functional groups or packing behavior, at least two or three different s 6 should be used for preliminary structural optimizations. Determining which range of s 6 may be best will allow for a better understanding of the relative forces that stabilize a complex structure.
Our analysis of the dielectric response of BPE, BPEth, and SA shows that functional group identity, degree of hydrogen bonding, and crystal packing/arrangement will influence the dielectric constant, ϵ . This work could be extended to analyze and compare the effects of fluoride, cyano, aldehyde, and amide substituents on aromatic heterocycles similar to what is presented here, or on different types of heterocycles (pyrazine, imidazole, furan, thiophene) and 3D bonding architectures (fullerenes, substituted adamantanes) on dielectric behavior, to be used as a design parameter for co-crystal APIs.
Here, we establish a baseline of DFT computed dielectric response and THz spectroscopy for the organic co-crystals 2(SA)·BPEth and 2(SA)·BPE forms I and II, as well as their constituent components, using dispersion corrected planewave DFT methods that employ the GBRV set of ultrasoft pseudopotentials. Further development of the work presented here could include using our Grimme-D2 results as a starting point for calculations that take into account molecular environments, such as the parameter free approach developed by Tkatchenko and Scheffler [73] or the Grimme-D3 methodology, which employs a three-body dispersion term [74] and fractional coordination numbers. Recently developed dispersion corrected hybrid functionals [75,76] may also improve the agreement between DFT modeling and THz measurements beyond what is presented here using a combination of DFT-GGA and Grimme-D2 dispersion correction.

Author Contributions

The authors contributions are as follows: conceptualization, J.W.B., M.E.R., S.M.O., L.R.M., M.A.A., and S.E.M.; methodology, J.W.B., M.E.R., S.M.O., L.R.M., M.A.A., and S.E.M.; formal analysis, J.W.B. and M.E.R.; investigation, J.W.B., M.E.R. and S.M.O.; data curation, J.W.B., M.E.R., and S.M.O.; writing—original draft preparation, J.W.B., M.E.R., M.A.A., and S.E.M.; writing—review and editing, J.W.B., M.E.R., S.M.O., L.R.M., M.A.A., and S.E.M.; visualization, J.W.B. and M.E.R.; supervision, L.R.M., M.A.A., and S.E.M.; project administration, L.R.M., M.A.A., and S.E.M.; funding acquisition, L.R.M and S.E.M.

Funding

This research was supported in part through computational resources provided by The University of Iowa, Iowa City, Iowa and the National Science Foundation grant CHE-0840494. This work used the Extreme Science and Engineering Discovery Environment (XSEDE [77]), which is supported by the National Science Foundation Grant No. ACI-1548562 through allocation ID TG-GEO160006. L.R.M. also acknowledges the National Science Foundation (DMR-1708673) for funding.

Acknowledgments

This work was supported by the University of Iowa, College of Liberal Arts and Sciences.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

References

  1. Steed, J.W. The Role of Co-Crystals in Pharmaceutical Design. Trends Pharm. Sci. 2013, 34, 185–193. [Google Scholar] [CrossRef] [PubMed]
  2. Chieng, N.; Rades, T.; Aaltonen, J. An overview of recent studies on the analysis of pharmaceutical polymorphs. J. Pharm. Biomed. Anal. 2011, 55, 618–644. [Google Scholar] [CrossRef] [PubMed]
  3. Vishweshwar, P.; McMahon, J.A.; Bis, J.A.; Zaworotko, M.J. Pharmaceutical Co-Crystals. J. Pharm. Sci. 2006, 95, 499–516. [Google Scholar] [CrossRef] [PubMed]
  4. 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. [Google Scholar] [CrossRef] [PubMed]
  5. Wallace, V.P.; Taday, P.F.; Fitzgerald, A.J.; Woodward, R.M.; Cluff, J.; Pye, R.J.; Arone, D.D. Terahertz Pulsed Imaging and Spectroscopy for Biomedical and Pharmaceutical Applications. Faraday Discuss. 2004, 126, 255–263. [Google Scholar] [CrossRef] [PubMed]
  6. Berry, D.J.; Steed, J.W. Pharmaceutical Cocrystals, Salts and Multicomponent Systems; Intermolecular Interactions and Property Based Design. Adv. Drug Deliv. Rev. 2017, 117, 3–24. [Google Scholar] [CrossRef]
  7. Braga, D.; Maini, L.; Grepioni, F. Mechanochemical preparation of co-crystals. Chem. Soc. Rev. 2013, 42, 7638–7648. [Google Scholar] [CrossRef] [PubMed]
  8. Urbanus, J.; Roelands, C.P.M.; Verdoes, D.; Jansens, P.J.; ter Horst, J.H. Co-Crystallization as a Separation Technology: Controlling Product Concentrations by Co-Crystals. Cryst. Growth Des. 2010, 10, 1171–1179. [Google Scholar] [CrossRef]
  9. McIntosh, A.I.; Yang, B.; Goldup, S.M.; Watkinson, M.; Donnan, R.S. Terahertz Spectroscopy: A Powerful New Tool for the Chemical Sciences? Chem. Soc. Rev. 2012, 41, 2072–2082. [Google Scholar] [CrossRef]
  10. Ruggiero, M.T.; Sutton, J.J.; Fraser-Miller, S.J.; Zaczek, A.J.; Korter, T.M.; Gordon, K.C.; Zeitler, J.A. Revisiting the Thermodynamic Stability of Indomethacin Polymorphs with Low-Frequency Vibrational Spectroscopy and Quantum Mechanical Simulations. Cryst. Growth Des. 2018, 18, 6513–6520. [Google Scholar] [CrossRef]
  11. Taday, P.F.; Bradley, I.V.; Arnone, D.D.; Pepper, M. Using Terahertz Pulse Spectroscopy to Study the Crystalline Structure of a Drug: A Case Study of the Polymorphs of Ranitidine Hydrochloride. J. Pharm. Sci. 2003, 92, 831–838. [Google Scholar] [CrossRef] [PubMed]
  12. Larkin, P.J.; Wasylyk, J.; Raglione, M. Application of Low- and Mid-Frequency Raman Spectroscopy to Characterize the Amorphous-Crystalline Transformation of Indomethacin. Appl. Spectrosc. 2015, 69, 1217–1228. [Google Scholar] [CrossRef] [PubMed]
  13. Zeitler, J.A.; Newnham, D.A.; Taday, P.F.; Threlfall, T.L.; Lancaster, R.W.; Berg, R.W.; Strahan, C.J.; Pepper, M.; Gordon, K.C.; Rades, T. Characterization of Temperature-Induced Phase Transitions in Five Polymorphic Forma of Sulfathiazole by Terahertz Pulsed Spectroscopy and Differential Scanning Calorimetry. J. Pharm. Sci. 2006, 95, 2486–2498. [Google Scholar] [CrossRef] [PubMed]
  14. Inoue, M.; Hisada, H.; Koide, T.; Carriere, J.; Heyler, R.; Fukami, T. In Situ Monitoring of Crystalline Transformation of Carbamazepine Using Probe-Type Low Frequency Raman Spectroscopy. Org. Process Res. Dev. 2017, 21, 262–265. [Google Scholar] [CrossRef]
  15. Zeitler, J.A.; Taday, P.F.; Gordon, K.C.; Pepper, M.; Rades, T. Solid-State Transition Mechanism in Carbamazepine Polymorphs by Time-Resolved Terahertz Spectroscopy. Chem. Phys. Chem. 2007, 8, 1924–1927. [Google Scholar] [CrossRef] [PubMed]
  16. Neu, J.; Nemes, C.T.; Regan, K.P.; Williams, M.R.C.; Schmuttenmaer, C.A. Exploring the solid state phase transition in DL-norvaline with terahertz spectroscopy. Phy. Chem. Chem. Phys. 2018, 20, 276–283. [Google Scholar] [CrossRef]
  17. Nyugen, K.L.; Friscic, T.; Day, G.M.; Gladden, L.F.; Jones, W. Terahertz time-domain spectroscopy and the quantitative monitoring of mechanochemical cocrystal formation. Nat. Mater. 2007, 6, 206–209. [Google Scholar]
  18. Delaney, S.P.; Korter, T.M. Terahertz Spectroscopy and Computational Investigation of the Flufenamic Acid/Nicotinamide Cocrystal. J. Phys. Chem. A 2015, 119, 3269–3276. [Google Scholar] [CrossRef] [PubMed]
  19. King, M.D.; Ouellette, W.; Korter, T.M. Noncovalent Interactions in Paired DNA Nucleobases Investigated by Terahertz Spectroscopy and Solid-State Density Functional Theory. J. Phys. Chem. A 2011, 115, 9467–9478. [Google Scholar] [CrossRef] [PubMed]
  20. Druzbicki, K.; Mielcarek, J.; Kiwilsza, A.; Troupet, L.; Collet, E.; Pajzderska, A.; Wasicki, J. Computationally Assisted (Solid State Density Functional Theory) Structural (X-ray) and Vibrational Spectroscopy (FT-IR, FT-RS, TDs-THz) Characterization of the Cardiovascular Drug Lacidipine. Cryst. Growth Des. 2015, 15, 2817–2830. [Google Scholar] [CrossRef]
  21. Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, B864–B871. [Google Scholar] [CrossRef]
  22. Kohn, W.; Sham, L.J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, A1133. [Google Scholar] [CrossRef]
  23. Stillinger, F.H.; Weber, T.A. Computer simulation of local order in condensed phases of silicon. Phys. Rev. B 1985, 31, 5262–5271. [Google Scholar] [CrossRef]
  24. Bernard, J.E.; Zunger, A. Electronic structure of ZnS, ZnSe, ZnTe, and their pseudobinary alloys. Phys. Rev. B 1987, 36, 3199–3228. [Google Scholar] [CrossRef]
  25. Rappe, A.M.; Rabe, K.M.; Kaxiras, E.; Joannopoulos, J.D. Optimized Pseudopotentials. Phys. Rev. B Rapid Commun. 1990, 41, 1227–1230. [Google Scholar] [CrossRef]
  26. Ramer, N.J.; Rappe, A.M. Designed Nonlocal Pseudopotentials for Enhanced Transferability. Phys. Rev. B 1999, 59, 12471–12478. [Google Scholar] [CrossRef]
  27. Vanderbilt, D. Soft Self-Consistent Pseudopotentials in a Generalized Eigenvalue Formalism. Phys. Rev. B Rapid Commun. 1990, 41, 7892–7895. [Google Scholar] [CrossRef]
  28. Hamann, D.R. Optimized norm-conserving Vanderbilt pseudopotentials. Phys. Rev. B 2013, 88, 085117. [Google Scholar] [CrossRef] [Green Version]
  29. Garrity, K.F.; Bennett, J.W.; Rabe, K.M.; Vanderbilt, D. Pseudopotentials for high-throughput DFT calculations. Comput. Mater. Sci. 2014, 81, 446–452. [Google Scholar] [CrossRef] [Green Version]
  30. Lejaeghere, K.; Bihlmayer, G.; Bjorkman, T.; Blaha, P.; Blugel, S.; Blum, V.; Caliste, D.; Castelli, I.E.; Clark, S.J.; Corso, A.D.; et al. Reproducibility in density functional theory calculations of solids. Science 2016, 351, aad3000. [Google Scholar] [CrossRef]
  31. Perdew, J.P.; Chevary, J.A.; Vosko, S.H.; Jackson, K.A.; Pederson, M.R.; Singh, D.J.; Fiolhais, C. Atoms, molecules, solids and surfaces: Applications of the generalized gradient approxiamtion for exchange and correlations. Phys. Rev. B 1992, 46, 6671–6687. [Google Scholar] [CrossRef]
  32. Perdew, J.P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865–3868. [Google Scholar] [CrossRef] [PubMed]
  33. Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158–6170. [Google Scholar] [CrossRef]
  34. Gonze, X.; Amadon, B.; Anglade, P.; Beuken, J.M.; Bottin, F.; Boulanger, P.; Bruneval, F.; Caliste, D.; Caracas, R.; Cote, M.; et al. ABINIT: First-Principles Approach to Material and Nanosystem Properties. Comp. Phys. Commun. 2009, 180, 2582–2615. [Google Scholar] [CrossRef]
  35. Giannozzi, P.; Baroni, S.; Bonini, N.; Calandra, M.; Car, R.; Cavazzoni, C.; Ceresoli, D.; Chiarotti, G.L.; Cococcioni, M.; Dabo, I.; et al. QUANTUM-ESPRESSO: A Modular and Open-Source Software Project for Quantum Simulation of Materials. J. Phys. Condens. Matter 2009, 21, 395502. [Google Scholar] [CrossRef] [PubMed]
  36. Giustino, F. Materials Modelling Using Density Functional Theory: Properties and Preditions; Oxford University Press: Oxford, UK, 2014. [Google Scholar]
  37. Roy, A.; Bennett, J.W.; Rabe, K.M.; Vanderbilt, D. Half-Heusler Semiconductors as Piezoelectrics. Phys. Rev. Lett. 2012, 109, 037601. [Google Scholar] [CrossRef] [PubMed]
  38. Bennett, J.W.; Garrity, K.F.; Rabe, K.M.; Vanderbilt, D. Orthorhombic ABC Semiconductors as Antiferroelectrics. Phys. Rev. Lett. 2013, 110, 017603. [Google Scholar] [CrossRef]
  39. Lee, J.H.; Rabe, K.M. Large spin-phonon coupling and magnetically induced phonon anisotropy in SrMO3 perovskites (M=V, Cr, Mn, Fe, Co). Phys. Rev. B 2011, 84, 104440. [Google Scholar] [CrossRef]
  40. Lee, J.H.; Rabe, K.M. Coupled Magnetic-Ferroelectric Metal-Insulator Transition in Epitaxially Strianed SrCoO3 from First-Principles. Phys. Rev. Lett. 2011, 107, 067601. [Google Scholar] [CrossRef]
  41. Yu, L.; Kokenyesi, R.S.; Kaszler, D.A.; Zunger, A. Inverse Design of High Absorption Thin-Film Photovoltaic Materials. Adv. Energy Mater. 2013, 3, 43–48. [Google Scholar] [CrossRef]
  42. Dion, M.; Rydberg, H.; Schroder, E.; Langreth, D.C.; Lundquist, B.I. Van der Waals Density Functional for General Geometries. Phys. Rev. Lett. 2004, 92, 246401. [Google Scholar] [CrossRef] [PubMed]
  43. Grimme, S. Semiempirical GGA-type density functional constructed with a long-range dispersion correction. J. Comput. Chem. 2006, 27, 1787–1799. [Google Scholar] [CrossRef] [PubMed]
  44. Schwabe, T.; Grimme, S. Double-hybrid density functionals with long-range dispersion corrections: Higher accuracy and extended applicability. Phys. Chem. Chem. Phys. 2007, 9, 3397–3406. [Google Scholar] [CrossRef] [PubMed]
  45. Chai, J.D.; Head-Gordon, M. Long-range corrected hybrid density functionals with damped atom-atom dispersion corrections. Phys. Chem. Chem. Phys. 2008, 10, 6615–6620. [Google Scholar] [CrossRef] [PubMed]
  46. Grimme, S.; Antony, J.; Ehrlich, S.; Krieg, H. A consistent and accurate ab initio parameterization of density functional dispeersion correction (DFT-D) for the 94 elements H-Pu. J. Chem. Phys. 2010, 132, 154104. [Google Scholar] [CrossRef] [PubMed]
  47. Berland, K.; Cooper, V.R.; Lee, K.; Schroder, E.; Thonhauser, T.; Hyldgaard, P.; Lundqvist, B.I. van der Waals forces in density functional theory: a review of the vdW-DF method. Rep. Prog. Phys. 2015, 78, 066501. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  48. Thanthiriwatte, K.S.; Hohenstein, E.G.; Burns, L.A.; Sherrill, C.D. Assessment of the Performance of DFT and DFT-D methods for describing distance dependence of hydrogen-bonded interactions. J. Chem. Theory Comput. 2011, 7, 88–96. [Google Scholar] [CrossRef] [PubMed]
  49. Peverati, R.; Baldridge, K.K. Implementation and Performance of DFT-D with Respect to Basis Set and Functional for Study of Dispersion Interactions in Nanoscale Aromatic Hydrocarbons. J. Chem. Theory Comput. 2008, 4, 2030–2048. [Google Scholar] [CrossRef] [PubMed]
  50. Caldeweyher, E.; Bannwarth, C.; Grimme, S. Extension of the D3 Dispersion Coefficient. J. Chem. Phys. 2017, 147, 034112. [Google Scholar] [CrossRef] [PubMed]
  51. Grimme, S. Density functional theory with London dispersion corrections. Wiley Interdiscip. Rev. 2011, 1, 211–228. [Google Scholar] [CrossRef]
  52. Vansant, J.; Smets, G.; Declercq, J.P.; Germain, G.; Meerssche, M.V. Azastilbenes. 1. Synthesis, characterization, and structure. J. Org. Chem. 1980, 45, 1557–1565. [Google Scholar] [CrossRef]
  53. Ide, S.; Karacan, N.; Tufan, Y. 1,2-Bis(4-pyridyl)ethane. Acta Cryst. C 1995, 51, 2304–2305. [Google Scholar] [CrossRef] [Green Version]
  54. Sundaralingham, M.; Jensen, L.H. Refinement of the structure of salicylic acid. Acta Cryst. 1965, 18, 1053–1058. [Google Scholar] [CrossRef] [Green Version]
  55. 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. [Google Scholar] [CrossRef] [PubMed]
  56. Raglione, M.E.; Zhang, T.; Arnold, M.A. Stability of Polytetrafluoroethylene, Polyethylene, and Polystyrene Pellets as a Matrix for Analytical Terahertz Spectroscopy. Int. J. Exp. Spectrosc. Tech. 2018, 3. in press. [Google Scholar]
  57. Rupasinghe, T.P.; Hutchins, K.M.; Bandaranayake, B.S.; Ghorai, S.; Karunatilake, C.; Bucar, D.K.; Swenson, D.C.; Arnold, M.A.; MacGillivray, L.R.; Tivanski, A.V. Mechanical Properties of a Series of Macro- and Nanodimensional Organic Cocrystals Correlate with Atomic Polarizability. J. Am. Chem. Soc. 2015, 137, 12768–12771. [Google Scholar] [CrossRef] [PubMed]
  58. Naftaly, M.; Miles, R.E. Terahertz time-domain spectroscopy of silicate glasses and the relationship to material properties. J. Appl. Phys. 2007, 102, 043517. [Google Scholar] [CrossRef]
  59. Looyenga, H. Dielectric constants of homogeneous mixture. Mol. Phys. 1965, 9, 501–511. [Google Scholar] [CrossRef]
  60. Sun, J.; Lucyszyn, S. Extracting Complex Dielectric Properties from Reflection-Transmission Mode Spectroscopy. IEEE Access 2018, 6, 1–20. [Google Scholar] [CrossRef]
  61. Monkhorst, H.J.; Pack, J.D. Special Points for Brillouin-Zone Integrations. Phys. Rev. B 1976, 13, 5188–5192. [Google Scholar] [CrossRef]
  62. Cockayne, E.; Burton, B.P. Phonons and Static Dielectric Constant in CaTiO3 from First Principles. Phys. Rev. B 2000, 62, 3735–3743. [Google Scholar] [CrossRef]
  63. Cockayne, E. First-principles calculations of the dielectric properties of perovskite-type materials. J. Eur. Ceram. 2003, 23, 2375–2379. [Google Scholar] [CrossRef]
  64. Bennett, J.W.; Grinberg, I.; Rappe, A.M. Effect of symmetry lowering on the dielectric response of BaZrO3. Phys. Rev. B 2006, 73, 180102R. [Google Scholar] [CrossRef]
  65. Bennett, J.W.; Grinberg, I.; Rappe, A.M. Nonmonotonic Composition Dependence of the Dielectric Response of Ba1−xCaxZrO3. Chem. Mater. 2008, 20, 5134–5138. [Google Scholar] [CrossRef]
  66. Bennett, J.W.; Grinberg, I.; Rappe, A.M. Effect of substituting S for O: The sulfide perovskite BaZrS3 investigated with density functional theory. Phys. Rev. B 2009, 79, 235115. [Google Scholar] [CrossRef]
  67. Gonze, X.; Lee, C. Dynamical matrices, Born effective charges, dielectric permittivity tensors, and interatomic force constants from density-functional perturbation theory. Phys. Rev. B 1997, 55, 10355–10368. [Google Scholar] [CrossRef]
  68. Ghosez, P.; Michenaud, J.P.; Gonze, X. Dynamical atomic charges: The case of ABO3 compounds. Phys. Rev. B 1998, 58, 6224–6240. [Google Scholar] [CrossRef]
  69. Zhang, F.; Wang, H.W.; Tominaga, K.; Hayashi, M. Characteristics of Low-Frequency Molecular Phonon Modes Studied by THz Spectroscopy and Solid-State ab initio Theory: Polymorphs I and III of Diflunisal. J. Phys. Chem. B 2016, 120, 1698–1710. [Google Scholar] [CrossRef]
  70. Ruggiero, M.T.; Zeitler, J.A.; Erba, A. Intermolecular anharmonicity in molecular crystals: Interplay between experimental low-frequency dynamics and quantum quasi-harmonic simulations of solid purine. Chem. Commun. 2017, 53, 3781–3784. [Google Scholar] [CrossRef]
  71. Brandenburg, J.G.; Potticary, J.; Sparkes, H.A.; Price, S.L.; Hall, S.R. Thermal Expansion of Carbamazepine: Systematic Crystallographic Measurements Challenge Quantum Chemical Calculations. J. Phys. Chem. Lett. 2017, 8, 4319–4324. [Google Scholar] [CrossRef]
  72. Ruggiero, M.T.; Zhang, W.; Bond, A.D.; Mittleman, D.M.; Zeitler, J.A. Uncovering the Connection Between Low-Frequency Dynamics and Phase Transformation Phenomena in Molecular Solids. Phys. Rev. Lett. 2018, 120, 196002. [Google Scholar] [CrossRef] [PubMed]
  73. Tkatchenko, A. Accurate Molecular Van der Waals Interactions from Ground State Electron Density and Free-Atom Reference Data. Phys. Rev. Lett. 2009, 102, 073005. [Google Scholar] [CrossRef] [PubMed]
  74. Moellmann, J.; Grimme, S. DFT-D3 Study of Some Molecular Crystals. J. Phys. Chem. C 2014, 118, 7615–7621. [Google Scholar] [CrossRef]
  75. Wang, C.W.; Hui, K.; Chai, J.D. Short- and long-range corrected hybrid density functionals with the D3 dispersion corrections. J. Chem. Phys. 2016, 145, 204101. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  76. Goerigk, L. Treating London-Dispersion Effects with the Latest Minnesota Density Functionals: Problems and Possible Solutions. J. Phys. Chem. Lett. 2015, 6, 3891–3896. [Google Scholar] [CrossRef] [PubMed]
  77. Towns, J.; Cockerill, T.; Dahan, M.; Foster, I.; Gaither, K.; Grimshaw, A.; Hazlewood, V.; Lanthrop, S.; Lifka, D.; Peterson, G.D.; et al. XSEDE:Accelerating Scientific Discovery. Comp. Sci. Eng. 2014, 16, 62–74. [Google Scholar] [CrossRef]
Figure 1. Shown here are the co-crystal components (a) BPE, (b) BPEth, and (c) SA. Atoms of carbon, hydrogen, nitrogen, and oxygen are depicted as black, light pink, light blue, and red spheres, respectively. Dashed lines represent the crystallographic axes of the primitive lattice unit cell for each of the components.
Figure 1. Shown here are the co-crystal components (a) BPE, (b) BPEth, and (c) SA. Atoms of carbon, hydrogen, nitrogen, and oxygen are depicted as black, light pink, light blue, and red spheres, respectively. Dashed lines represent the crystallographic axes of the primitive lattice unit cell for each of the components.
Molecules 24 00959 g001
Figure 2. Comparisons of the computed phonon modes (black solid curves) to experimentally obtained THz spectra (solid blue curves) of BPE using (a) DFT and (b) DFT-D. Dashed black vertical lines are the normalized phonon mode intensity used to obtain the phonon density of states (DOS). All measurements are given in arbitrary units.
Figure 2. Comparisons of the computed phonon modes (black solid curves) to experimentally obtained THz spectra (solid blue curves) of BPE using (a) DFT and (b) DFT-D. Dashed black vertical lines are the normalized phonon mode intensity used to obtain the phonon density of states (DOS). All measurements are given in arbitrary units.
Molecules 24 00959 g002
Figure 3. Depicted here are the DFT-D computed IR active vibrational modes of BPE that occur at (a) 105.6 cm−1 (along the x-direction) and (b) 106.4 cm−1 (along the y-direction. Color scheme is the same as before, but red arrows indicate relative atomic displacements. The displacements of H atoms are not shown for clarity.
Figure 3. Depicted here are the DFT-D computed IR active vibrational modes of BPE that occur at (a) 105.6 cm−1 (along the x-direction) and (b) 106.4 cm−1 (along the y-direction. Color scheme is the same as before, but red arrows indicate relative atomic displacements. The displacements of H atoms are not shown for clarity.
Molecules 24 00959 g003
Figure 4. Same description as in Figure 2, but for 1,2-bis(4-pyridyl)ethane (BPEth).
Figure 4. Same description as in Figure 2, but for 1,2-bis(4-pyridyl)ethane (BPEth).
Molecules 24 00959 g004
Figure 5. Depicted here are the DFT-D computed IR active vibrational modes of BPEth that occur at (a) 61.2 cm−1 (along the y-direction) and (b) 70.9 cm−1 (along the x-direction).
Figure 5. Depicted here are the DFT-D computed IR active vibrational modes of BPEth that occur at (a) 61.2 cm−1 (along the y-direction) and (b) 70.9 cm−1 (along the x-direction).
Molecules 24 00959 g005
Figure 6. Same description as in Figure 2, but for salicylic acid (SA).
Figure 6. Same description as in Figure 2, but for salicylic acid (SA).
Molecules 24 00959 g006
Figure 7. Shown here are the DFT-D computed IR active modes of SA that occur at (a) 37.9 cm−1 (along y); (b) 76.1 −1 (along xz); and (c) 92.3 −1 (along y).
Figure 7. Shown here are the DFT-D computed IR active modes of SA that occur at (a) 37.9 cm−1 (along y); (b) 76.1 −1 (along xz); and (c) 92.3 −1 (along y).
Molecules 24 00959 g007
Figure 8. Shown here are the co-crystals (a) 2(SA)·BPE-I; (b) 2(SA)·BPE-II; and (c) 2(SA)·BPEth. Color scheme is as before in Figure 1.
Figure 8. Shown here are the co-crystals (a) 2(SA)·BPE-I; (b) 2(SA)·BPE-II; and (c) 2(SA)·BPEth. Color scheme is as before in Figure 1.
Molecules 24 00959 g008
Figure 9. Shown here are comparisons of the computed phonon modes (black solid curves) to experimentally obtained THz spectra (solid blue curves) of (a) 2(SA)·BPE-I; (b) 2(SA)·BPE-II; and (c) 2(SA)·BPEth using DFT-D. Dashed black vertical lines are the normalized phonon mode intensity used to obtain the phonon density of states (DOS). All measurements are given in arbitrary units.
Figure 9. Shown here are comparisons of the computed phonon modes (black solid curves) to experimentally obtained THz spectra (solid blue curves) of (a) 2(SA)·BPE-I; (b) 2(SA)·BPE-II; and (c) 2(SA)·BPEth using DFT-D. Dashed black vertical lines are the normalized phonon mode intensity used to obtain the phonon density of states (DOS). All measurements are given in arbitrary units.
Molecules 24 00959 g009
Figure 10. Comparisons of the Experimental THz spectra (top) and DFT-D computed phonon modes (bottom) of (left) 2(SA)·BPE-I (solid yellow line) and 2(SA)·BPE-II (dash-dotted purple line), and (right) 2(SA)·BPE-II (dash-dotted purple line) and 2(SA)·BPEth (solid orange line).
Figure 10. Comparisons of the Experimental THz spectra (top) and DFT-D computed phonon modes (bottom) of (left) 2(SA)·BPE-I (solid yellow line) and 2(SA)·BPE-II (dash-dotted purple line), and (right) 2(SA)·BPE-II (dash-dotted purple line) and 2(SA)·BPEth (solid orange line).
Molecules 24 00959 g010
Figure 11. Shown here are the DFT-D computed IR active modes of 2(SA)·BPEth that occur at (a) 77.7 cm−1 (along xz); (b) 79.9 cm−1 (along y); (c) 89.8 cm−1 (along xz); (d) 91.3 cm−1 (along y); and (e) 93.2 cm−1 (along xz).
Figure 11. Shown here are the DFT-D computed IR active modes of 2(SA)·BPEth that occur at (a) 77.7 cm−1 (along xz); (b) 79.9 cm−1 (along y); (c) 89.8 cm−1 (along xz); (d) 91.3 cm−1 (along y); and (e) 93.2 cm−1 (along xz).
Molecules 24 00959 g011
Figure 12. Plotted here are the DFT-D computed IR active modes of 2(SA)·BPE-I that occur at (a) 28.9 cm−1 (along xz); (b) 33.9 cm−1 (along xz); (c) 34.9 cm−1 (along y); (d) 69.8 cm−1 (along y); and (e) 70.3 cm−1 (along xz).
Figure 12. Plotted here are the DFT-D computed IR active modes of 2(SA)·BPE-I that occur at (a) 28.9 cm−1 (along xz); (b) 33.9 cm−1 (along xz); (c) 34.9 cm−1 (along y); (d) 69.8 cm−1 (along y); and (e) 70.3 cm−1 (along xz).
Molecules 24 00959 g012
Figure 13. Depicted here are the DFT-D computed IR active modes of 2(SA)·BPE-II that occur at (a) 42.3 cm−1 (along xz); (b) 68.5 cm−1 (along z); (c) 69.1 cm−1 (along y); (d) 92.3 cm−1 (along xz); and (e) 102.9 cm−1 (along xz).
Figure 13. Depicted here are the DFT-D computed IR active modes of 2(SA)·BPE-II that occur at (a) 42.3 cm−1 (along xz); (b) 68.5 cm−1 (along z); (c) 69.1 cm−1 (along y); (d) 92.3 cm−1 (along xz); and (e) 102.9 cm−1 (along xz).
Molecules 24 00959 g013
Figure 14. For co-crystal components (a) BPE; (b) BPEth; and (c) SA, we vary s6 from 0 to 1 to determine how DFT-D computed lattice parameters a, b, c, and β will be affected when compared to experimentally determined lattice parameters. The same tests for co-crystal 2(SA)·BPEth yield a similar range of s6 in which lattice constant underestimation may be minimized. Brown dashed lines denote ±4% error with respect to experimentally determined lattice parameters.
Figure 14. For co-crystal components (a) BPE; (b) BPEth; and (c) SA, we vary s6 from 0 to 1 to determine how DFT-D computed lattice parameters a, b, c, and β will be affected when compared to experimentally determined lattice parameters. The same tests for co-crystal 2(SA)·BPEth yield a similar range of s6 in which lattice constant underestimation may be minimized. Brown dashed lines denote ±4% error with respect to experimentally determined lattice parameters.
Molecules 24 00959 g014
Table 1. Lattice parameters of co-crystal components BPE, BPEth, and SA, using DFT (top) and DFT-D (bottom). All lattice constants are reported in units of Å, and both the experimental value and % deviation from experimentally determined values are given in parentheses.
Table 1. Lattice parameters of co-crystal components BPE, BPEth, and SA, using DFT (top) and DFT-D (bottom). All lattice constants are reported in units of Å, and both the experimental value and % deviation from experimentally determined values are given in parentheses.
a DFT b DFT c DFT β DFT
(Å)(Å)(Å)
BPE11.00 (7.82, +40.66)10.38 (10.56, −1.70)6.89 (5.77, +19.41)116.89 (92.68, +26.12)
BPEth5.69 (5.56, +2.33)9.64 (8.16, +18.14)11.33 (11.35, −0.18)99.26 (100.73, −1.45)
SA5.20 (4.89, +6.38)11.02 (11.20, −1.64)14.83 (11.24, +31.94)75.39 (92.49, −18.49)
BPE7.30 (7.82, −6.66)10.51 (10.56, −0.50)5.56 (5.77, −2.87)90.39 (92.68, −2.47)
BPEth5.31 (5.56, −4.53)7.80 (8.16, −4.38)11.11 (11.35, −2.16)98.73 (100.73, −1.99)
SA4.85 (4.89, −0.78)11.02 (11.20, −1.64)11.27 (11.24, +0.23)93.54 (92.49, +1.14)
Table 2. Comparison of experimentally determined dielectric constant, ϵ , to the isotropically averaged values computed using DFT and DFT-D for the co-crystal components BPE, BPEth, and SA. Percent deviation from experiment is shown in parenthesis for DFT and DFT-D.
Table 2. Comparison of experimentally determined dielectric constant, ϵ , to the isotropically averaged values computed using DFT and DFT-D for the co-crystal components BPE, BPEth, and SA. Percent deviation from experiment is shown in parenthesis for DFT and DFT-D.
ExperimentDFTDFT-D
BPE3.43 ± 0.082.68 (−21.87)3.92 (+14.29)
BPEth3.13 ± 0.152.68 (−14.38)3.25 (+3.83)
SA3.09 ± 0.043.70 (+19.74)2.79 (−9.71)
Table 3. Mode-by-mode analysis of the directional components of the IR active response of BPE for DFT (top) and DFT-D (bottom) for the THz frequency range of 20–110 cm 1 . ω are given in units of cm 1 and all ϵ are unitless. Asterisks (*) are next to the frequencies with high contributions to ϵ μ , and are shown in Figure 3. The two final rows in each of the DFT and DFT-D calculations are the total ionic portion of the dielectric response per direction x, y, z per mode μ (denoted as ϵ μ ), and the directionally decomposed electronic contribution of the dielectric response ϵ .
Table 3. Mode-by-mode analysis of the directional components of the IR active response of BPE for DFT (top) and DFT-D (bottom) for the THz frequency range of 20–110 cm 1 . ω are given in units of cm 1 and all ϵ are unitless. Asterisks (*) are next to the frequencies with high contributions to ϵ μ , and are shown in Figure 3. The two final rows in each of the DFT and DFT-D calculations are the total ionic portion of the dielectric response per direction x, y, z per mode μ (denoted as ϵ μ ), and the directionally decomposed electronic contribution of the dielectric response ϵ .
Method ω ϵ x ϵ y ϵ z
DFT31.1-0.29-
47.2-0.01-
47.5-0.03-
55.20.02-0.03
65.8-0.07-
70.50.01-0.03
81.90.02-0.00
92.0-0.06-
97.00.02-0.00
ϵ μ 0.120.550.15
ϵ 1.972.862.40
DFT-D67.6-0.02-
69.50.01-0.00
105.6 *0.15-0.00
106.4 *-0.15-
ϵ μ 0.440.300.14
ϵ 2.634.633.62
Table 4. Same description as Table 3, but for BPEth. Asterisks (*) are next to the vibrational modes shown in Figure 5.
Table 4. Same description as Table 3, but for BPEth. Asterisks (*) are next to the vibrational modes shown in Figure 5.
Method ω xyz
DFT31.8-0.11-
32.5-0.10-
41.60.10-0.16
54.8-0.03-
60.5-0.06-
68.80.10-0.01
70.8-0.11-
ϵ μ 0.260.480.22
ϵ 2.202.182.69
DFT-D45.4-0.06-
53.90.03-0.00
61.2 *-0.08-
70.9 *0.13-0.04
93.8-0.01-
108.90.00-0.01
ϵ μ 0.240.430.13
ϵ 2.972.553.43
Table 5. Same description as Table 3, but for SA. Asterisks (*) denote the modes depicted in Figure 7.
Table 5. Same description as Table 3, but for SA. Asterisks (*) denote the modes depicted in Figure 7.
Method ω xyz
DFT30.8-0.04-
37.00.02-0.00
42.2-0.02-
42.30.02-0.00
58.90.03-0.00
63.0-0.22-
71.30.03-0.00
79.2-0.04-
79.50.09-0.00
ϵ μ 0.791.010.44
ϵ 3.003.112.73
DFT-D37.9 *-0.29-
49.70.05-0.00
66.5-0.04-
76.1 *0.07-0.03
78.4-0.02-
92.3 *-0.09-
100.20.03-0.01
ϵ μ 0.550.690.28
ϵ 2.192.322.34
Table 6. Lattice parameters of co-crystals using DFT (top) and DFT-D (bottom). All lattice constants are reported in units of Åand % deviation from experimentally determined values are given in parentheses after the experimentally determined values.
Table 6. Lattice parameters of co-crystals using DFT (top) and DFT-D (bottom). All lattice constants are reported in units of Åand % deviation from experimentally determined values are given in parentheses after the experimentally determined values.
Co-Crystal a DFT b DFT c DFT β DFT
(Å)(Å)(Å)
2(SA)·BPE-I13.48 (11.93, +12.99)4.89 (4.87, +0.45)20.71 (20.25, −2.27)99.81 (106.92, −6.65)
2(SA)·BPE-II9.12 (8.76, +4.11)6.48 (6.81, −4.85)24.65 (19.66, +25.38)110.86 (105.34, +5.24)
2(SA)·BPEth9.25 (8.63, +7.18)6.49 (6.86, −5.39)22.35 (19.55, +14.32)91.94 (101.36, −9.29)
2(SA)·BPE-I11.76 (11.93, −1.41)4.89 (4.87, +0.45)19.89 (20.25, −1.08)108.07 (106.92, +1.08)
2(SA)·BPE-II8.73 (8.76, −0.34)6.87 (6.81, +0.91)18.06 (19.66, −8.15)105.91 (105.34, +0.54)
2(SA)·BPEth8.61 (8.63, −0.26)6.88 (6.86, +0.34)17.98 (19.55, −8.02)102.92 (101.36, +1.54)
Table 7. Experimentally determined THz frequencies for the co-crystals 2(SA)·BPE form I and II, and 2(SA)·BPEth, given in units of cm 1 . Values are partitioned to highlight similarities and differences.
Table 7. Experimentally determined THz frequencies for the co-crystals 2(SA)·BPE form I and II, and 2(SA)·BPEth, given in units of cm 1 . Values are partitioned to highlight similarities and differences.
2(SA)·BPE-I29.533.538.844.954.6 68.173.7 92.799.6102.4
2(SA)·BPE-II24.935.4 42.8 60.8 84.0 108.7
2(SA)·BPEth 35.2 43.6 61.6 83.4 110.7
Table 8. Mode-by-mode analysis of the directional components of the IR active response of 2(SA)·BPEth co-crystal using DFT-D for the THz frequency range of 20–110 cm 1 . ω are given in units of cm 1 and all ϵ are unitless. Asterisks are next to the frequencies that have high contributions to ϵ μ and are plotted in Figure 11. The final two rows are the total ionic portion of the dielectric response per direction (x, y, z) per mode μ (denoted as ϵ μ ), and the directionally decomposed electronic contribution of the dielectric response ϵ .
Table 8. Mode-by-mode analysis of the directional components of the IR active response of 2(SA)·BPEth co-crystal using DFT-D for the THz frequency range of 20–110 cm 1 . ω are given in units of cm 1 and all ϵ are unitless. Asterisks are next to the frequencies that have high contributions to ϵ μ and are plotted in Figure 11. The final two rows are the total ionic portion of the dielectric response per direction (x, y, z) per mode μ (denoted as ϵ μ ), and the directionally decomposed electronic contribution of the dielectric response ϵ .
System ω xyz
2(SA)·BPEth38.50.00-0.05
44.5-0.01-
56.0-0.01-
77.7 *0.09-0.05
79.9 *-0.15-
89.8 *0.22-0.06
91.3-0.08-
93.2 *0.02-0.10
96.6 *0.12-0.01
98.9-0.01-
ϵ μ 0.991.570.85
ϵ 3.023.192.71
Table 9. Same description as Table 8, but for 2(SA)·BPE polymorphs I (top) and II (bottom) using DFT-D. Asterisks are next to the frequencies that have high contributions to ϵ μ and are plotted in Figure 12 and Figure 13 for 2(SA)·BPE polymorphs I and II, respectively.
Table 9. Same description as Table 8, but for 2(SA)·BPE polymorphs I (top) and II (bottom) using DFT-D. Asterisks are next to the frequencies that have high contributions to ϵ μ and are plotted in Figure 12 and Figure 13 for 2(SA)·BPE polymorphs I and II, respectively.
2(SA)·BPE-I28.9 *0.00-0.08
33.9 *0.01-0.12
34.9 *-0.15-
48.5-0.01-
50.90.01-0.02
51.1-0.01-
61.50.02-0.02
69.8 *-0.13-
70.3 *0.06-0.02
87.80.01-0.03
109.50.06-0.02
ϵ μ 0.582.191.66
ϵ 2.863.163.14
2(SA)·BPE-II41.5-0.01-
42.3 *0.02-0.28
47.3-0.04-
68.5 *0.00-0.11
69.1 *-0.19-
88.40.07-0.03
90.9-0.02-
92.3 *0.19-0.02
102.9 *0.24-0.02
103.9-0.04-
110.2-0.02-
ϵ μ 0.791.010.64
ϵ 3.443.722.83
Table 10. Experimentally determined ϵ of co-crystals compared to the isotropically averaged ϵ computed using DFT-D. Percent deviation from experiment is in parentheses.
Table 10. Experimentally determined ϵ of co-crystals compared to the isotropically averaged ϵ computed using DFT-D. Percent deviation from experiment is in parentheses.
ExperimentDFT-D
2(SA)·BPE-I6.13 ± 0.194.53 (−26.10)
2(SA)·BPE-II4.89 ± 0.054.14 (−15.34)
2(SA)·BPEth7.82 ± 0.154.11 (−47.44)

Share and Cite

MDPI and ACS Style

Bennett, J.W.; Raglione, M.E.; Oburn, S.M.; MacGillivray, L.R.; Arnold, M.A.; Mason, S.E. DFT Computed Dielectric Response and THz Spectra of Organic Co-Crystals and Their Constituent Components. Molecules 2019, 24, 959. https://doi.org/10.3390/molecules24050959

AMA Style

Bennett JW, Raglione ME, Oburn SM, MacGillivray LR, Arnold MA, Mason SE. DFT Computed Dielectric Response and THz Spectra of Organic Co-Crystals and Their Constituent Components. Molecules. 2019; 24(5):959. https://doi.org/10.3390/molecules24050959

Chicago/Turabian Style

Bennett, Joseph W., Michaella E. Raglione, Shalisa M. Oburn, Leonard R. MacGillivray, Mark A. Arnold, and Sara E. Mason. 2019. "DFT Computed Dielectric Response and THz Spectra of Organic Co-Crystals and Their Constituent Components" Molecules 24, no. 5: 959. https://doi.org/10.3390/molecules24050959

APA Style

Bennett, J. W., Raglione, M. E., Oburn, S. M., MacGillivray, L. R., Arnold, M. A., & Mason, S. E. (2019). DFT Computed Dielectric Response and THz Spectra of Organic Co-Crystals and Their Constituent Components. Molecules, 24(5), 959. https://doi.org/10.3390/molecules24050959

Article Metrics

Back to TopTop