1. Introduction

Terahertz (THz) spectroscopy has been utilized in a diverse number of fields, ranging from security applications¹ to the pharmaceutical industry.² These experimental THz investigations have demonstrated that many compounds have distinct absorption spectra that can be used for their detection, identification, and characterization. While the origins of these absorption features are generally attributed to intermolecular vibrations, intramolecular torsions, or even crystal lattice vibrations, the actual assignment of particular observed spectral features to specific atomic motions is rare. The assignment of these spectral features is crucial to understanding these characteristic THz spectra.

The difficulty in assigning low-frequency THz spectra derives from the very nature of the vibrational motions. Broadly speaking, these motions can no longer be considered as localized atomic motions, like those encountered in the mid-infrared (e.g., O-H stretching modes). The THz frequency vibrational motions must be considered in a global sense where all of the atoms in the entire molecule are participating in the intramolecular or intermolecular motion. This global motion, particularly the intermolecular coordinates, leads to a complete failure of familiar single-molecule-based modeling approaches. A typical calculation treats the molecule of interest in isolation with no environmental interactions present. This is, of course, largely appropriate for gas-phase molecules and some weakly interacting solute-solvent solutions but generally inappropriate for molecular solids.

To address the limitations of isolated-molecule calculations, we have applied periodic boundary condition (PBC) calculations to interpret the THz spectrum of the -crystal form (the most stable of the -, -, -, and -HMX crystal forms)³ of the high explosive octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine (HMX). This energetic material, among others, is a threat agent that can be detected using THz spectroscopy in, for instance, the nondestructive and noninvasive screening of opaque containers.⁴ The experimental THz spectrum of -HMX has been previously reported,^4-7 but there has never been a complete assignment of the spectral features. For the rational use of the -HMX THz spectrum in threat detection databases, a thorough understanding of the spectrum is required.

2. Methods

The experimental room-temperature THz spectrum of crystalline -HMX was obtained from TeraView Limited (Cambridge, U.K.).⁵ Density functional theory (DFT) calculations were performed using the programs Gaussian03⁸ (isolated-molecule) and DMol³ (version 3.2)^9,10 (isolated-molecule and solid-state). Gaussian03 was employed for HMX geometry optimizations and normal-mode analyses with the B3LYP^11,12 hybrid density functional and the 6-31G(d,p)¹³ Gaussian-type basis set. DMol³ isolated-molecule and solid-state calculations were performed with "fine" grid sizes (corresponding to a k-point separation of 0.04 Å^-1) and convergence criteria (program options), the DNP (double numerical with d and p polarization) basis set (comparable to a 6-31G(d,p) Gaussian-type basis set),⁹ and one functional from among the PWC¹⁴ and VWN¹⁵ local density approximation (LDA) and BLYP,^16,17 BOP,¹⁸ PW91,¹⁹ HCTH,²⁰ BP,^16,21 PBE,^22,23 RPBE,²⁴ and VWN-BP^15,16,21 generalized gradient approximation (GGA) density functionals. These functionals represent those available in the current version of the DMol³ program. The presented calculations are referred to only by their functional in the text, with the DNP basis set used in all cases. Unlike many plane-wave DFT packages, DMol³ does not optimize lattice constants, meaning the molecules within the unit cell were optimized within the cell parameters specified by the room-temperature X-ray diffraction study.³ Crystal cell parameters are as follows: space group P21/c (Z = 2), a = 6.54 Å, b = 11.05 Å, c = 8.70 Å, = 90.0, = 124.3, = 90.0.

3. Results

3.1. Comparison of Experimental and Theoretical HMX Structural Data. Crystallographic and calculated (isolated-molecule and solid-state) HMX structural data are provided in Tables 1 (bond lengths) and 2 (bond angles). The calculated bond length differences from experiment are shown diagrammatically in Figure 1 for comparative purposes. Both the isolated-molecule and solid-state optimizations yield C_i symmetry HMX minima (with the exception of the PWC and VWN isolated-molecule calculations, which deviate only slightly to C₁ symmetry), consistent with the crystal data and previous relative energy comparisons of HMX conformations.²⁵ The best bond length agreements occur in the PWC, VWN, HCTH, and B3LYP optimizations by isolated-molecule root-mean-square deviation (RMSD) comparisons. This same trend continues in the solid-state calculations, where the LDA results are second only to the HCTH results in the accurate reproduction of crystal bond lengths. As expected from the inclusion of neighboring molecules, considerable improvement in bond angle agreement is found in the solid-state calculations relative to their isolated-molecule counterparts, with the solid-state angle agreement surpassing the accuracy of the isolated-molecule B3LYP results in all LDA and GGA cases. It is noteworthy that the functionals with the best experimental bond length agreement are also the ones with the worst bond angle RMSD agreement, although these differences occur over a narrow range.

