(TATB) Using Symmetry Adapted Perturbation ... - ACS Publications

Article pubs.acs.org/JPCA

Intermolecular Forces and Molecular Dynamics Simulation of 1,3,5Triamino-2,4,6-trinitrobenzene (TATB) Using Symmetry Adapted Perturbation Theory DeCarlos E. Taylor* U.S. Army Research Laboratory, Aberdeen Proving Ground, Maryland 21005, United States ABSTRACT: The dimer potential energy surface (PES) of 1,3,5-triamino2,4,6-trinitrobenzene (TATB) has been explored using symmetry adapted perturbation theory based on a Kohn−Sham density functional theory description of the monomers [SAPT(DFT)]. An intermolecular potential energy function was parametrized using a grid of 880 ab initio SAPT(DFT) dimer interaction energies, and the function was used to identify stationary points on the SAPT(DFT) dimer PES. It is shown that there exists a variety of minima with a range of bonding configurations and ab initio analyses of the interaction energy components, along with radial cross sections of the PES near each minimum, are presented. Results of isothermal-isostress molecular dynamics simulations are reported, and the simulated structure, thermal expansion, sublimation enthalpy, and bulk modulus of the TATB crystal, based on the SAPT(DFT) interaction potential, are in good agreement with experiment.

1.. INTRODUCTION 1,3,5-Triamino-2,4,6-trinitrobenzene1 (TATB) is a crystalline energetic material with extraordinary stability. The unit cell (Figure 1) contains two molecules in the P-1 space group and the crystal consists of graphitic-like sheets with extensive hydrogen bonding within each layer and weaker van der Waals interactions between layers. Experimental results2 suggest an unusually high insensitivity to impact and shock loading and as a result, there have been numerous investigations on the origin of this stability with phenomena such as intramolecular conjugation3 and intermolecular hydrogen bonding4 being suggested as stabilizing mechanisms within the crystal. Experimental studies5,6 have shown significant anisotropy in the thermal and mechanical response of TATB. The anisotropy has been attributed to differences in intermolecular force between coplanar and interlayer molecules and as a result, the intermolecular interactions within TATB have been studied extensively both experimentally and theoretically. Davidson et al.4 measured Raman spectra of compressed TATB and concluded that a spectral shift at 28 GPa was due to changes in the intralayer hydrogen bond network due to reorientation of the nitro groups. Manaa and Fried,7 using density functional theory (DFT), computed the TATB crystal structure up to 250 GPa and found that as the pressure was increased, the strength of the hydrogen bond network increased accordingly. At a pressure of 30 GPa, the intra- and intermolecular hydrogen bonds, using a distance criterion, were essentially equivalent. Ojeda and Cagin8 showed a similar increase in hydrogen bond character with pressure and found that the total number of hydrogen bonds between layers increased from 0 to 8, per monomer, at a pressure of 1.6 GPa. Valenzano et al.9 computed This article not subject to U.S. Copyright. Published XXXX by the American Chemical Society

second-order elastic constants of TATB, using DFT, with thermal effects incorporated by fixing the volume of the computational unit cell to that obtained experimentally at 295K. The Cij values reported in that work indicate that the stiffness of TATB in the a/b direction (Figure 1) exceeds that along c by a factor of 4. Theoretical studies such as those presented above provide information on the magnitude and relative strength of the intermolecular interactions as a function of geometry however additional insight can be obtained through the application of energy decomposition techniques, such as symmetry adapted perturbation theory (SAPT),10 which enable determination of the magnitude of each distinct contribution (electrostatics, induction, dispersion, etc.) to the total interaction energy. Taylor11 applied SAPT to dimers extracted from a series of energetic molecular crystals and showed that the ordering of interaction energies for compounds of similar structure (nitroaromatics, nitroanilines, etc.) correlated with experimentally measured impact sensitivities. The SAPT interaction energy of the TATB dimer extracted from the experimental unit cell was dominated by the dispersion contribution. An analysis of TATB interaction energies based on SAPT was also conducted by Song et al.12 Roszak et al.13 used DFT and second order many body perturbation theory with a Moller− Plesset Hamiltonian partitioning (MP2) to identify stable gas phase configurations of TATB dimers and applied a variationalperturbational energy decomposition method to separate the Received: January 16, 2013 Revised: April 4, 2013

A

dx.doi.org/10.1021/jp4005289 | J. Phys. Chem. A XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry A

Article

Figure 1. TATB monomer, unit cell, and arrangement of layers in crystal.

ologies have been used extensively for analysis of TATB dimers.13,16 However, given that the interlayer interaction within TATB has been shown to be dominated by the dispersion contribution,11,13 the application of interaction energy analyses based on DFT methods without van der Waals correction for dispersion is uncertain. In the results presented by Roszak et al.,13 the dispersion contribution to the interaction energy for a π-stacked TATB configuration obtained with B3LYP was only 30% of that given by MP2. Although accurate interaction energies, including those with a significant contribution from dispersion, can be obtained using explicitly correlated wave function based methods such as coupled cluster theory,18 the computational expense of such approaches renders them difficult to apply to large systems such as TATB. As a result, dimer potential energy surfaces have been studied using a combined DFT and symmetry adapted perturbation theory approach known as SAPT(DFT).19 In SAPT(DFT), the intramonomer correlation perturbation is included implicitly via the exchange-correlation density functional, therefore the SAPT double perturbation series for the intra- and intermonomer perturbations can be reduced to a single series for the intermolecular interaction only which corresponds to truncation of the intramonomer perturbation series at zeroth order. In SAPT(DFT), due to the incorrect behavior of exchange-correlation potentials at large r, asymptotically corrected20 potentials are necessary and the Fermi−Amaldi21 long-range form with the splicing scheme of Tozer and Handy22 is applied. Accurate dispersion energies, of particular importance for the interlayer interactions in TATB, are obtained via application of the generalized Casimir-Polder expression with frequency dependent density susceptibilities computed within the coupled Kohn−Sham approach.23 Using resolution of the identity techniques24(density fitting), SAPT-

