J. Phys. Chem. C 2007, 111, 2787-2796
2787
Ab Initio Study of Compressed 1,3,5,7-Tetranitro-1,3,5,7-tetraazacyclooctane (HMX), Cyclotrimethylenetrinitramine (RDX), 2,4,6,8,10,12-Hexanitrohexaazaisowurzitane (CL-20), 2,4,6-Trinitro-1,3,5-benzenetriamine (TATB), and Pentaerythritol Tetranitrate (PETN) Edward F. C. Byrd* and Betsy M. Rice US Army Research Laboratory, AMSRD-ARL-WM-BD, Aberdeen ProVing Ground, Maryland 21005-5066 ReceiVed: March 22, 2006; In Final Form: NoVember 15, 2006
Using the PW91, PBE, and LDA density functional theories (DFT), we have calculated crystal structures for five energetic molecular crystals over a range of experimental pressures. These crystals are 1,3,5,7-tetranitro1,3,5,7-tetraazacyclooctane (HMX), cyclotrimethylenetrinitramine (RDX), 2,4,6,8,10,12-hexanitrohexaazaisowurzitane (CL-20), 2,4,6-trinitro-1,3,5-benzenetriamine (TATB), and pentaerythritol tetranitrate (PETN). Both PW91 and PBE generally overestimate volumes relative to experimental values, while LDA underestimates crystal volumes when compared to experiment. However, the inaccuracy diminishes as pressures are increased. In particular, PW91 and PBE volumes approach experimental values at pressures greater than 6-7 GPa. Furthermore the PW91 and PBE volumes appear to converge to the same value as pressures increase, regardless of the size of the planewave basis set. We have also demonstrated that for systems such as these, care should be taken to ensure convergence in DFT calculations. We emphasize that the mistreatment of van der Waals forces by DFT will have unforeseen consequences on all crystal calculations for weakly bound organic molecules, and therefore caution should be employed whenever interpreting results obtained from the current DFT functionals available in solid-state codes.
1. Introduction The often prohibitive time and cost associated with the synthesis, testing, and fielding of a new energetic material1 have driven the advancement of computationally assisted design of new materials. Significant efforts have been directed toward the development of theoretical models that will predict properties of an energetic material that will allow an estimate of its potential performance or possible hazard.2 One of the key properties of an energetic material required to assess potential performance in a weapons system is the crystal density, since detonation velocities and pressures are proportion to powers in the density. Therefore, it is crucial that this property be accurately predicted. Various modeling schemes exist to predict crystal densities of materials, ranging from classical molecular simulation methods (e.g., molecular dynamics, molecular packing) to semiempirical and quantum mechanical methods. The classical molecular simulation methods are limited by the dependence on empirically based descriptions of interatomic interactions and are parametrized to describe specific chemical systems. Quantum mechanical theories, however, are inherently transferable and thus are significantly more desirable in predicting novel systems. Indeed, numerous scientists within the energetic materials community have embraced quantummechanically based methods in their research efforts.3-41 There are a variety of quantum mechanical theories that can treat condensed phase materials, including Hartree-Fock (which does not include electron correlation), or more sophisticated treatments, including Density Functional Theory (DFT), MøllerPlesset (MP) Theory, and Quantum Monte Carlo (QMC). Unfortunately, the computational resources required for MP and QMC are so much greater than those required for DFT * Corresponding author. Phone: 410-306-0729.
10.1021/jp0617930
calculations, studies using these methods are usually limited in scope. For this reason, applications of DFT to condensed phase materials have significantly increased over the past few years.42-45 DFT has been routinely applied to the study of large polyatomic molecules over the past decade due to its demonstration of reasonable accuracies coupled with computational efficiency. Likewise, DFT properly describes metallic and covalently bound condensed phase systems. However, DFT does not demonstrate similar successes for condensed phase systems composed of weakly interacting species (such as organic molecular crystals). It is well-known that conventional formulations of DFT do not include dispersion interactions,46-47 although efforts are being made to remedy this deficiency.48-60 Unfortunately, dispersion interactions appear to be major contributors to the binding forces within many conventional energetic materials, particularly neutral CHNO crystals, as we and others have demonstrated in earlier investigations.10,12,61 Although it is disappointing that at this time DFT appears to be inadequate for studies of energetic molecular crystals at low pressures, there is evidence that DFT can better represent weakly bound molecular systems under sufficient degrees of compression.10,62 Since we are also interested in the theoretical characterization of energetic materials at extreme pressure, it is useful to establish the predictive capability of various DFT methods for such materials under large degrees of compression. Toward this end, we present calculations of the crystallographic and intramolecular parameters for five energetic molecular crystals under varying degrees of compression. Different forms of density functional theories, pseudopotentials, and planewave basis set will be employed in order to eliminate possible sources of error.
This article not subject to U.S. Copyright. Published 2007 by the American Chemical Society Published on Web 01/24/2007
2788 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Byrd and Rice
TABLE 1: Experimental Lattice Parameters and Space Groups at Ambient Pressure a (Å)
b (Å)
c (Å) R (deg) β (deg) γ (deg) space group
RDX71 13.182 11.574 10.709 90.0 90.0 90.0 -CL-2072 8.852 12.556 13.386 90.0 106.82 90.0 β-HMX73 6.54 11.05 8.7 90.0 124.3 90.0 (i) PETN74 9.386 9.386 6.715 90.0 90.0 90.0 75 TATB 9.01 9.028 6.812 108.58 91.82 119.97
Pbca P21/n P21/c P421c P1
2. Computational Details Calculations of the crystal and intramolecular parameters for HMX, RDX, CL-20, 2,4,6-trinitro-1,3,5-benzenetriamine (TATB), and pentaerythritol tetranitrate (PETN) were performed using the two Generalized Gradient Approximation (GGA) density functional theories Perdew-Wang 91 (PW91)63 and the PerdewBurke-Ernzerhof (PBE)64 functionals as implemented in the Vienna Ab-Initio Package (VASP).65 The Local Density Approximation LDA66 functional was also employed as an alternate method for TATB, PETN, and RDX. The Vanderbilt ultrasoft pseudopotentials (USP)67 and Monkhorst-Pack k-point generation method were used for all VASP calculations unless otherwise stated. Projected augmented wavefunction (PAW)68,69 pseudopotentials were used to check for possible pseudopotential-induced errors for the TATB and PETN crystals. The influence of pseudopotentials was also tested using normconserving Troullier-Martins pseudopotentials on PETN with the PBE functional using the PWSCF program package.70 K-point meshes employed were: 2 × 1 × 2 for HMX, 1 × 1 × 1 for both RDX and CL-20, 2 × 2 × 2 for TATB, and 2 × 2 × 2 for PETN. The Troullier-Martins PBE PWSCF calculations used solely the gamma point. We also explored the effects of basis set size on results for the compressed crystals, similar in spirit to our earlier study investigating the performance of GGA DFT when applied to four energetic molecular crystals (three of which are studied here) at zero pressure.61 That study showed that the predicted volumes were very sensitive to the size of the planewave basis set employed (ranging from 280 to 800 eV). Lattice vector lengths deviated from experiment by as much as 1.0 Å with the degree of error dependent on the basis set size. Interestingly, the degree of error did not necessarily decrease with increasing basis set size. Rather, Byrd et al.61 observed a kinetic energy cutoff (∼400 eV) that fortuitously produced near-exact agreement with experimental lattice vector lengths for the four systems. In the current study, we have employed planewave kinetic energy cut-offs (Ecut) of either 396 eV and/or 545 eV for both the Vanderbilt ultrasofts and PAW pseudopotentials. The smaller cutoff was chosen since it is very near the kinetic energy cutoff that produced results in near exact agreement with experimental lattice vectors in the earlier study (and is the default kinetic energy cutoff in VASP). The 545 eV kinetic energy cutoff was selected since the lattice vectors were shown to be near-converged with respect to basis set size by this cutoff.61 Starting from the zero pressure (ambient conditions) experimental crystal structures, all ionic and cell degrees of freedom were allowed to relax concurrently and without enforcing symmetry restrictions at various imposed isotropic pressures. Table 1 gives the ambient pressure experimental values used for initial cell parameters in all calculations. Except for the Trouillier-Matins PBE calculations, all theoretical equilibrium crystal structures were obtained through analytic gradient techniques implemented in VASP. Using the default convergence criteria in VASP, the electronic energies were converged to 2.0 × 10-6 eV, while the structures were converged once the difference in free energy between gradient steps was less
than 2.0 × 10-5 eV. We also calculated changes in the mean atomic displacements and residual forces and changes in lattice parameters between the next-to-last and final steps in the geometry optimizations, in order to better establish the degree of convergence of the results. The largest mean error of the components of the atomic forces are ∼1.0 × 10-7 eV/Å or smaller, and the average errors in the atomic displacements are ∼1.0 × 10-4 Å or smaller. Changes in lattice vector lengths between the two optimization steps are no greater than 2.2 × 10-4 Å. This suggests that the convergence criteria used in these calculations are adequate. As these molecular crystals are insulators, setting the smearing width to 0.0001 eV in the Methfessel-Paxton method76minimized electron smearing. To accelerate convergence, the RMM-DIIS method was employed, except for certain calculations where it was necessary to use conjugate gradients instead. The Troullier-Martins PBE calculations of equilibrium crystals structures of PETN utilized PWSCF.70 In these calculations, the electronic energies were converged to 1.0 × 10-9 Rydberg, while the total energy convergence was 1.0 × 10-6 Rydberg. The final converged geometries were obtained through a damped molecular dynamics calculation as coded in PWSCF.70 3. Results and Discussion Although we have performed calculations on the five crystals, only three have extensive experimental information over a large pressure range. The majority of the experimental data available are for the TATB,77 HMX,78 and PETN79 crystals. RDX78 and CL-2080 exhibit pressure-induced structural changes at relatively lower pressures than the aforementioned systems. Since these systems lack measured crystallographic information against which to compare at the higher pressures, the majority of the results in this study shall focus on the TATB, HMX, and PETN crystals. 3.1. Ultrasoft Pseudopotential PW91. The first method investigated was the GGA-based PW91 using Vanderbilt ultrasoft pseudopotentials (USP) with both the 396 and 545 eV planewave basis sets. The 396 eV basis set is the default basis set from the VASP program package, as this is the energy cutoff for the oxygen atom ultrasoft pseudopotential. It is also very close to the kinetic energy cutoff that was used to produce the zero-pressure results reported in Byrd et al.,61 wherein the errors appear to fortuitously cancel, thus producing good agreement with experiment. The 545 eV basis set was also chosen as it was shown in Byrd et al. to be the point where convergence due to basis set size began, yet exhibited large errors in predicted volumes.61 Figures 1-5 compare 545 eV predicted volumes with experimental values for the five crystals as functions of pressure. Due to the unavailability of detailed crystallographic information for RDX and CL-20 at pressures greater than 4 and 2.5 GPa, respectively, we did not perform calculations at higher pressures for these two systems. Overall, the predicted volumes of the five crystals are larger than experimental values over the entire pressure ranges explored. While the behavior at low pressure is similar to that predicted for molecular crystals at zero pressure, the theoretical results more closely agree with experiment as pressure increases. The 396 eV Ecut yields volumes which are smaller than experiment at zero pressure; however, as the pressure increases the theoretical volumes, for the most part, are larger than experimental values. We shall investigate this reduction in error with increasing pressure in more detail for HMX, PETN, and TATB. For HMX, PETN, and TATB, the 396 and 545 eV basis sets converge to the same values at pressures greater than 6 GPa
Structures of HMX, RDX, CL-20, PETN, and TATB
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2789
Figure 1. Variation of volume and lattice dimensions with pressure for RDX. Experimental values taken from ref 78.
Figure 2. Variation of volume and lattice dimensions with pressure for CL-20. Experimental values taken from ref 80.
(approximately 80% compression). For the 396 eV basis set, the zero pressure theoretical crystal volumes undershoot the experimental results and then proceed to overpredict the experimental volumes with increasing pressure. The zero pressure 545 eV results are substantially larger than the experimental values but the deviations of the predictions from experiment decrease with increasing pressure. Linear and volume compression data are often presented in the form of X/X0, where X denotes volume or cell axis. This type of representation using theoretical results could be misleading if one incorrectly assumes that the X0 is in good agreement with measured values. To illustrate, Figure 6 shows V/V0 versus pressure for PETN where V0 is either the corresponding theoretical volume at zero pressure (Figure 6a) or the experimental volume at zero pressure (Figure 6b). As evident in Figure 6a, employing the poorly predicted theoretical V0 in the volume compression curve produces extremely poor agreement with the experimental curve. However, using the experimental volume at zero pressure in such a curve (Figure 6b) produces fairly good agreement with experiment at the higher pressure, with the larger error at the lowest pressure. Since we have shown that there are such large errors in the theoretical volume at zero pressure, we urge caution in representing linear and volume compression data in the form of X/X0.
In order to ensure ourselves that the 545 eV kinetic energy cutoff was sufficient for convergence, the 645 eV kinetic energy cutoff was used at two pressures for both TATB and PETN. The differences in volume between the 545 and 645 eV kinetic energy cutoffs for TATB at 3.5 and 7.0 GPa were approximately 1.0 Å.3 For the PETN crystal at 5.6 and 10.4 GPa, the differences between the basis set cutoffs were approximately 1.9 Å.3 With these small differences, it is apparent the 545 eV basis set is adequate for convergence. In order to test k-point convergence, we performed optimizations for β-HMX and PETN at 0 GPa using a 4 × 4 × 4 k-point grid and PW91 at 545 eV. We did not perform the same tests for RDX and CL-20 due to the exorbitant computational demands for these systems using the larger k-point grid. The resulting volumes were smaller than those calculated using the smaller k-point grids by 12.4 and 9.5 Å3 for PETN and β-HMX, respectively. As for the earlier results, these correspond to unsatisfactorily large deviations from experiment (16.3% and 12.0%, respectively). As the results were not significantly improved with these increases in k-point grid, we expect that only similarly small improvements would occur for predictions at higher pressures using larger k-point grids. 3.2. Ultrasoft Pseudopotential LDA. The second method investigated was LDA with Vanderbilt ultrasoft pseudopotentials
2790 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Byrd and Rice
Figure 3. Variation of volume and lattice dimensions with pressure for HMX. Experimental values taken from ref 78.
Figure 4. Variation of volume and lattice dimensions with pressure for TATB. Experimental values taken from ref 77.
and the 545 eV planewave basis set. The systems studied using this method are RDX, TATB, and PETN. Due to the computational demands for the RDX crystal, only two points were calculated, and as we will show, it was unnecessary to compute volumes at additional pressures. We did not expect LDA to work well, since it is well-known that LDA underpredicts lattice constants; however, we wanted to determine if its performance would improve at high degrees of compression. Additionally, we sought to compare a non-GGA method with the various GGA theories. Figures 1, 4, and 5 illustrate the errors LDA exhibit compared to experimental crystal volumes for RDX, TATB, and PETN respectively. In all cases, LDA volumes are significantly smaller than experimental values. LDA-predicted volumes for RDX at 0.73 and 1.75 GPa are smaller than the experimental values by 11% and 9%, respectively. For the TATB and PETN crystals at zero pressures, the LDA volumes are smaller than experiment by 11% and 15%, respectively. For TATB, the LDA errors decrease with increasing pressure, with a 4% volumetric error at 7 GPa. Volume predictions for PETN exhibit a similar trend, with the error decreasing to approximately 7% by 10 GPa.
In an earlier DFT study, Miao et al.62 investigated the performance of LDA and GGA (PW92) of the I2, Cl2, and Br2 systems under varying degrees of compression. The results indicated that at low pressures, both LDA and GGA did not adequately predict the volumes of these systems. LDA underpredicted the volumes, while GGA overpredicted the volumes, similar to what is seen in this study. However, at higher compression, the LDA and GGA volumes converged to the same values. However, convergence was not observed until ∼80 GPa. Unlike the results of Miao et al.,62 we do not observe the convergence of LDA and GGA volumes at higher pressures, but this study does not explore pressures greater than ∼11 GPa. At the relatively low pressures investigated, LDA is unsuited for use in the prediction of volumes of weakly bound energetic molecular crystals. 3.3. Projected Augmented Planewave PBE. In order to ascertain whether the error seen for PW91 was indicative of GGA methods, PBE calculations were performed on the TATB and PETN systems. The PBE GGA with the PAW pseudopotential and the 545 eV planewave basis set were utilized for this purpose.
Structures of HMX, RDX, CL-20, PETN, and TATB
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2791
Figure 5. Variation of volume and lattice dimensions with pressure for PETN. Experimental values taken from ref 79.
Figure 6. Variation of V/V0 with pressure for PETN where V0 corresponds to (a) the theoretical value for each individual theory; (b) the experimental value.
At zero pressure, the PBE results were more accurate than PW91. The errors in the TATB and PETN zero-pressure volumes were reduced by 4% and 6%, respectively. For TATB, the calculated PW91 and PBE a and b lattice vector lengths were approximately the same, while PBE exhibited greater accuracy in the c lattice vector length at zero pressure. For PETN at zero pressure, however, PBE predictions of all lattice vector lengths were in better agreement with experiment. At increased pressures (>2 GPa), the PBE volumes converge to those produced using PW91. Crystal volumes calculated using the two methods at the elevated pressures differ by at most 4.0 Å.3 Despite changing both GGA methods and pseudopotentials, these results demonstrate that, at least for PW91 and PBE at
545 eV, there is no influence on predicted crystal volumes due to differences in GGA methods or pseudopotentials at pressures greater than 2 GPa. While the GGA results converge relative to one another, however, they still overpredict the experimental values. We next explored the dependence of predicted lattice vector lengths on algorithms used in the crystal optimizations. The aforementioned calculations were all performed using the VASP program package. We used the PWSCF package to perform similar PBE ultrasoft pseudopotential calculations for comparison. While we recognize that there can be no exact comparison between the PAW and Vanderbilt pseudopotentials, in principle there should be a miniscule effect due to pseudopotentials employed for calculations at these pressures. The use of the PWSCF program package not only allows for the comparison of a different combination of theories (namely PBE with ultrasoft pseudopotentials), but also different convergence algorithms, as PWSCF uses damped molecular dynamics to compute optimized geometries. Using PWSCF, we calculated the volume of PETN at zero pressure (where the error is greatest) using the Vanderbilt ultrasoft pseudopotential from VASP and PBE at 545 eV. The PWSCF volume (639.9 Å3) differs by less than a 1% compared to the VASP PBE PAW (645.3 Å3) volume. Therefore we are confident that the choice of pseudopotential has little effect on the computed volumes for calculations at these pressures. It is also heartening to note the convergence method utilized does not seriously affect the final calculated volumes. An additional benefit of using PWSCF was the molecular dynamics algorithm used for determining optimized geometries allows for the ability to explore the nature of the potential energy surface near the minimum. Figures 7a and 7b show the calculated energies, pressures, and crystal volumes versus damped molecular dynamic integration step for PETN using PW91 USP at 545 eV. We have provided insets that illustrate the very small decreases in the internal energies and pressures after only the first hundred integration steps; however, the crystal volume continues to expand. The larger graphs that contain the insets show the energy and pressure decreases at a greatly magnified scale. We have included these figures because the default convergence criteria for such a calculation is 10-4 Rydberg in the energy and 0.5 kbar in pressure. The point at
2792 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Byrd and Rice TABLE 2: Experimental and PBE Calculated PETN Crystalline Volumes at Zero Pressure
Figure 7. Time history of PWSCF damped molecular dynamics calculations for PETN geometry optimizations using PW91 USP at 545 eV. (a) Energy and volume per iteration; (b) the pressure and volume per iteration. Insets in the graphs show the entire ranges of calculated values for the duration of the damped MD geometry optimization.
which the default energy criterion was satisfied was at iteration step 66 (corresponding to a volume of 598 Å3). However, the pressure was not converged at this point; the pressure criterion was met at iteration step 82 (corresponding to a volume of 602 Å3). At these points, however, the PBE 545 eV predicted volumes do not agree with those generated using VASP. Therefore, we tightened the convergence criteria to 10-6 in energy and allowed the optimization to progress. The fluctuations in the pressure profile at later steps indicate that the pressure has not completely met the convergence criterion. At the end of the optimization using the tightened convergence criteria, not only is the pressure converged to under 0.5 kbar, the fluctuations in pressure are also under the convergence limit. Also at this point, the volume is within 1% of the value generated using the VASP program package. Although the volume does not appear to be completely converged at the end of the optimization, we did not continue the calculation, since the almost converged volume is already in significant disagreement with experiment. Continuing the optimization would only produce a volume in even worse agreement with experiment and would provide no additional useful information. The purpose of this exercise, rather, was to illustrate the necessity of precise pressure and energy convergence in this optimization method. Should a user terminate the optimization after the first few hundred iterations, erroneous volumes might be reported. Therefore, a user is strongly cautioned to explore possible sources of errors in predicted volumes that can be attributed to the nature of the potential energy surface near equilibrium or, as previously discussed, the fortuitous choice of the planewave basis set. 3.4. Troullier-Martins PBE. Finally, we investigated PETN with PBE and a Troullier-Martins (TM) very hard pseudopo-
method
volume (Å3)
experiment PBE PAW 545 eV PBE USP 545 eV PBE TM 1890 eV PBE TM 2160 eV PBE TM 2430 eV PBE TM 2700 eV
591.5713 645.2671 639.9162 620.2293 647.3377 646.7022 647.2075
tential with a series of planewave basis sets. The TroullierMartins pseudopotentials were used in order to verify that the non-norm conserving basis sets employed to date (Vanderbilt ultrasoft pseudopotentials and PAW) introduced no additional error. The Troullier-Martins pseudopotential employed was designed specifically for the purposes of this study to be on the opposite end of the spectrum of hardness from that of the Vanderbilt ultrasofts. To this end, the cutoff radii (rc) were chosen to be much smaller than required for the pressures investigated in this study. The rc chosen was set equal to the largest rpeak, which corresponded to where the maximum of the wavefunction (Ψ) occurs. In constructing this pseudopotential, the following rc values were used: carbon (1.2 a0), hydrogen (1.0 a0), nitrogen (1.0 a0), and oxygen (0.9 a0). These very hard radii obviously demand the use of a large number of planewaves, more than are reasonably expected in a standard calculation. However, these calculations should remove any suspicion of error due to the non-norm conserving pseudopotentials. Table 2 illustrates all the PBE theoretical crystal volumes for PETN at zero pressure. Similar to the results seen in the early work of Byrd et al.,61 the volume for the Troullier-Martins pseudopotential at 1890 eV is too small and has not converged. However, at the larger basis sets (2160 eV and larger), the Troullier-Martins results converge both within themselves and to similar solutions as those seen for the Vanderbilt ultrasoft and PAW pseudopotentials. We therefore are confident the choice of pseudopotential does not affect the final converged volumes. 3.5. Structural Analysis. We have performed an extensive structural analysis of the contents of the unit cell for all crystals to determine any pressure-dependent rotational and translational disorder of the molecules within the unit cell. Because we have already established that there are no significant differences in the results predicted by the other GGA methods, we have limited this analysis to results generated using PW91 at 545 eV. Unfortunately, detailed experimental data about atomic positions in the crystals are lacking for pressures above 1 atm. Therefore, when discussing changes in molecular structure (or arrangement of the molecules within the unit cells), we wish to emphasize that we are comparing predicted information about the crystals under compression with those of the experimental structure at room conditions. However, since there are experimental data for the six cell constants (a, b, c, R, β, γ) at the various pressures, we will provide a direct comparison of theoretical values with the corresponding experimental data. Unit cell degrees of freedom that will be compared are the six cell constants (a, b, c, R, β, γ) and the orientational and translational degrees of freedom of each of the molecules in the unit cells. The molecular translational degrees of freedom will be given in fractional coordinates (sx, sy, sz) of the molecular mass centers. Euler angles for rigid-body rotation about inertial axes (θ, φ, ψ) for all molecules in the unit cell are used to describe the molecular orientational degrees of freedom. Symmetry constraints were not imposed in the
Structures of HMX, RDX, CL-20, PETN, and TATB geometry optimizations; therefore, the molecules within the unit cell are not required to be symmetry equivalents. Thus, molecular translational and orientational parameters will be given for each molecule within the unit cell. The contents of the predicted unit cells for the five systems at various pressures are compared with those of the experimental cell at ambient conditions (in lieu of available experimental data). The mass centers of the predicted and experimental cells are located at the origin. The six unit cell constants (a, b, c, R, β, γ) as functions of pressure are given in Supplementary Table 1S in Supporting Information. Fractional values of the mass centers and orientations of all molecules in both experimental and predicted unit cells are given in Supplementary Tables 2S6S in Supporting Information. Pressure-dependent changes in the molecular structure of each of the molecules in the unit cells were determined through atomto-atom comparison with its experimental counterpart (at room conditions) using the method proposed by Kearsley.81 In this procedure, the coordinates of both the experimental and predicted molecules are first shifted such that the molecular centers of mass are located at the origin. The predicted molecule is then superimposed onto the experimental molecule using a rotational matrix that minimizes the sum of the squared distances between corresponding atoms for the two molecules. The rootmean-square (rms) and maximum deviations of predicted atomic displacements from the experimental positions are given in the far right columns in Tables 2S-6S. We emphasize that we are comparing theoretical structures calculated under compression with experimental values measured at 1 atm. For RDX, fractional coordinates of the molecular mass centers are approximately the same as those in the ambient crystal, with differences in the thousandths. Euler angle deviations are less than 3°. The largest rms deviation of the predicted atomic displacements from the experimental positions is 0.064 Å; the maximum atomic displacement is 0.229 Å. The differences in fractional coordinates of the molecular mass centers for β-HMX are similar to those of RDX up to 4.03 GPa, where differences are in the thousandths. For greater pressures, deviations of the fractional coordinates eventually become as large as 0.011. Deviations in the Euler angles are no greater than 2°. Rms deviations of atomic displacements increase with increasing pressure, with the largest rms deviation being 0.106 Å; the maximum atomic displacement is 0.229 Å. For -CL-20, differences in fractional coordinates of the molecular mass centers are also similar to RDX and β-HMX over the pressure range studied (0.35 to 2.5 GPa), with differences in the thousandths. Differences in Euler angles are all smaller than 1° except for a single Euler angle at 2.5 GPa, where the difference is 2.6°. The largest rms deviation of atomic displacements is 0.094 Å, and corresponds to the lowest pressure simulated. The maximum atomic displacement is 0.192 Å. Differences in fractional coordinates of the molecular mass centers of TATB are similar to the RDX, HMX, and CL-20 crystals in that changes are in the thousandths up to ∼4 GPa; at greater pressures the differences are less than 0.015. The differences in Euler angles θ and φ are all less than 1° over the entire pressure range; however, the differences of predicted values of ψ from experimental values range from ∼14° at the lowest pressure up to 56°. At first glance, this large difference between theoretical and experimental values of this Euler angle suggests that the molecules within the unit cell undergo significant rotational reorientation with pressure. However, graphical depictions of the superimposed unit cells (from two
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2793
Figure 8. (a, b) Pictorial views of the 7.02 GPa PW91 (545 eV) unit cell structure of TATB superimposed onto the ambient pressure experimental unit cell from two perspectives. (c) Two-dimensional representation of TATB molecule and inertial axes ωA and ωB projected onto the plane perpendicular to the C3 axis of rotation.
perspectives) at the highest pressure (Figure 8) demonstrate that little rotational reorientation occurs with pressure. The source of this misleading result is due to the alignment of the inertial axes in the molecules. TATB is a nearly planar molecule, as seen in the depiction of the superimposed unit cells in Figure 8b. In this figure, the unit cell calculated at the highest pressure (7.02 GPa) is superimposed onto the experimental unit cell at room conditions. One of the inertial axes lies along the C3 rotational axis that is perpendicular to the plane of the molecule. The remaining two inertial axes (denoted as ωA and ωB) lie within the plane of the molecule. Figure 8c is a two-dimensional representation of the molecule, with the C3 inertial axis coming out of the plane of the paper. In this figure, the arrangement of two inertial axes ωA and ωB at the lowest pressure (and for the experimental system) are depicted in red. The two axes in blue depict the approximate arrangement of ωA and ωB at all other pressures. The geometric difference between the two sets of axes is an approximately 30-42° rotation about the C3 axis. As the Euler angles that
2794 J. Phys. Chem. C, Vol. 111, No. 6, 2007
Byrd and Rice
Figure 9. (a) Pictorial view of the 10.45 GPa PW91 (545 eV) unit cell structure of PETN superimposed onto the ambient pressure experimental unit cell. Frames b and c are the superpositions of each of the two predicted molecules at the same condition superimposed onto its experimental counterpart.
describe the orientation of the body-fixed coordinate system relative to that of the space-fixed system are dependent on the direction of the inertial axes, a change in the directions of these axes will result in a corresponding change in the Euler angles (in this case, ψ). However, the differences in relative orientation of these two inertial axes clearly do not affect the overall molecular orientation within the unit cell, as evident in the graphic depicting the superimposed crystals. The rms in atomic displacements are all ∼0.1 Å and are attributed to the difference in the orientation of the amine groups relative to the ring. In experiment, the amine hydrogens are out of the plane of the molecule, whereas in the theoretical calculations, all atoms lie in approximately within a plane. The analysis of PETN was more challenging than the other systems due to its high-symmetry. PETN is a prolate symmetric top molecule. The unit cell contains two molecules composed of 29 atoms, and the asymmetric unit is composed of eight atoms. Extraction of Euler angles using inertial axes for either of the two complete molecules in the unit cell produced an Euler angle θ ) 90°. For these cases, the remaining Euler angles (ψ and φ) are interrelated, and no unique solution for them exists.82 For purposes of comparison with experiment, however, we are interested in finding a single solution. Therefore, to avoid this complication, we have chosen, instead, to compare Euler angles for the asymmetric unit of PETN. As indicated earlier, the asymmetric unit in the PETN system is composed of eight atoms; there are eight symmetry equivalent fragments within the crystal that can be used as the asymmetric unit; from any of these, positions for all atoms for the two PETN molecules in the unit cell can be generated using the crystal symmetry operations. We will report a comparison of the translational and orientational degrees of freedom for each of these eight fragments. Differences in fractional coordinates of the fragment
mass centers are all less than 0.03. Differences in the Euler angles are no greater than 3.5°. In evaluating the rms and maximum deviation in atomic displacements, we used the full molecules rather than the fragments. Maximum rms of the atomic displacement is 0.117 Å at 10.45 GPa, and the maximum deviation is 0.170 Å. A graphical depiction of the predicted unit cell at 10.45 GPa superimposed on the experimental cell at ambient conditions is given in Figure 9. In this same figure, we have also provided graphical depictions of the superimposed molecules at these same pressures (where the individual molecules were superimposed on the experimental counterpart using the method of Kearsley81 described above). It is clear that the contents of the unit cell, even at this high degree of compression, does not show significant difference in molecular structure, orientation, or arrangement within the unit cell compared to that of the ambient crystal. This small degree of deviation from the low-pressure molecular structure is somewhat surprising, since PETN is an acyclic, branched-chain molecule that one might expect to undergo molecular structural deformation with compression In fact, one experimental study suggests a conformational change to a molecular structure with lower symmetry than that of the highly symmetry S4 state of the ambient crystal and corresponding change in crystalline space group symmetry from P421c to P21212 for pressures higher than 5 GPa.83 In another study,84 we reported several attempts to obtain an optimized unit cell with P21212 symmetry and the molecular structure proposed by Grudzkov et al.83 However, the cell repeatedly converged to the P-421c structure. Since there is conflicting experimental information regarding this phase transition,85-87 none of which agree with one another, we concluded that further experimental measurements were required to resolve discrepancies among
Structures of HMX, RDX, CL-20, PETN, and TATB measured results before we could proceed with further evaluation of other high-pressure phases of PETN. All of the aforementioned results show that there are no significant displacements of the molecular (or fragment) centers of mass nor changes in the orientations of the molecules. Additionally we have utilized the space group symmetry evaluator88 in the Crystal Builder module of the Cerius2 molecular modeling software package89 to determine the space group symmetry for all calculated unit cells. This tool checks to see if the atoms in the system are symmetry copies of each other, according to a specified tolerance. In these calculations, the Cerius2 default tolerance (0.1 Å) was used. For each case, the space group symmetry corresponding to the experimental crystal at ambient conditions was maintained. All of these results indicate that for these calculations, the main structural effect of the compression on the lattices is due almost entirely to a reduction of the intermolecular distances. 4. Conclusion Using the PW91, PBE, and LDA DFT functionals, we have studied five common energetic molecular crystals: HMX, RDX, CL-20, TATB, and PETN over a wide range of pressures. Calculated results near zero pressures are poor for GGA methods, with increasing accuracy as the external pressure increases. GGA methods overpredict volumes by approximately 5% for pressures ranging from ∼4-10 GPa, while LDA underestimates crystal volumes when compared to experiment at all pressures over this range. PW91 and PBE results appear to converge to the same volumes as pressures increase (pressures over 6 GPa), including when utilizing different size planewave basis sets. The lack of agreement with experiment of the predicted lattice vector lengths at lower pressures support the hypothesis that the error is due to a lack of van der Waals forces for non-overlapping densities in current DFT functionals. Inaccurate van der Waals forces cause an underestimation in the attractive forces between the molecules in the crystals and thus an overestimation in the overall lattice vector lengths. This deficiency has unforeseen consequences on all molecular crystal calculations, particularly for those at low pressure, and therefore caution should be employed whenever interpreting results obtained from any solid-state DFT code. We reemphasize to properly describe the electronic structure of the crystals, possible solutions are (1) development of new functionals that include accurate dispersion forces for non- or nearly overlapping densities,90,91 (2) employment of linear response theory to calculate the frequency dependent susceptibility and use it to obtain the dispersion forces, (3) or use perturbation theory for the intermolecular correlation and DFT for the intramolecular correlation, taking care to avoid doublecounting. Acknowledgment. Calculations were performed at the DOD High Performance Computing site at the U.S. Army Research Laboratory, Aberdeen Proving Grounds, MD. The authors thank Dr. Nichols A. Romero and Dr. William D. Mattson for helpful discussions. Supporting Information Available: The six unit cell constants (a, b, c, R, β, γ) as functions of pressure are given in Supplementary Table 1S. Fractional values of the mass centers and orientations of all molecules in both experimental and predicted unit cells are given in Supplementary Tables 2S-6S.
J. Phys. Chem. C, Vol. 111, No. 6, 2007 2795 This material is available free of charge via the Internet at http:// pubs.acs.org. References and Notes (1) Miller, M. S.; Rice, B. M.; Kotlar, A. J.; Cramer, R. J. Clean Products Process. 2000, 2, 37. (2) See, for example, Politzer, P.; Murray, J. S., Ed., Energetic Materials: Part 1. Decomposition, Crystal and Molecular Properties and Part 2. Detonation, Combustion (Theoretical and Computational Chemistry); Elsevier Science: Amsterdam, 2003. (3) Gan, C. K.; Sewell, T. D.; Challacombe, M. Phys. ReV. B 2004, 69, 035116. (4) Brand, H. V. Chem. Phys. Lett. 2006, 418, 428. (5) Brand, H. V. J. Phys. Chem. B 2005, 109, 13668. (6) Reed, E. J.; Joannopoulos, J. D.; Fried, L. E. Phys. ReV. B 2000, 62, 16500. (7) Margetis, D.; Kaxiras, E.; Elstner, M.; Frauenheim, Th.; Manaa, M. R. J. Chem. Phys. 2002, 117, 788. (8) Tuckerman, M. E.; Klein, M. L. Chem. Phys. Lett. 1998, 283, 147. (9) Manaa, M. R.; Reed, E. J.; Fried, L. E. J. Chem. Phys. 2004, 120, 10146. (10) Liu, H.; Zhao, J.; Wei, D.; Gong, Z. J. Chem. Phys. 2006, 124, 124501. (11) Perger, W. F.; Pandey, R.; Blanco, M. A.; Zhao, J. Chem. Phys. Lett. 2004, 388, 175. (12) Perger, W. F.; Zhao, J.; Winey, J. M.; Gupta, Y. M. AIP Conf. Proc. 2006, 845 (Pt. 1, Shock Compression of Condensed Matter-2005, Part 1), 551. (13) Perger, W. F.; Zhao, J.; Winey, J. M.; Gupta, Y. M. Chem. Phys. Lett. 2006, 428, 394. (14) Perger, W. F.; Vutukuri, S.; Dreger, Z. A.; Gupta, Y. M.; Flurchick, K. Chem. Phys. Lett. 2006, 422, 397. (15) Vutukuri, S.; Perger, W. F.; Dreger, Z. A.; Gupta, Y. M. AIP Conf. Proc. 2004, 706 (Pt. 1, Shock Compression of Condensed Matter-2003, Part 1), 275. (16) Zerilli, .F. J.; Kuklja, M. M. AIP Conf. Proc. 2006, 845 (Pt. 1, Shock Compression of Condensed Matter-2005, Part 1), 183. (17) Zerilli, F. J.; Kuklja, M. M. J. Phys. Chem. A 2006, 110, 5173. (18) Zerilli, F. J.; Kuklja, M. M. AIP Conf. Proc. 2004, 706 (Pt. 1, Shock Compression of Condensed Matter-2003, Part 1), 123. (19) Kuklja, M. M.; Zerilli, F. J.; Sushko, P. Mater. Res. Soc. Symp. Proc. 2003, 800 (Synthesis, Characterization and Properties of Energetic/ Reactive Nanomaterials), 211. (20) Kuklja, M. M.; Zerilli, F. J.; Peiris, S. M. J. Chem. Phys. 2003, 118, 11073. (21) Kuklja, M. M. AIP Conf. Proc. 2002, 620 (Shock Compression of Condensed Matter, Pt. 1), 454. (22) Qiu, L.; Xiao, H.-M.; Zhu, W.-H.; Xiao, J.-J.; Zhu, W. J. Phys. Chem. B 2006, 110, 10651. (23) Kunz, A. B.; Kuklja, M. M.; Botcher, T. R.; Russell, T. P. Thermochim. Acta 2002, 384, 279. (24) Kuklja, M. M. J. Phys. Chem. B 2001, 105, 10159. (25) Kuklja, M. M.; Kunz, A. B. J. Appl. Phys. 2001, 89, 4962. (26) Kuklja, M. M.; Aduev, B. P.; Aluker, E. D.; Krasheninin, V. I.; Krechetov, A. G.; Mitrofanov, A. Yu. J. Appl. Phys. 2001, 89, 4156. (27) Kuklja, M. M.; Kunz, A. B. J. Appl. Phys. 2000, 87, 2215. (28) Kuklja, M. M.; Stefanovich, E. V.; Kunz, A. B. J. Chem. Phys. 2000, 112, 3417. (29) Kuklja, M. M.; Kunz, A. B. J. Appl. Phys. 1999, 86, 4428. (30) Kuklja, M. M.; Kunz, A. B. J. Phys. Chem. B 1999, 103, 8427. (31) Kuklja, M. M.; Kunz, A. B. Mater. Res. Soc. Symp. Proc. 1999, 538 (Multiscale Modelling of Materials), 347. (32) Kunz, A. B. Phys. ReV. B 1996, 53, 9733. (33) Kunz, A. B. Mater. Res. Soc. Symp. Proc. 1996, 418 (Decomposition, Combustion, and Detonation Chemistry of Energetic Materials), 287. (34) Zhang, C.; Shu, Y.; Zhao, X.; Wang, X. Sci. Technol. Energ. Mater. 2005, 66, 261. (35) Kuklja, M. M.; Rashkeev, S. N.; Zerilli, F. J. AIP Conf. Proc. 2006, 845 (Pt. 1, Shock Compression of Condensed Matter-2005, Part 1), 535. (36) Kuklja, M. M.; Rashkeev, S. N.; Zerilli, F. J. Appl. Phys. Lett. 2006, 89, 071904/1. (37) Zhurova, E. A.; Tsirelson, V. G.; Zhurov, V.; Stash, A. I.; Pinkerton, A. A. Acta Crystallogr. Sect B 2006, 62, 513. (38) Kuklja, M. M.; Rashkeev, S. N.; Zerilli, F. J. AIP Conf. Proc. 2004, 706 (Pt. 1, Shock Compression of Condensed Matter-2003, Part 1), 363. (39) Peiris, S. M.; Wong, C. P.; Zerilli, F. J. J. Chem. Phys. 2004, 120, 8060. (40) Rashkeev, S. N.; Kuklja, M. M.; Zerilli, F. J. Appl. Phys. Lett. 2003, 82, 1371. (41) Peiris, S. M.; Pangilinan, G. I.; Zerilli, F. J.; Russell, T. P. AIP Conf. Proc. 2002, 620 (Shock Compression of Condensed Matter, Pt. 1), 181.
2796 J. Phys. Chem. C, Vol. 111, No. 6, 2007 (42) Hafner, J.; Wolverton, C.; Ceder, G. MRS Bull. 2006, 31, 659. (43) Oxtoby, D. W. Ann. ReV. Mater. Sci. 2002, 32, 39. (44) Wu, J.; Li, Z. Ann. ReV. Phys. Chem. 2007, in press. (45) Martin, R. M. Electronic Structure: Basic Theory and Practical Methods; Cambridge University Press: Cambridge, UK, 2004. (46) Wu, X.; Vargas, M. C.; Nayak, S.; Lotrich, V.; Scoles, G. J. Chem. Phys. 2001, 115, 8748. (47) Reimers, J. R.; Cai, Z.-L.; Bilic, A.; Hush, N. S. Ann. N. Y. Acad. Sci. 2003, 1006 (Molecular Electronics III), 235. (48) Kamiya, M.; Tsuneda, T.; Hirao, K. J. Chem. Phys. 2002, 117, 6010. (49) Baer, R.; Neuhauser, D. Phys. ReV. Lett. 2005, 94, 043002/1. (50) Sato, T.; Tsuneda, T.; Hirao, K. Mol. Phys. 2005, 103, 1151. (51) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. J. Chem. Phys. 2001, 115, 3540. (52) Williams, R. W.; Malhotra, D. Chem. Phys. 2006, 327, 54. (53) Ortmann, F.; Bechstedt, F.; Schmidt, W. G. Phys. ReV. B 2006, 73, 205101/1. (54) Angyan, J. G.; Gerber, I. C.; Savin, A.; Toulouse, J. Phys. ReV. A 2005, 72, 012510/1. (55) Neumann, M. A.; Perrin, M.-A. J. Phys. Chem. B 2005, 109, 15531. (56) Kleis, J.; Schroeder, E. J. Chem. Phys. 2005, 122, 164902/1. (57) Grimme, S. J. Comput. Chem. 2004, 25, 1463. (58) Kamiya, M.; Tsuneda, T.; Hirao, K. J. Chem. Phys. 2002, 117, 6010. (59) Wu, Q.; Yang, W. J. Chem. Phys. 2002, 116, 515. (60) Sato, T. I.; Tsuneda, T.; Hirao, K. J. Chem. Phys. 2005, 123, 104307. (61) Byrd, E. F. C.; Scuseria, G. E.; Chabalowski, C. F. J. Phys. Chem. B, 2004, 108, 13100. (62) Miao, M. S.; Doren, V. E.; Martins, J. L. Phys ReV B 2003, 68, 094106. (63) Perdew, J. P. In Electronic Structures of Solids ’91; Ziesche, P. Eschrig, H., Eds.; Akademie-Verlag: Berlin, 1991. (64) Perdew, J. P.; Berke, K.; Ernzerhof, M. Phys. ReV. Lett. 1996, 77, 3865. (65) Kresse, G.; Furthmuller, J. Vienna Ab-initio Simulation Package (VASP): The Guide; VASP Group, Institut fur Materialphysik, Universitat Wien: Sensengasse 8, A-1130 Wien, Vienna, Austria, 2003. (66) Ceperley, D. M.; Alder, B. J. Phys. ReV. Lett. 1980, 45, 566. (67) Vanderbilt, D. Phys ReV B. 1990, 41, 7892. (68) Blo¨chl, P. E. Phys. ReV. B. 1994, 50, 17953. (69) Kresse, G.; Joubert, J. Phys. ReV. B. 1999, 59, 1758.
Byrd and Rice (70) Baroni, S.; Dal Corso, A.; de Gironcoli, S.; Giannozzi, P.; http:// www.pwscf.org. (71) Choi, C. S.; Prince, E. Acta Crystallogr. Sect. B 1972, 28, 2857. (72) Nielson, A. T.; Chafin, A. P.; Christian, S. L.; Moore, D. W.; Nadler, M. P.; Nissan, R. A.;Vanderah, D. J.; Gilardi, R. D.; George, C. F.; Flippen-Anderson, J. L. Tetrahedron 1998, 54, 11793. (73) Choi, C. S.; Boutin, H. P. Acta Crystallogr. Sect. B 1970, 26, 123. (74) Cady, H. H.; Larson, A. C. Acta Crystallogr. 1975, 31, 1864. (75) Cady, H. H.; Larson, A. C. Acta Crystallogr. 1965, 18, 485. (76) Methfessel, M.; Paxton, A. T. Phys. ReV. B 1989, 40, 3616. (77) Olinger, B. W.; Cady, H. H. The Hydrostatic Compression of ExplosiVes and Detonation Products to 10 GPa (100 kbar) and their Calculated Shock Compression: Results for PETN, TATB, CO2 and H2O; Sixth Symposium on Detonation, San Diego, CA, Aug 24-27, 1976; Los Alamos LA-UR-76-1174. (78) Olinger, B. W.; Roof, B.; Cady, H. H. Symposium international Sur Le Comportement Des Milieus Denses Sous Hautes Pressions Dynamiques; Commissariat a l’Energie Atomique Cnetre d’Etudes de Vajours: Paris, France, 1978; p 3. (79) Olinger, B. W.; Halleck, P. M.; Cady, H. H. J. Chem. Phys. 1975, 62, 4480. (80) Pinkerton, A. Private communication, 1999. (81) Kearsley, S. K. Acta Crystallogr. 1989, A45, 208. (82) Mukundan, R. Proceedings of the 7th Asian Technology Conference in Mathematics; ATCM, Inc.: Malaysia, 2002; p 97. (83) Gruzdkov, Y. A.; Gupta, Y. M. J. Phys. Chem. A 2001, 105, 6197. (84) Rice, B. M.; Byrd, E. F. C. J. Mater. Res. 2006, 21, 2444. (85) Lipinska-Kalita, K. E., Pravica, M. G.; Nicol, M. J. Phys. Chem. B 2005, 109, 19223. (86) Gruzdkov, Y. A., Dreger, Z. A., and Gupta, Y. M. J. Phys. Chem. A 2004, 108, 6216. (87) Ciezak, J. U. S. Army Research Laboratory, 2006, private communication. (88) Cerius2 Builder. http://www.chem.cmu.edu/courses/09-560/docs/ msi/builder/ 3_Crystal.html#10784 (accessed Nov 2006). (89) Accelrys Cerius2. http://www.accelrys.com/products/cerius2/(accessed Nov 2006). (90) Rydberg, H.; Dion, M.; Jacobson, N.; Schro¨der, E.; Hyldgaard, P.; Simak, S. I.; Langreth, D. C.; Lundqvist, B. I. Phys. ReV. Lett. 2003, 91, 126402. (91) Dion, M.; Rydberg, H.; Schro¨der, E.; Langreth, D. C.; Lundqvist, B. I. Phys. ReV. Lett. 2004, 92, 246401.