Figure 1 Calculated bond length differences from experiment for the isolated-molecule (left) and solid-state (right) DMol3 minimum energy geometries. The PWC/VWN results were obtained by taking the average bond lengths for the C1 symmetry structures. Some density functionals yield nearly identical trends and are grouped in the figure for clarity.

The solid-state GGA and LDA -HMX geometry optimizations, through their inclusion of crystal packing interactions and weak intermolecular hydrogen-bonding interactions, improve the agreement of the calculated HMX geometry with experiment. All but three (BLYP, BOP, RPBE) functionals either match or surpass the isolated-molecule B3LYP calculations in accuracy by RMSD comparisons. The solid-state calculations still generally overestimate bond lengths as compared to observed crystal distances. The GGA bond length overestimation is larger than in the LDA results. The only cases where the LDA bond lengths are shorter than observed are the C-N distances in the eight-membered ring.

Overlay of the isolated-molecule and solid-state HMX structures reveals the directions and magnitudes of structural changes in the molecule due to crystal cell interactions. The overlay of the HCTH structures (those with the best crystal bond length agreement) is shown in Figure 2. The local features (bond lengths, angles) are found to change very little, while the pendant NO₂ groups are all shifted toward the isolated-molecule inversion center in the solid-state structure. The features of the crystal cell not accounted for by this structural comparison are the intermolecular hydrogen-bonding interactions that couple each HMX molecule to its neighbors. Assigned hydrogen-bonding interactions are provided in Table 3. Intermolecular hydrogen bonding in the -HMX crystal occurs well within the range of weak hydrogen bonds observed in other molecular crystals.²⁶ Additionally, five symmetry-unique intramolecular hydrogen-bonding interactions are found in the isolated-molecule HMX calculations in the range of 2.10-2.60 Å. Of note in the isolated-molecule calculations are the reductions in the O₃-H_b and O₄-H_c distances in the LDA calculations relative to the other methods. These result from the reduced N₁-N₂ and N₃-N₄ bond lengths in these same LDA calculations (Table 1).

Figure 2 The packing arrangement of the -HMX crystal cell (left), an overlay of the HCTH/DNP isolated-molecule and solid-state geometries (center), and the labeling scheme for the HMX molecule used in bond length and angle assignments (right). Vertical (V) and horizontal (H) NO2 groups are differentiated in the overlay and labeled figures. Crystal and overlay figures were rendered with VMD.31

In the crystal cell, intermolecular hydrogen-bonding interactions form that approach and, in some instances, occur within the range of the intramolecular distances (Table 3). Among the GGA solid-state calculations, all but the O₃-H_b distances are found to increase upon geometry optimization relative to the isolated-molecule distances. Atoms O₃ and H_b reside below the plane of the four methylene carbons and are brought into closer contact due to the bending of the N₄-bound NO₂ group toward the molecular center in the crystal environment. In the LDA calculations, both the O₃-H_b and O₂-H_a distances are found to be reduced slightly. In the crystal cell, intermolecular distances to the H_a and H_b atoms (from oxygen atoms in nearest-neighbor HMX molecules) occur at distances within the intramolecular regime, which is noteworthy both for the shorter distances and for the fact that the H_b atom is found to interact more closely with its intramolecular acceptor (O₃) in the GGA cases (both H_a and H_b reside closer to their acceptors in the LDA cases) than the isolated-molecule geometry optimizations. The crystal environment also positions other HMX molecules near H_d, the only atom to which short (<2.60 Å) intramolecular hydrogen-bonding distances could not be assigned in the isolated-molecule calculations.