total interaction energy into its constituent parts. For planar dimer complexes with orientations similar to those found within layers of the crystal, it was found that the dominant contribution resulted from the electrostatic term (which they attributed to intermolecular hydrogen bonding) and π-stacked complexes, similar to those found between layers in the crystal, were dispersion bound. The energy decompositions presented by Roszak et al.13 were applied to dimer orientations similar to those present in the experimental crystal structure. However, Zhang et al.14 presented potential energy surface (PES) cross sections using a combined DFT and Dreiding15 forcefield approach indicating the presence of a diversity of additional minima including planar, π-stacked nitro to nitro, π-stacked nitro to amino, and T-shaped dimer geometries. In that work, the minima were identified by selectively configuring TATB monomers relative to one another, based on chemical intuition, and performing static scans of the PES by variation of the intermonomer distance. At each separation, energies were computed using the Dreiding forcefield and the Dreiding interaction energies were compared to another parametrized forcefield developed by Gee et al.16 However, it has been shown17 that the Gee forcefield underestimates the experimental sublimation enthalpy significantly therefore this forcefield may not be a good metric for validation of the dimer interaction energies obtained using the Dreiding forcefield. Given the variety of minima present on the TATB dimer PES and the importance of the intermolecular interactions for this system, a more rigorous analysis of the relevant minima on the dimer PES, based on first principles quantum mechanics, is warranted. Due to the large size of the TATB dimer (48 atoms) and large basis sets required to obtain accurate interaction energies, computationally affordable MP2 and DFT based methodB



Article

(DFT) has a very economical N5 scaling which renders it applicable to large systems and it has been used to obtain accurate dimer potential energy surfaces for other energetic crystals such as cyclotrimethylene trinitramine25 (RDX) and 1,1-diamino-2,2-dinitroethylene26 (FOX-7). In this work, an exploration of the TATB dimer PES using SAPT(DFT) is presented. Interaction energies for 880 TATB dimer configurations were computed and fitted to a classical potential energy function (forcefield). Using this forcefield, minima on the dimer PES were located and it is shown that there exists a diversity of minima including planar, π-stacked, and T-shaped orientations, some of which are similar to arrangements found in the crystal. Further, it is shown that the MP2 and B3LYP π-stacked configurations considered by Roszak et al.13 constitute local minima on the SAPT(DFT) surface and a similar π-stacked configuration of lower energy exists . For selected minima, radial cross sections of the PES, computed quantum mechanically, are presented and, in addition to the ab initio analyses, molecular dynamics (MD) simulations for a range of temperatures and pressures are reported. The SAPT(DFT) fitted forcefield yields MD results in good agreement with experiment for several properties including the crystal structure, bulk modulus, and sublimation enthalpy, all of which were obtained without empirical fitting of experimental results. The well-known anisotropic thermal expansion of TATB is also demonstrated and compared to the recent (2010) experiments of Sun et al.6

dispersion correction reduced the error in the TATB crystal density from 6.3% to 4.1%. However, it has been determined that accurate prediction of performance properties (detonation velocity, detonation pressure) of an energetic material requires a theoretically derived crystal density that is within 3% of the experimental value.33 Whereas the Reax forcefield exceeds that metric of error for TATB, as will be shown, the SAPT(DFT) forcefield developed herein yields a crystal density within 0.6% of the experimental value. The forcefield in this work was parametrized from first principles by fitting to a grid of TATB dimer interaction energies computed using the density fitting version of SAPT(DFT) as implemented in the SAPT200834 software program. The geometry of the TATB monomer used in the SAPT(DFT) calculations was taken from the experimental1 unit cell (the only experimental information used in this work). This geometry was chosen since, in addition to exploration of the dimer PES, condensed phase molecular dynamics simulations are an end goal of this work. The SAPT(DFT) energies (and forcefield parameters determined using those energies) are a function of the intramolecular geometry and may not be transferable to other configurations. Since crystal field effects may alter the geometry of the TATB monomer (relative to the gas phase), in order to have a forcefield with parameters applicable to condensed phase configurations near the ambient state, the experimental geometry was chosen for the gas phase calculations. All calculations were done using an aug-cc-pVDZ basis in the MC+BS35 approach. In the MC+BS approach, the basis set of each monomer is augmented via additional s/p (for first row atoms) basis functions taken from the aug-cc-pVDZ basis set. These additional functions are centered on the nuclei of the adjacent monomer and this basis yields monomer energies comparable to the full dimer basis with only a modest increase in computational cost. The basis was also supplemented by a set of 3s(α=0.9,0.3,0.1) 3p(α=0.9,0.3,0.1) 2d(α=0.6,0.2) 2f(α=0.6,0.2) “midbond functions”. The total number of contracted gaussians in the MC+BS basis was 756, and the location of the midbond functions, as done in previous work,25,26 was determined using a weighted average of atom−atom midpoints

2. FORCE FIELD DEVELOPMENT 2.1. Ab Initio Calculations. Current SAPT(DFT) implementations do not have analytic energy gradients therefore an efficient location of local minima is not possible at the ab initio level. Gee et al.16 developed an analytic TATB forcefield by fitting to MP2 interaction energies and the model yielded a crystal structure in good agreement with experiment. However, the sublimation enthalpy, which is indicative of the strength of the intermolecular interactions within the crystal, underestimated the experimental values of Rosen,27 as well as those of Garza,28 by 23 and 28%, respectively. Therefore, the intermolecular energetics provided by that forcefield may not accurately characterize all of the relevant minima on the dimer surface. The TATB forcefields from Bedrov et al.29 and Filippini and Gavezzotti30 yield accurate structures and sublimation enthalpies. However, these models were parametrized, in part, using experimental reference data with the equilibrium unit cell parameters and sublimation enthalpy included as part of the fit. A goal of the present work is the development of a predictive forcefield, purely from quantum mechanically derived reference data, without regard to fitting of experimentally measured material properties such as crystal structure, sublimation enthalpy, elastic moduli etc. A predictive capability is critical since the fitting of experimental data is not possible for materials where no such data is available. This same strategy has been applied in development of the Reax31 forcefield, where only quantum mechanically derived reference data is used for parametrization. Liu32 et al. augmented the original Reax functional form with a dispersion correction suitable for application to several energetic molecular crystals, including TATB. In their work, low temperature crystal densities and heats of sublimation were used to fit parameters of a long-range correction which was included to compensate for the inadequate treatment of dispersion interactions using the original Reax functional form. The inclusion of the

rmb =

1 2

∑ ∑ ωab(ra + rb) a∈A b∈B

(1)

where the weights ωab are given by

ωab =

rab‐6 ∑ab rab−6

(2)

and indices (a,b) run over atoms of monomer A and B, respectively. The density fitting basis, consistent with the augcc-pVDZ basis used on each monomer, was taken from 37 Weigend et al.36 and the δEHF to include terms int correction, higher than second order, was also used. The Kohn−Sham orbitals (PBE0 functional38) were obtained using a modified version of the DALTON39 program system which includes the Fermi-Amaldi21 asymptotic correction and the splicing scheme of Tozer and Handy22 with switching values of 3 and 4 (Bragg− Slater radii). The ionization potential (0.3257 au) required for the asymptotic correction was determined using Kohn−Sham energies computed for the N and N − 1 electron states and all calculations were run on 128 cores of an IBM P6 (4.7 GHz) distributed memory computer at the Navy Defense Shared Resource Center. C



Article

2.2. Parameterization. The functional form of the forcefield is the sum of the terms Vab(r ) =

qaqb r

⎛ r ⎞ C ab + A ab exp⎜ − ⎟ − 66 ⎝ Bab ⎠ r

Table 1. Monomer Coordinates (Å), Charges, and Pair Parameters

(3)

where a and b enumerate atoms of monomer A and B, respectively. The atomic charges were determined by fitting to the PBE electrostatic potential of a single TATB molecule using the CHELPG40 procedure as implemented in the Gaussian41 quantum chemistry program. In previous work,25,26 the charges were determined by fitting to a large (more than 59 000 data points in some cases) grid of electrostatic energies computed using asymptotic multipole expansions based on the quantum mechanical density. However, for the FOX-7 energetic, a comparison of the long-range electrostatic energies, computed using CHELPG derived charges, already showed good agreement with the electrostatic component of the ab initio SAPT(DFT) interaction energies for the same asymptotic points. The CHELPG fitted charges determined for TATB yield a dipole moment value of 1.21D which compares well with the ab initio value of 1.18D. Using this fixed set of charges, the remaining pair parameters Aab, Bab, and Cab were fitted to SAPT(DFT) dimer interaction energies sampling the repulsive wall, potential well, and dissociative tail. In order to improve the level of agreement between the fitted energies and SAPT(DFT) reference values, the three carbons bonded to nitro groups (see Figure 1) were assigned different pair parameters from carbons bonded to amine groups. A similar distinction was made between nitrogen atoms within −NO2 and −NH2 groups. As a result, there were 63 parameters total (not including the 24 atomic charges). The ab initio grid was generated using randomly oriented configurations between TATB monomers and for each angular orientation, a radial grid was computed by variation of the center of mass separation from the repulsive region, out to dissociation, in steps of 0.5 Å with center of mass separations ranging from ∼3.0 to ∼16 Å (depending on the orientation). Points were continually computed in this way until there was good agreement between ab initio and forcefield interaction energies for a set of points not included in the fit (indicating convergence of the parametrization). Parameter saturation required 880 SAPT(DFT) points total. Using the PIKAIA42 genetic algorithm (GA), the parameters were optimized by minimizing the root-mean-square deviation between the forcefield and the SAPT(DFT) interaction energy reference data set. Although each contribution to the SAPT(DFT) total interaction energy could have been fit separately (as done in previous work), it was found that with the relatively small number of parameters in the current potential, a better fit was obtained using total energies with all parameters varied simultaneously. During the GA run, each data point was assigned a weight given by e[(1/3)(10.0−Ei)] which weighs the bound configurations most heavily and exponentially decreases to one as the interaction energy approaches 10 kcal/mol. The values of 1/3 and 10 kcal/mol in the weighting factor were chosen empirically to maximize the quality of the fit. The root-mean-square deviation for all points with an energy less than 30 kcal/mol is 0.75 kcal/mol and for all bound configurations the deviation is 0.31 kcal/mol. Cartesian coordinates of the monomer, fitted atomic charges, and all additional potential parameters are given in Table 1.

atom

atomic charge

X

Y

Z

C-nitro C-amine C-nitro C-amine C-nitro C-amine H H H H H H N-nitro N-amine N-nitro N-amine N-nitro N-amine O O O O O O

−0.38383 0.48795 −0.41075 0.53015 −0.37935 0.43779 0.45446 0.43947 0.48086 0.42375 0.42803 0.37691 0.74437 −0.81518 0.76648 −0.86762 0.73146 −0.71383 −0.43213 −0.47212 −0.46401 −0.45584 −0.45294 −0.45406

1.42707 0.72099 −0.72649 −1.4513 −0.69328 0.74121 2.22418 0.93593 −3.19544 −3.29101 0.97153 2.34591 2.84158 1.35994 −1.4479 −2.76959 −1.39978 1.40956 3.48782 3.47323 −0.85059 −2.68622 −2.64218 −0.78975 pair parameters

−0.01928 −1.26886 −1.2321 0.00823 1.23939 1.25798 −2.37943 −3.13335 −0.59144 0.94927 3.13412 2.32644 −0.02765 −2.40995 −2.45158 0.03362 2.47191 2.38575 1.02344 −1.09742 −3.54316 −2.47153 2.50434 3.57437