3.2. IR Intensity Calculations by Difference-Dipole Methods. The assignment of THz features in the -HMX spectrum is complicated by (1) the differences in the treatment of crystal packing interactions among the various density functionals, (2) the lack of lattice constant optimizations in DMol³, which also affects the treatment of intermolecular interactions in the calculations, and (3) the presence of both internal and external modes in the THz region. For these reasons, both frequency and intensity data are important for use in accurate assignments. Currently, DMol³ cannot calculate solid-state IR intensities, nor can static electric fields be applied in solid-state calculations for calculating IR intensities by atomic polar tensor methods. A means for indirectly calculating IR intensities is available in the solid state through the calculation of unit cell dipole moment differences between the ground and vibrational states based on Mulliken²⁷ and Hirshfeld²⁸ atomic charges. The solid-state IR intensities reported in this work are calculated from the change in dipole moments for the unit cell that result from atomic displacements along each normal-mode coordinate (d/dQ). A previous study reports the use of this difference-dipole method for the assignment of -sheet vibrations in oligopeptides based on Mulliken charges.²⁹ While this approach, in principle, is expected to yield reasonable IR intensities due to the dependence of the intensity on the square of the transition dipole moment, it is important to first validate this approach external to the solid-state calculations. It is also known that Hirshfeld analyses are more useful than Mulliken analyses for explaining molecular properties,³⁰ and so the differences between difference-dipole calculations using both charge partitioning methods is of some importance. To test this method, the difference-dipole technique was first applied to the isolated HMX molecule, where IR intensities and molecular dipoles are directly calculated in the normal-mode analysis.

The scaling of difference dipoles (all between 0.0 and 0.1 D) to intensities (0.0-30.0 km/mol) was performed by applying a constant multiplicative term to the square of the difference dipoles and adjusting this term as to minimize their RMSDs with respect to the calculated molecular intensities of the sub-200 cm^-1 modes (in the PWC and VWN calculations, this range was increased to 220 cm^-1 to include two modes that occur in the GGA functionals below 200 cm^-1). The RMSDs of the difference-dipole calculations are provided in Table 4. The histograms in Figure 3 show the IR intensity fits for the C₁ symmetry PWC results (the PWC and VWN results are nearly superimposable) and the C_i symmetry VWN-BP results, which are representative of the eight GGA functionals. The two additional modes in the PWC/VWN cases are a result of the absence of inversion symmetry in the optimized molecules and the reporting of only those modes with intensities larger than 0.00 km/mol.

Figure 3 Isolated-molecule DMol3 IR intensities (km/mol, black) and scaled calculated IR intensities (Debye2) for the PWC/DNP and VWN-BP/DNP structures as determined by the difference-dipole method using the calculated molecule dipoles (red) and dipoles calculated from Mulliken (green) and Hirshfeld (blue) charges. The two additional modes in the PWC/DNP histogram arise from the lowered symmetry in the LDA calculations.

RMSD minimization leads to reasonable agreements in both the absolute intensities and trends in intensities in the sub-200 cm^-1 region. The Mulliken dipole RMSDs are smaller than the Hirshfeld values for all but the BLYP and BOP functionals. The Mulliken and Hirshfeld dipole differences for the lowest-frequency IR-active modes are calculated too small in all cases, but it is clear that trends in intensity are correct across all methods. The Mulliken and Hirshfeld difference dipoles are noteworthy for their good agreement with one another. This agreement can be attributed to the use of the dipole differences between two states and not the actual charges themselves in any one state. Over small displacements, the changes in atomic charge by either partitioning method would be expected to be similar, leading to the consistent difference-dipole results throughout.

The results obtained from the dipole treatment of molecular intensities indicate that similar treatments in the solid-state should yield results accurate enough to prove useful in mode assignments. The results of these calculations are provided in the next section.

3.3. Comparison of Experimental and Simulated THz Spectra. The experimental THz spectrum of -HMX (3-120 cm^-1) is provided in Figure 4 with the calculated isolated-molecule and solid-state spectra. The intensities in the solid-state spectra result from the Mulliken difference-dipole analysis. The experimental spectrum in this region contains three prominent absorptions at 58.5, 82.2, and 95.4 cm^-1. The isolated-molecule normal-mode analyses with the GGA functionals yield five IR-active modes in the 0-120 cm^-1 region, with all eight functionals producing the same relative positions and intensities of the lowest three and the highest fifth modes. The frequency of the fourth mode varies between 74 (RPBE) and 105 (HCTH) cm^-1, with the other functionals generating a tighter grouping between 82 and 89 cm^-1. The LDA functionals are in good agreement with one another in the isolated-molecule calculations, with the fifth dominant IR-active mode (all modes have some IR activity in the C₁ structure) calculated at 140 cm^-1 by both methods and, therefore, beyond the reported experimental results. The calculated isolated-molecule normal modes (both IR-active and inactive) in this region are provided in Table 5. All isolated-molecule normal-mode analyses display little agreement with the THz spectrum of the crystalline material, as apparent from the lack of frequency agreement with the three experimental maxima and the calculated positions of the two most intense IR-active modes below the first apparent THz peak at 58.5 cm^-1.