0.04210 −0.00877 −0.00877 0.00630 −0.03703 0.011960 −0.07596 0.06220 0.02452 0.02452 −0.03828 −0.16388 0.12186 −0.07031 −0.02133 0.04838 −0.0948 0.01007 0.20036 0.08794 −0.07785 0.0151 −0.12369 −0.09041

a

b

Aab (kcal/mol)

Bab (Å)

Cab (Å6 kcal/mol)

C-nitro C-nitro C-nitro C-nitro C-nitro C-nitro C-amine C-amine C-amine C-amine C-amine H H H H O O O N-nitro N-nitro N-amine

C-nitro C-amine H O N-nitro N-amine C-amine H O N-nitro N-amine H O N-nitro N-amine O N-nitro N-amine N-nitro N-amine N-amine

63659.99859 217198.5658 67520.01563 60230.69865 219341.4172 5173.14212 76572.06124 43345.16654 36710.52593 56611.45465 96237.17353 68609.90494 81993.70178 37599.94741 121564.6467 47790.66867 8924.37424 13948.30802 72136.29346 76381.72781 34106.92661

0.29036 0.22791 0.07126 0.27173 0.23291 0.40641 0.28458 0.24712 0.25792 0.15626 0.24854 0.20192 0.1978 0.22782 0.07956 0.26222 0.31077 0.3263 0.05537 0.23922 0.30933

337.15867 535.83737 63.3226 209.38151 117.66595 1038.67756 184.23356 569.26353 134.34444 550.44004 160.32828 94.67707 228.10028 274.52737 116.30439 297.70704 351.84632 133.03496 23.25206 278.01648 405.14849

3. RESULTS 3.1. Analysis of Minima on Dimer Potential Energy Surface. Starting from random monomer orientations, local minima on the dimer potential energy surface were located using the L-BFGS43 optimization algorithm where the internal geometry of each monomer was kept fixed, consistent with the ab initio calculations, and the radial separation and angular orientation between the rigid monomers were optimized. Clearly one cannot locate all minima on the surface nor can it be guaranteed that all minima on the fitted surface correspond to minima on the ab initio surface. However, for each D



Article

Figure 2. Minima on PES. Left and right columns are different perspectives of the same configuration.

crystal that range from ∼2.2 to ∼2.4 Å. The dimer B center of mass separation, 4.73 Å, is shorter than that for the similar interaction occurring in the experimental unit cell (5.09 Å). Further, the monomers in dimer B have a 3.5 Å translational offset (see side perspective in Figure 2) compared to a 4 Å offset in the unit cell. Therefore, the gas phase structure favors a larger overlap between monomers than what is observed in the crystal. Roszak et al.13 report a similar gas phase structure using B3LYP; however, the intermonomer overlap at the B3LYP level is less than that obtained using SAPT(DFT) or MP2.13 As discussed below, for dimer B, dispersion constitutes the largest contribution to the total interaction energy. The B3LYP functional, without van der Waals correction, likely underestimates this attractive interaction which results in a larger intermonomer separation. The increase in overlap observed using SAPT(DFT) validates the findings reported by Roszak et al. using MP2.

configuration discussed in the following, radial cross sections were computed ab initio to verify that each is indeed a true stationary point on the SAPT(DFT) surface. Images of selected local minima are shown in Figures 2 and 3, and for each configuration, two perspectives are given for clarity of representation. As shown, there is a diversity of structures including planar (dimer A), π-stacked (dimers B and C), T-shaped (dimers E and F), and other bonding configurations with center of mass separations ranging in value from 3.34 to 9.07 Å. Of particular interest are dimers A and B whose geometries are very similar to the orientations found in the experimental crystal with dimer A resembling interactions within layers and dimer B resembling those between layers. The center of mass separation in dimer A (9.04 Å) is essentially equivalent to that found in the crystal (9.02 Å), and the amine-hydrogen to nitro-oxygen distance between the two monomers is 2.42 Å which compares with values found between coplanar monomers in the experimental E



Article

Figure 3. Minima on PES. Left and right columns are different perspectives of the same configuration.

Table 2. Interaction Energies (kcal/mol) Using ab Initio SAPT(DFT) and Fitted Forcefielda dimer

R (Å)

electrostatic

exchange

induction

dispersion

ETotal

forcefield

error

A B C D E F G H

9.04 4.73 3.34 7.36 5.89 6.52 9.07 5.89

−5.395 −4.892 −5.428 −3.867 0.439 −4.634 −4.065 −4.467

7.290 16.616 24.117 6.348 7.161 9.361 5.508 11.302

−2.385 −5.301 −7.876 −2.161 −2.681 −3.083 −1.84 −3.711

−4.893 −16.468 −22.465 −5.117 −8.523 −7.569 −4.059 −9.215

−5.384 −10.046 −11.652 −4.798 −3.603 −5.926 −4.465 −6.091

−4.834 −9.798 −11.521 −4.885 −4.090 −5.135 −4.113 −5.893

0.550 0.248 0.131 −0.087 −0.487 0.791 0.352 0.198

Exchange value is sum of first order exchange, and second order exchange-induction/exchange-dispersion energies. Induction energy includes δEHF int contribution. a

Dimer C is equivalent to the π-stacked structure discussed in the work by Zhang14 and, as discussed below, corresponds to the SAPT(DFT) minimum energy configuration among all of

the geometries analyzed in this study. The center of mass separation, 3.34 Å, compares with Zhang’s value of ∼3.4 Å, and in this configuration, the amine nitrogens of the upper F



Article

however, in the former, there is a significant out of plane tilt to the upper monomer compared to the T-shaped configurations. Dimers D and G both represent hydrogen bonded structures where the interaction occurs via amine-hydrogens to nitro oxygens. Both hydrogens of dimer G have a ∼2.50 bond length to the same nitro-oxygen of the adjacent monomer. The hydrogen to oxygen length in dimer D is shorter than that in dimer G with a value of 2.38 Å. 3.2. Interaction Energies. SAPT(DFT) interaction energy components (along with total interaction energies computed using the forcefield) for each structure in Figures 2 and 3 are presented in Table 2. The largest difference between the ab initio and forcefield total interaction energies for these minima is 0.791 kcal/mol (dimer F). Dimer C, a π-stacked structure, represents the global minimum for all points considered in this work with a value of −11.652 kcal/mol. The corresponding interaction energy using the DREIDING forcefield (−12.1 kcal/mol, reported by Zhang14) compares well to the SAPT(DFT) result. Dimer B, another π-stacked structure whose orientation is similar to that found between layers in the experimental crystal, has an energy of −10.046 kcal/mol and is +1.606 kcal/mol less stable than the eclipsed π-stacked structure dimer C. Therefore the eclipsed conformation is more energetically favorable in the gas phase yet the staggered conformation of dimer B is adopted in the crystal. As suggested by Roszak,13 this may be a result of many body interactions or “three-way ring contacts” present in the crystal. The contribution from dispersion in these two minima is a factor of 1.8 (or more) times larger than that found in the remaining configurations. In these π-stacked orientations, there is maximal interaction of the individual electron distributions due to π overlap resulting in an increase in the dispersion energy contribution compared to other configurations. The translational shift in dimer B reduces this overlap and results in a smaller dispersion contribution (−16.468 kcal/mol) than in dimer C (−22.465 kcal/mol). This loss of dispersion interaction may be compensated for in the experimental crystal

Figure 4. Electrostatic potential of an unperturbed TATB monomer. Red mesh indicates region of negative potential and blue mesh signifies positive character.

monomer are oriented directly above the nitro groups of the lower monomer. The shortest O−H bond length between the two monomers is 3.31 Å. Dimers E and F are T-shaped configurations with the upper monomer of dimer E interacting through an amine group and that of dimer F bonding through a nitro group. The center of mass separation of dimer E, 5.89 Å, is significantly shorter than that of dimer F, 6.52 Å; however, (as discussed later) the interaction energy of dimer F is stronger despite the larger intermonomer separation. Dimer H is similar to dimer F in that both contain a nitro group oriented over the center of the carbon ring of the adjacent monomer;

Figure 5. PES cross section of dimer A. Exchange value is the sum of the first order exchange and the second order exchange-induction and exchange-dispersion energies. Induction energy includes δEHF int contribution. G



Article

Figure 6. PES cross section of dimer B. Exchange value is the sum of the first order exchange and the second order exchange-induction and exchange-dispersion energies. Induction energy includes δEHF int contribution.

Figure 7. PES cross section of dimer C. Exchange value is the sum of the first order exchange and the second order exchange-induction and exchange-dispersion energies. Induction energy includes δEHF int contribution.

by an increased interaction with additional coplanar neighbors.13 Dimer A is representative of coplanar dimer configurations in the crystal structure where hydrogen bonding (an electrostatic interaction) between monomers is prevalent. In all configurations listed in Table 2, the dispersion energy exceeds that of the electrostatic term. However, for dimer A, the electrostatic component (−5.395 kcal/mol) has the largest magnitude (when neglecting the exchange repulsion contribution resulting from antisymmetrization of the wave function) indicating that the nitro-oxygen to amine-hydrogen bond interaction between monomers in this configuration leads to stabilization of the structure.

The T-shaped dimers E and F, though similar in structure, have interaction energies differing by ∼1 kcal/mol with the nitro group to carbon-ring interaction of dimer F being preferred over the amine group to carbon-ring bonding configuration in dimer E. Interestingly, the electrostatic contribution for dimer E is slightly repulsive at the minimum. An electrostatic potential contour of an unperturbed TATB monomer is shown in Figure 4. The electrostatic potential is positive in regions covering the 6-member carbon ring and amine hydrogens with negative values near the more electronegative nitro groups. This results in repulsive and attractive electrostatic interactions between the amine (dimer E) and nitro groups (dimer F) when oriented over the positive potential associated with the carbon ring and results in the H



Article

Figure 8. PES cross section of dimer F. Exchange value is the sum of the first order exchange and the second order exchange-induction and exchange-dispersion energies. Induction energy includes δEHF int contribution.

Table 3. Comparison of Experimental and Predicted Crystallographic Parameters property

MDa

experimentb

deviation (%)

a (Å) b (Å) c (Å) α (°) β (°) γ (°) density (g/cm3) ΔHsub (kcal/mol)

9.053 ± 0.051 9.071 ± 0.043 6.818 ± 0.042 111.480 ± 0.528 89.310 ± 0.494 120.15 ± 0.998 1.950 38.20

9.010 9.028 6.812 108.58 91.82 119.97 1.938 40.227 43.128

0.48 0.48 0.09 2.6 2.7 0.15 0.62 4.9 11.4

Figure 9. Simulated TATB unit cell configuration superimposed onto the experimental structure.

Uncertainties computed using the root-mean-square fluctuations of the cell vector components averaged over the final 100 ps of the MD simulation. bReference 1. a

still in very good agreement with the ab initio surface in the repulsive and dissociative regions (see Figure 8). For hydrogen bound dimer A (Figure 5), the attractive interaction is dominated by the electrostatic energy at all values of r; however, for the remaining configurations, the dispersion energy represents the largest attractive contribution for each intermonomer separation. Interestingly, dimer C (Figure 7), which represents the lowest energy configuration obtained in this work, is dipole−dipole repulsive at large intermonomer separation as evidenced by its positive electrostatic energy at r values beyond ∼4 Å. The interaction in the planar dimer A configuration (Figure 5) is much longer ranged than that observed in dimer B (Figure 6) with interaction energies going to zero at values of 12.5 and 8.5 Å, respectively. This is

repulsive electrostatic contribution observed in dimer E at the minimum. In addition to values at the stationary point, SAPT(DFT) energies as a function of center-of-mass separation (at fixed angular orientation) for each configuration have been computed. Plots of radial cross sections for selected configurations are shown in Figures 5−8. In each plot, the forcefield surface is in good agreement with the ab initio surface and stationary points on the fitted surface do indeed correspond to ab initio minima. Dimer F, which contains the largest error (0.791 kcal/mol) at the minimum between the ab initio and forcefield interaction energies reported in Table 2, is