Figure 4 Comparison of the experimental and simulated THz spectra of -HMX from 3 to 120 cm-1. Isolated-molecule B3LYP/6-31G(d,p) (green) results are provided for comparison. Isolated-molecule [functional]/DNP predictions are shown in red. Solid-state [functional]/DNP predictions are shown in blue. The experimental THz spectrum (black) for -HMX is provided for each data set.

The results of the solid-state normal-mode analyses for all GGA and LDA functionals are shown (with their calculated IR intensities) overlaid with the measured THz spectrum and isolated-molecule peaks in Figure 4 and are listed in Table 6. While variety in peak positions and calculated intensities is greater among the solid-state calculations, it is evident that inclusion of the crystal cell interactions results in significantly improved agreement with the observed THz spectrum in all cases. Both LDA functionals predict IR-active modes within close proximity of THz peaks. A range of agreement can be found from among the GGA functionals, with the VWN-BP, BP, and PW91 calculations yielding excellent THz peak position agreement. The eighth IR-active crystal mode occurs above the 120 cm^-1 THz cutoff in the BLYP, BOP, PBE, RPBE, and HCTH calculations and is grouped within a narrow range in the GGA series (116-126 cm^-1) where the experimental THz spectrum rises to what may be a nearby absorption maximum. The LDA functionals place this eighth mode at 99 cm^-1. The ninth IR-active mode is predicted to lie above the 120 cm^-1 THz cutoff in all the GGA calculations, while the LDA calculations place the ninth mode approximately 3 cm^-1 above the eighth crystal mode.

3.4. THz Mode Assignments From the VWN-BP/DNP Simulated Spectrum. From among the density functional comparisons, the VWN-BP simulated solid-state THz spectrum agrees best in terms of both peak positions and relative calculated IR intensities by the difference-dipole method. It is with these results that the assignment of the THz spectrum to specific molecular motions is performed. The nine lowest-frequency molecular vibrational modes are shown in Figure 5. These low-frequency modes include NO₂ pendulum motions and ring-twisting motions (modes 1 and 2), large-amplitude ring stretches and twisting modes (modes 3, 4, and 5), and in- and out-of-phase NO₂ rotor modes for both the vertically (modes 6 and 7) and horizontally (modes 8 and 9) oriented NO₂ fragments (for orientation, see Figure 5). The distinct lack of agreement between the THz peaks and isolated-molecule modes in Figure 6 indicates, as expected, that these large-amplitude motions are altered in the crystal environment.

Figure 5 The nine lowest-energy normal modes (and equilibrium structure) of the HMX molecule with frequency assignments based on the isolated-molecule VWN-BP/DNP calculations. The numerical labeling scheme is used for the identification of solid-state mode components in Figure 6. Displacements along the normal mode are scaled by a factor of 5 for clarity. The vertical (V) and horizontal (H) NO2 groups discussed in the text are identified in the equilibrium illustration. Figures were rendered with VMD.31

Figure 6 The isolated-molecule (top) and solid-state VWN-BP/DNP normal modes (0-120 cm-1) with solid-state intensities provided by the Hirshfeld (middle) and Mulliken (bottom) difference-dipole results. The isolated-molecule results show the Mulliken and Hirshfeld intensities offset from the actual mode frequencies (shown with the IR intensities) for clarity. The internal mode assignments (numerical labels) are as provided in Figure 5. External mode labels are as follows: OT1 = optical translation along the AC crystal plane; OT2 = optical translation along the AB crystal plane; OT3 = optical translation along the BC crystal plane.

A crystal unit cell containing M molecules with N atoms contains 3N-6M internal modes (those modes associated with intramolecular motions), 6M-3 external modes (those modes associated with relative motions between the M molecules, such as rotations and translations), and 3 acoustic modes. In the -HMX crystal (Z = 2), the nine external modes correspond to three optical translational (IR-active) modes (the two molecules moving in opposite directions along the [x,y,z] axes) and six optical rotational (IR-inactive) modes, which occur as in-phase and out-of-phase motions of the two molecules in the crystal unit cell. Because the energies and relative motions of large-amplitude (low-frequency) internal modes and external modes are dependent on crystal packing and intermolecular interactions (such as ionic and hydrogen-bonding interactions), these two different types of modes can be intermixed in the low-frequency region. This intermixing, coupled with the changes to molecular mode frequencies due to crystallographic forces, is the primary reason why isolated-molecule calculations are of limited use in low-frequency vibrational assignments.