Table 4. Experimental and SAPT(DFT) Center of Mass Orientational Parameters for the Two Symmetry-Equivalent Molecules in the TATB Unit Cell molecule 1

molecule 2

fractional

expt.

SAPT(DFT)

diff

expt.

SAPT(DFT)

diff

sx sy sz

0.37476 0.16531 0.2501

0.37419 0.1643 0.24886

−0.00056 −0.00101 −0.00123

0.62524 0.83469 0.7499

0.62614 0.83671 0.75299

0.00089 0.00203 0.00309

I



Article

Figure 10. Unit cell dimensions as a function of temperature from MD and experiment.

error than that obtained using the polarizable forcefield of Bedrov29 (0.5%); however, in that work, the enthalpy of sublimation was directly parametrized using the NO2−NO2 repulsion and dispersion parameters. The SAPT(DFT) value is a significant improvement over that obtained using the Gee16 potential (31.1 kcal/mol) based on MP2 energetics which is in error by 23% with respect to the experimental value reported by Rosen.27 The thermal expansion of TATB was studied by Kolb and Rizzo5 in the late seventies and recently (2010) by Sun et al.6 using X-ray diffraction. Over a 303−513 K temperature range, Sun observed significant anisotropy in the thermal response of the crystal with expansion proceeding primarily along the c axis. Plots of the unit cell dimensions as a function of temperature are shown in Figure 10. Coefficients of thermal expansion (CTE) were obtained by fitting the temperature dependent lattice vectors to a second-degree polynomial. The resulting CTE values at 298 K are 1.43 × 10−4, 1.52 × 10−4, and 3.93 × 10−4 K−1 for a, b, and c, respectively, indicating larger expansion along the c direction (factor of 2.5), in line with the experiment. The structural response of TATB under hydrostatic load was simulated from 0 to 10 GPa at 298 K. The resulting pressure dependent unit cell dimensions, angles, and volume are compared to experiments of Olinger and Cady45 and Stevens et al.46 in Figures 11−13. The simulation results overestimate the a vector length (Figure 11) at pressures above ∼2 GPa; however, b and c are in agreement with the experimental data of Olinger and Stevens, respectively, over a larger pressure range. The experimental pressure dependence of the cell angle γ (Figure 12) is well reproduced by the SAPT(DFT) potential, however the α and β angles show deviations of ∼3−6% as the pressure increases (depending on which set of experimental data is used as a reference). Similar deviations are reported by Bedrov29 et al. and as they discuss, there seems to be a systematic discrepancy between the two sets of experimental data which leads to a scatter in the reference values used to assess the performance of the potential. The PV data in Figure 13 was fitted to a third order Birch−Murnaghan equation of state and the theoretical bulk modulus and pressure derivative, along with both sets of experimental values, are presented in Table 5.

expected due to the faster decay of the dispersion interaction (R−6) which dominates the dimer B interaction energy. 3.3. MD Simulation. Using the DL_POLY44 version 2 MD program, the isothermal-isostress (NsT) MD simulation technique (rigid-molecule approximation) was applied to a 6 × 6 × 6 supercell (10 368 atoms) of TATB with the initial supercell configuration taken from ref 1. The interaction potential cutoff was 15 Å, and the long-range Coulombic interactions were handled using Ewald summations. The total simulation time was 2 ps (0.5 fs time step) with velocity scaling at every fifth step (during the first 100 ps) to equilibrate the system. Thermodynamic and crystallographic properties were accumulated over the final 200 000 time steps, and from this data, thermally averaged molecular orientational information (i.e., center-of-mass fractional coordinates and Euler angles) for each molecule was generated. The Euler angles describe the orientations of molecules within the unit cell and were determined using a rotation matrix that transforms the principal axes of inertia of each molecule to the space-fixed coordinate system of the molecule (i.e., the x, y, and z Cartesian axes). Molecular parameters for each of the symmetry-equivalent molecules in the unit cell were averaged over time and over all unit cells in the supercell. Table 3 provides a comparison of time-averaged lattice parameters and the crystal density determined from MD simulations of TATB at 298 K and 1 atm. The cell edge lengths and angles have max errors of 0.48 and 2.7%, respectively, and the crystal density is within 0.62% of the experimental value. This level of agreement is comparable to that obtained empirically by Bedrov et al.29 In order to evaluate the arrangement of the molecules within the unit cell, a comparison of time- and spatially averaged center-of-mass fractional coordinates with corresponding experimental values was performed, and the results are shown in Table 4. There is minimal disagreement between the simulated and experimental orientational parameters and this is further exemplified in Figure 9 where the simulated unit cell is superimposed on the experimental structure with only a slight variation between the two structures. The computed sublimation enthalpy (38.20 kcal/mol) predicted using the SAPT(DFT) forcefield is in good agreement with the experimental values of 40.2 and 43.1 reported by Rosen27 and Garza,28 respectively. The computed error with respect to the Rosen’s value is 4.9%. This is a larger J



Article

Figure 11. Unit cell dimensions as a function of pressure from MD and experiment.

4. CONCLUSIONS This study is additional evidence that the SAPT(DFT) method represents a good quantum mechanical method for generation of intermolecular pair potentials of molecular crystals, and as demonstrated, accurate results that are in good agreement with experiment can be obtained with minimal empiricism. The results presented herein serve as quantum mechanical validation of many of the findings reported previously that were based on energetics obtained using classical forcefields such as Dreiding14 or Lennard-Jones models. The gas phase structure suggested by Roszak et al.13 using MP2 was verified at

the SAPT(DFT) level. The potential parameters developed in this work are specific to the experimental monomer geometry given in Table 1 and may not be transferable to other configurations. Therefore, application of this intermolecular forcefield in simulations where significant intramonomer reorientation occurs is not advised. TATB, much like FOX-7, is a challenging system to model due to the variety of intermolecular forces present within the crystal such as hydrogen bonding and dispersion. In general, wave function based methodologies, which may be too expensive, or dispersion corrected density functionals, which may be of K



Article

Figure 12. Unit cell angles as a function of pressure from MD and experiment.

Figure 13. PV isotherm from MD and experiment.

L



Article

(17) Rai, N.; Bhatt, D.; Siepmann, J. I.; Fried, L. E. Monte Carlo Simulations of 1,3,5-Triamino-2,4,6-Trinitrobenzene (TATB): Pressure and Temperature Effects for the Solid Phase and Vapor-Liquid Phase Equilibria. J. Chem. Phys. 2008, 129, 194510. (18) Bartlett, R. J. In Modern Electronic Structure Theory − Part II; Yarkony, D.R., Ed.; World Scientific: Singapore, 1995. (19) Williams, H. L.; Chabalowski, C. F. Using Kohn-Sham Orbitals in Symmetry-Adapted Perturbation Theory to Investigate Intermolecular Interations. J. Phys. Chem. A 2001, 105, 646−659. (20) Misquitta, A.; Szalewicz, K. Intermolecular Forces from Asymptotically Corrected Density Functional Description of Monomers. Chem. Phys. Lett. 2002, 357, 301−306. (21) Fermi, E.; Amaldi, G. Le Orbite Oos Degli Element. Mem. Acad. Italia 1934, 6, 117−149. (22) Tozer, D.; Handy, N. Improving Virtual Kohn-Sham Orbitals and Eigenvalues: Application to Excitation Energies and Static Polarizabilities. J. Chem. Phys. 1998, 109, 10180−10189. (23) Misquitta, A.; Jeziorski, B.; Szalewicz, K. Dispersion Energy from Density-Functional Theory Description of Monomers. Phys. Rev. Lett. 2003, 91, 033201. (24) Podeszwa, R.; Bukowski, R.; Szalewicz, K. Density-Fitting Method in Symmetry-Adapted Perturbation Theory Based on KohnSham Description of Monomers. J. Chem. Theory Comput. 2006, 2, 400−412. (25) Podeszwa, R.; Bukowski, R.; Rice, B. M.; Szalewicz, K. Potential Energy Surface for Cyclotimethylene Trinitramine Dimer from Symmetry-Adapted Perturbation Theory. Phys. Chem. Chem. Phys. 2007, 9, 5561−5569. (26) Taylor, D. E.; Rob, F.; Rice, B. M.; Podeszwa, R.; Szalewicz, K. A Molecular Dynamics Study of 1,1-diamino-2,2-dinitroethylene (FOX7) Crystal Using a Symmetry Adapted Perturbation Theory-Based Intermolecular Force Fiel. Phys. Chem. Chem. Phys. 2011, 13, 16629− 16636. (27) Rosen, J. M.; Dickinson, C. Vapor Pressures and Heats of Sublimation of Some High-Melting Organic Explosives. J. Chem. Eng. Data 1969, 14, 120−124. (28) Garza, R. G. Thermogravimetric Study of TATB and Two TATB-Based Explosives. Lawrence Livermore Natl. Lab Rep. 1979, 82723, 1−13. (29) Bedrov, D.; Borodin, O.; Smith, G. D.; Sewell, T. D.; Dattelbaum, D. M.; Stevens, L. L. A Moleculary Dynamics Simulation Study of Crystalline 1,3,5-Triamino-2,4,6-Trinitrobenzene as a Function of Pressure and Temperature. J. Chem. Phys. 2009, 131, 224703. (30) Filippini, G.; Gavezzotti, A. The Crystal Structure of 1,3,5Triamino-2,4,6-Trinitrobenzene. Centrosymmetric or Non-Centrosymmetri. Chem. Phys. Lett. 1994, 231, 86−92. (31) van Duin, A.; Dasgupta, S.; Lorant, F.; Goddard, W. A. ReaxFF: A Reactive Force Field for Hydrocarbons. J. Phys. Chem. A 2001, 105, 9396−9409. (32) Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A. ReaxFFlg: Correction of the ReaxFF Reactive Force Field for London Dispersion, with Applications to the Equations of State for Energetic Materials. J. Phys. Chem. A 2011, 115, 11016−11022. (33) Rao, S. ARDEC, Picatinny Arsenal, NJ. Private communication. (34) SAPT2008: “An Ab Initio Program for Many-Body SymmetryAdapted Perturbation Theory Calculations of Intermolecular Interaction Energies” by Bukowski, R.; Cencek, W.; Jankowski, P.; Jeziorski, B.; Jeziorska, M.; Kucharski, S. A.; Lotrich, V. F.; Misquitta, A. J.; Moszynski, R.; Patkowski, K. et al. (35) Williams, H.; Mas, E.; Szalewicz, K.; Jeziorski, B. On the Effectiveness of Monomer, Dimer, and Bond-Centered Basis Functions in Calculations of Intermolecular Interaction Energies. J. Chem. Phys. 1995, 103, 7374−7391. (36) Weigend, F.; Köhn, A.; Hättig, C. Efficient Use of the Correlation Consistent Basis Sets in Resolution of the Identity MP2 Calculations. J. Chem. Phys. 2002, 116, 3175−3183. (37) Misquitta, A.; Podeszwa, R.; Jeziorski, B.; Szalewicz, K. Intermolecular Potentials Based on Symmetry-Adapted Perturbation

Table 5. Comparison of Predicted and Experimental Bulk Modulus (GPa) and Pressure Derivative EOS variable

SAPT(DFT)

expt. (ref 46)

expt. (ref 45)

K (Gpa) K′

12.4 11.0

13.6 ± 0.6 12.4 ± 1.1

16.7 5.7