In the absence of strong electrostatic interactions to couple neighboring molecules (and, subsequently, their vibrational modes), the low-frequency vibrations of a molecular crystal (those not due to external modes) can be assigned based on the relative motions of isolated molecules with respect to their internal modes. In this absence of strong interactions, a single molecular mode occurs as in-phase and out-of-phase combinations of the same normal mode displacements in crystal cells containing two identical molecules (such as HMX). The degree to which the crystal packing interactions alter one or both of these combinations is reflected in both the frequency shift relative to the isolated molecular mode positions and the relative splitting of these in- and out-of-phase solid-state combinations.

The solid-state vibrational modes in -HMX can be readily assigned to either internal (molecular) mode combinations or external modes (Table 7). The THz spectrum and VWN-BP spectra are shown and labeled in Figure 6 using both the Hirshfeld (middle) and Mulliken (bottom) difference-dipole intensities. The normal-mode assignments reveal considerable shifting of internal modes in the crystal cell and the prominence of the three optical translational modes in the THz region. It is these three modes that isolated-molecule calculations are incapable of reproducing. The six rotational modes are all predicted to be IR-inactive and occur in the 40-90 cm^-1 range (Table 7). The IR-active modes below 120 cm^-1 are found to include the in- and out-of-phase combinations of the vertical NO₂ pendulum motions (mode 1 in Figure 5) split by 25.5 cm^-1, the out-of- and in-phase combinations of the A_u symmetry ring-twisting modes (mode 5 in Figure 5) split by 1.0 cm^-1, and, at the very edge of this range, the out-of-phase combination of the horizontal NO₂ pendulum motions (mode 2 in Figure 5) that lies 7.4 cm^-1 below its in-phase combination. The remaining combinations of the low-frequency molecular modes occur beyond the 120 cm^-1 THz limit. From among the molecular modes available for assignment to this THz region, the isolated-molecule calculation includes six modes (including IR-active modes 7 and 9) that do not contribute to features in the THz spectrum according to the VWN-BP calculations.

4. Discussion

4.1. Density Functional Selection and Molecular Geometry. Inclusion of the crystal environment is found to significantly improve the HMX molecular geometry agreement between theory and experiment. In RMSD comparisons, solid-state DFT calculations are always found to improve over their molecular counterparts, with the HCTH calculations yielding the best overall bond length agreement and RPBE calculations yielding the best overall bond angle agreement with experiment. The only two bond lengths that are not corrected for in the solid-state calculations are the N₁-O₁ and N₁-O₂ distances of the vertical NO₂ groups (for orientation, see Figure 2). One of the useful insights to come from the functional comparisons is the realization that the degree of agreement between theory and experiment depends on the experimental feature being considered. While the HCTH results yield the best bond length agreement, it is clear from the vibrational analyses that this functional is less accurate than others in the GGA series at reproducing major spectral features. The sensitivity of the normal modes to the choice of functional is, from these results, more significant than the range in RMSDs in the LDA/GGA series.

The observed agreement in the solid-state DMol³ results with the crystal distances surpasses the results of the more computationally demanding isolated-molecule calculations with the B3LYP hybrid density functional. From a methodological perspective, it is not the better agreements of the GGA and LDA functionals over the isolated-molecule B3LYP results that are noteworthy here, but the improvement in structural agreement to experiment that comes from including intermolecular interactions in the GGA and LDA computational analyses. While crystal packing alone is an important factor in determining molecular geometry for flexible molecules, the -HMX crystal is found to contain numerous intermolecular weak hydrogen-bonding interactions between molecules that, with crystal packing, alter the intramolecular interactions calculated for the isolated molecule. These interactions are chemically weak but lead to significant changes in the low-frequency internal modes from their isolated-molecule origins.