variable accuracy, are required to study such systems. However, SAPT(DFT) has proven to be a good balance between accuracy and computational cost, as shown in this work and in previous studies.

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Cady, H. H.; Larson, A. C. Crystal Structure of 1,3,5-Triamino2,4,6-Trinitrobenzene. Acta Crystallogr. 1965, 18, 485−496. (2) Kennedy, J. E.; Lee, K.; Son, S. F.; Martin, E. S.; Asay, B. W.; Skidmore, C. B. Second-harmonic Generation and the Shock Sensitivity of TATB. AIP Conf. Proc. 2000, 505, 711−714. (3) Baldridge, K. K.; Siegel, J. S. Balancing Steric and Electronic Factors in Push-Pull Benzenes − An ab initio Study on the MolecularStructure of 1,3,5-Triamino-2,4,6-Trinitrobenzen and Related Compounds. J. Am. Chem. Soc. 1993, 115, 10782−10785. (4) Davidson, A. J.; Dias, R.; Dattelbaum, D.; Yoo, C. Stubborn” Triaminotrinitrobenzene: Unusually High Chemical Stability of a Molecular Solid to 150 GPa. J. Chem. Phys. 2011, 135, 174507. (5) Kolb, J. R.; Rizzo, H. F. Growth of 1,3,5-Triamino-2,4,6Trinitrobenzene (TATB) I. Anisotropic Thermal Expansion. Propell. Explos. 1979, 4, 10−16. (6) Sun, J.; Kang, B.; Xue, C.; Liu, Y.; Xia, Y.; Liu, X.; Zhang, W. Crystal State of 1,3,5-Triamino-2,4,6-Trinitrobenzene (TATB) Undergoing Thermal Cycling Proces. J. Ener. Mater. 2010, 28, 189−201. (7) Manaa, M. R.; Fried, L. E. Nearly Equivalent Inter- and Intramolecular Hydrogen Bonding in 1,3,5-Triamino-2,4,6-Trinitrobenzene at High Pressure. J. Phys. Chem. C 2011, 116, 2116−2122. (8) Ojeda, O.; Cagin, T. Hydrogen Bonding and Molecular Rearrangement in 1,3,5-Triamino-2,4,6-Trinitrobenzene Under Comperssion. J. Phys. Chem. B 2011, 115, 12085−12093. (9) Valenzano, L.; Slough, W. J.; Perger, W. Accurate Prediction of Second-Order Elastic Constants from First Principles PETN and TATB. AIP Conf. Proc. 2011, 1426, 1191−1194. (10) Jeziorski, B.; Moszynski, R.; Szalewicz, K. Perturbation-Theory Approach to Intermolecular Potential-Energy Surfaces of Van-derWaals Complexes. Chem. Rev. 1994, 94, 1887−1930. (11) Taylor, D. E. Prediction of the Impact Sensitivity of Energetic Materials Using Symmetry Adapted Perturbation Theory. Army Res. Lab. Tech. Rep. 2011, 5550, 1−17. (12) Song, H.; Xiao, H.; Hai-Shan, D. Intermolecular Forces and Gas Geometries of TATB Dimers. Acta Chim. Sin. 2007, 65, 1101−1109. (13) Roszak, S.; Gee, R. H.; Balasubramanian, K.; Fried, L. E. Molecular Interactions of TATB Clusters. Chem. Phys. Lett. 2003, 374, 286−296. (14) Zhang, C.; Kang, B.; Cao, X.; Xiang, B. Why is the Crystal Shape of TATB So Similar To Its Molecular Shape? Understanding By Only Its Root Molecule. J. Mol. Model. 2012, 18, 2247−2256. (15) Mayo, S. L.; Olafson, B. D.; Goddard, W. A. Dreiding − A Generic Force-Field for Molecular Simulations. J. Phys. Chem. 1990, 94, 8897−8909. (16) Gee, R. H.; Roszak, S.; Balasubramanian, K.; Fried, L. E. Ab Initio Based Force Field and Molecular Dynamics Simulations of Crystalline TATB. J. Chem. Phys. 2004, 120, 7059−7066. M



Article

Theory with Dispersion Energies from Time-Dependent DensityFuntional Calculations. J. Chem. Phys. 2005, 123, 214103. (38) Perdew, J. P.; Burke, K.; Ernzerhof, M. Generalized Gradient Approximation Made Simple. Phys. Rev. Lett. 1996, 77, 3865−3868. (39) Dalton: A Molecular Electronic Structure Program, Release 2.0, see http://www.kjemi.uio.no/software/dalton/dalton.html. (40) Breneman, C. M.; Wiberg, K. B. Determining Atom-Centered Monopoles From Molecular Electrostatic Potentials − The Need for High Sampling Density in Formamide Conformational Analysis. J. Comput. Chem. 1990, 11, 361−373. (41) Frisch, M. J. et al. Gaussian 09, revision A.1; Gaussian, Inc.: Wallingford, CT, 2009. (42) See the PIKAIA homepage at http://whitedwarf.org/parallel/. (43) Liu, D. C.; Nocedal, J. On the Limited Memory Method for Large Scale Optimization. Math. Program. B 1989, 45, 503−528. (44) Smith, W.; Forester, J. DL_POLY_2.0: A General-Purpose Parallel Molecular Dynamics Simulation Package. J. Mol. Graph. 1996, 14, 136−141. (45) Olinger, B.; Cady, H. H. The Hydrostatic Compression of Explosives and Detonation Products to 10 GPa (100 Kbars) and Their Calculated Shock Compression: Results for PETN, TATB, CO2 and H2O. Proc. 6th Int. Deton. Symp. 1976, 700−709. (46) Stevens, L. L.; Valisavljevic, N.; Hooks, D. E.; Dattelbaum, D. M. Hydrostatic Compression Curve for Triamino-Trinitrobenzene Determined to 13.0 GPa with Powder X-Ray Diffraction. Propellants, Explos., Pyrotech.. 2008, 33, 286−295.

N


(TATB) Using Symmetry Adapted Perturbation ... - ACS Publications

Recommend Documents