4.2. Normal-Mode Analyses and Solid-State Assignments. The assignment of low-frequency vibrational modes in molecular solids from isolated-molecule calculations is certainly demonstrated to be incorrect for HMX. The majority of the sub-120 cm^-1 isolated-molecule modes for HMX do not occur in the available THz range for -HMX according to the VWN-BP calculations. Despite variability among the available density functionals, the crystal cell normal-mode calculations are always found to be improvements over their isolated-molecule counterparts. This is certainly because the solid-state calculations include crystal packing, intermolecular hydrogen-bonding interactions, and restriction of motion in the normal modes due to these structural interactions. This restriction of motion is the time-dependent feature of normal-mode analyses and vibrational spectroscopy that combines with the time-independent quality of geometry optimizations and diffraction data to provide the more complete picture of the crystal environment. A detailed analysis of the VWN-BP solid-state calculations reveals the degree to which the crystal environment changes the isolated-molecule vibrations. Several absorption features in the THz spectrum are assigned to optical translational (IR-active) modes that cannot be accounted for in isolated-molecule calculations since these motions are solely the result of the relative motions of the two HMX molecules in the crystal cell. The isolated-molecule modes undergo both shifting to higher frequency and, due to the presence of in- and out-of-phase combinations of these individual modes in the Z = 2 crystal cell, significant mode splitting when the relative motions of one combination are less restricted than the oppositely phased motions in the other. Within the 120-cm^-1 range of the THz experiment, the in-phase combination of the lowest-frequency isolated-molecule pendulum motions (mode 1 in Figure 5) is shifted 40 cm^-1 to higher frequency, while the out-of-phase combination is 25.5 cm^-1 higher still in energy. In contrast, the asymmetric twisting mode at 62.6 cm^-1 (mode 5) is shifted by only 27 cm^-1 and split between its in- and out-of-phase components by 1 cm^-1.

4.3. Mulliken and Hirshfeld Difference-Dipole Solid-State Intensities. The calculation of crystal cell dipole differences along each normal coordinate yields IR intensity trends that can, provided the frequencies themselves accurately reproduce the positions of major peaks in the THz spectrum, aid greatly in solid-state mode assignments. Both the Mulliken and Hirshfeld difference-dipole intensities in the VWN-BP solid-state calculations reproduce the positions of major peaks up to 120 cm^-1 in -HMX and, when plotted with a Lorentzian line shape, accurately account for major shoulder features where the prominence of minor absorptions are diminished by local nearby maxima. The only notable difference between the Mulliken and Hirshfeld intensity spectra (Figure 6) is the relative calculated intensities of the THz peak at 82 cm^-1, which find the Hirshfeld results predicting this mode too high relative to both the Mulliken spectrum and experiment. These two peaks are found to contain two closely spaced crystal cell modes each, making the relative intensities of these two modes sensitive to the accuracy of the charge analyses in four modes. Here, the Mulliken intensities are found to reproduce the experiment through both a larger difference in the two most intense modes in this region ("OT3" and "5-5" in Figure 6) and the prediction of negligible intensity in the two other contributing modes ("1 - 1" and "5 + 5" in Figure 6). This difference-dipole approach has proven useful in the -HMX analysis and is being further considered in other solid-state studies of molecular crystals.

5. Conclusions

The THz spectrum of -HMX has been analyzed with isolated-molecule and solid-state DFT calculations. Isolated-molecule calculations are incapable of reproducing the observed spectral features. The inclusion of the crystal environment in the DFT calculations greatly improves both the crystallographic structural agreement relative to the isolated-molecule calculations and the theoretical agreement with the measured THz spectrum. The available THz spectrum of -HMX (3-120 cm^-1) is assigned based on the results of solid-state normal-mode analyses at the VWN-BP/DNP level of theory. Analysis of the individual modes reveals that major features in the THz spectrum are due to optical translational motions of molecules in the crystal cell and that the majority of isolated-molecule normal modes predicted to lie below 120 cm^-1 are, in fact, shifted to much higher frequency. This is due to crystal packing and intermolecular interactions, further indicating the limitation of employing isolated-molecule calculations for assignment purposes. Reasonable IR intensity predictions can be made in the solid-state by using difference-dipole calculations from Mulliken and Hirshfeld charge analyses in DMol³, with Mulliken charge analyses generally providing better agreement in the -HMX studies. Further investigations of the application of solid-state DFT methods to the analysis of the THz spectra of energetic materials are ongoing.

Acknowledgment

The authors gratefully acknowledge the support of the National Science Foundation (PHY-0442188) and the ACS Petroleum Research Fund (40610-G6). The authors thank TeraView Limited (Cambridge, UK) for providing the experimental THz spectrum of -HMX and the reviewers for their insightful comments.