Study of the Elastic Properties of the Energetic Molecular Crystals Using Density Functionals with van der Waals Corrections.

We studied the elastic properties of crystalline energetic materials within the framework of density functional theory with van der Waals interactions (DFT-D3(BJ) and rev-vdW-DF2). The full sets of elastic constants were computed. The computed parameters are in good agreement with the experimental data. Among the crystals studied in this work, FOX7 had the lowest compressibility value of 0.0034 GPa-1 and had the highest anisotropy. Crystalline pentaerythritol tetranitrate had almost isotropic mechanical properties.

Introduction

Energetic materials are of great interest for practical applications; there has been a steady interest in them for decades.[1−6] Moreover, the interest will be maintained as long as there are practical applications for these compounds, since there are still a lot of unanswered questions related to the stability and decomposition of energetic materials. This problem traditionally has been of great interest, thus experimental data for energy materials, obtained using precision measurements, are widely available.[7−14] Analysis of the elastic properties of energetic materials is of great interest. There is a lot of experimental data for some crystals of energetic materials, despite the fact that they have low symmetry, which means it is difficult to determine experimentally the full set of elastic constants. Furthermore, with the development of new measurement methods, the data are regularly updated.

The projection of 2,4,6-trinitrotoluene (TNT) unit cell (21/a)[15] on the ab-plane is shown in Figure . Figures S1–S7 show unit cells of the pentaerythritol tetranitrate (PETN), erythritol tetranitrate (ETN), 1,3,5-trimethylene trinitramine (RDX), 1,3,5,7-tetranitro-1,3,5,7-tetrazoctan (HMX), 1,3,5-triamino-2,4,6-trinitrobenzene (TATB), nitromethane (NM), and FOX7. Weak intermolecular forces play an important role in the low-pressure region. Determination of elastic constants makes it possible to calculate the mechanical properties of crystals, which allows us to study the crystal reaction to compression along various directions. Since molecules of energetic materials have both planar and nonplanar structures, this presents additional interest in understanding the role of van der Waals forces in molecular crystals.

Figure 1

Crystal structure of TNT viewed along the a-axis (a). TNT molecule in the crystal (b).

Crystals of energetic materials such as PETN,[7] RDX,[12−14,16] and HMX[8,9,17] were also studied in sufficient detail. The symmetry of their unit cells is higher than monoclinic ones. The absence of a monoclinic angle greatly simplifies the task, since the number of elastic constants decreases. Thus, low symmetry, defects, and weak van der Waals forces make experimental measurements of elastic constants for molecular crystals difficult. Nevertheless, for example, a series of studies are required for the RDX crystal to determine the correct values for a full set of elastic constants.[11−14] Thus, it is of interest to study the elastic properties of energetic molecular crystals within a single method. Therefore, theoretical calculations of the elastic constants of molecular crystals are very useful. Note that the theoretical values of the elastic constants strongly depend on the accuracy of the calculation.[18] The inclusion of intermolecular forces in the DFT theory allows such calculations to be performed.[19−21]

The interest in the structures, where van der Waals interactions play a key role, prompted the development of schemes that take into account dispersion forces within the framework of DFT: (i) fully ab initio approaches (i.e., nonempirical), (ii) the reparameterization of the existing functionals, and (iii) the inclusion of empirical terms. One of the goals of this work is to study the possibilities of DFT-D3(BJ)[22,23] and rev-vdW-DF2[24] for determining a full set of elastic constants. The generic type of energetic materials has been selected for this purpose. Elastic constants are extremely sensitive to calculation accuracy; thus, the obtained results will be able to demonstrate the schemes’ possibilities. DFT-D3(BJ) and rev-vdW-DF2 are often used to obtain relaxed structures.[25−27] Direction dependence of linear compression will demonstrate the methods’ possibilities, as it will allow us to demonstrate how van der Waals forces are taken into account along various directions in a crystal. Mechanical characteristics of energetic molecular crystals have also been determined. There is no data on the elastic constants on some of the compounds under study.

Computational Details

A plane-wave pseudopotential approach within DFT was used to compute the total energy. The Quantum ESPRESSO (QE)[28] with the functional of Perdew, Burke, and Ernzerhof (PBE),[29] was used to carry out the computations. The ultrasoft pseudopotentials of the Rab–Rape–Kaxiras–Joannopoulos type were used for calculations.[30] The crystal structures were optimized with the Broyden–Fletcher–Goldfarb–Shanno (BFGS) method.[31] The energy cutoff equals 65 Ry. The Monkhorst–Pack scheme[32] was used for the Brillouin zone sampling (Table ).

Table 1

k-Point Grid for the Molecular Crystals of Energetic Materials

crystal	formula	k-grid
TATB	C6H6N6O6	2 × 2 × 3
PETN	C5H8N4O12	2 × 2 × 3
ETN	C4H6N4O12	1 × 3 × 1
TNT	C7H5N3O6	2 × 3 × 1
RDX	C3H6N6O6	2 × 2 × 3
HMX	C4H8N8O8	3 × 2 × 3
NM	C5H8N4O12	3 × 3 × 2
FOX7	CH3NO2	3 × 3 × 2

Geometry relaxation was completed when all components of all forces are smaller than 0.1 mRy (a.u.)−1. Structural data were used as the starting point.[3,15,33−38] Dispersion interaction is caused by the correlation of electrons and many-body effects, which are not taken into account in the existing DFT approximations.[39] The possibility of including van der Waals interactions in DFT has been suggested by Kohn and co-workers,[40] and the van der Waals density functional (vdW-DF)[41] is the most well-known first principles scheme. We used the fully nonlocal exchange–correlation rev-vdW-DF2 functional in our research.[24,41−45] Also, we used DFT-D3(BJ).[22,23,46,47] This semiempirical scheme is notable for being simple and robust. DFT-D3(BJ) almost does not increase the computation time. In this scheme, the empirical potential was added to the exchange–correlation potential and the total energy is given bywhere EKS-DFT is the Kohn–Sham energy and Edisp is a dispersion correction.[22]

The dispersion energy iswithHere, the sum is over all atom pairs in the crystal. C6 and C8 are the isotropic dispersion coefficients for atom pair AB and RAB is the internuclear distance. s6 and s8 are the global (functional dependent) scaling factors. The detailed information and values of all parameters can be found in the original works.[22,23] The calculations of the elastic constants were carried out using QE and thermo_pw package.[48−50]

Apart from TATB (triclinic symmetry), the energetic molecular crystals under consideration have orthorhombic, tetragonal, and monoclinic symmetry. The linear compressibility for crystals with a monoclinic symmetry of the unit cell is written as followsl1, l2, and l3 are the directional cosines with respect to the x, y, and z axes and s is the elastic compliances. The components tensor A

Results and Discussion

Tables S1–S8 present the calculated lattice parameters. Experimental values are also listed. Table S9 presents the values of the elastic constants for TATB. Among other energetic materials under study, TATB is the only crystal that has a triclinic symmetry of unit cell, that is why a complete set of elastic constants contains 21 parameters. There is no data on the experimental study of the elastic constants for this crystal. This is due to the low crystal symmetry, which makes experimental study difficult. However, there is theoretical data that has been obtained with various schemes,[20,51] including molecular dynamics calculations.[52] There is a good agreement between the theoretical results obtained in this work and the existing literature data.[52] The obtained values of the elastic constants make it possible to determine the values of linear compressibility along various directions in a crystal (Figure ). Projections of the direction-dependent linear compressibility in the xy-, xz-, and yz-planes of the TATB crystal are shown in Figure . Direction-dependent linear compressibility and bulk modulus of the energetic molecular crystals are shown in Figures S8–S10.

Figure 2

Illustration of the direction-dependent linear compressibility of TATB (GPa–1).

Figure 3

Projections of the direction-dependent linear compressibility of TATB (GPa–1) in different planes.

The maximum compressibility is observed along the direction in which the molecules approach and slide relative to each other. Weak van der Waals forces are the main contributors to the retention of molecules. Thus, even a small force causes molecules to approach each other. Molecular layers are displaced in such a way that an α angle increases, which is confirmed by the experimental and theoretical data.[53,54] When compressed in the ab-plane, the crystal demonstrates a large resistance because the approaching molecules cause repulsive forces to increase. Thus, the shortest distance between the molecules, located in one layer, equals ∼2.4 Å (H–O). It means that the molecular attraction is accompanied by a strong repulsion caused by an overlap of atomic orbitals. The maximum compression equals 0.055 GPa and is almost the same as that of graphite (0.041 GPa).[19] The minimum compression can be observed in the molecular plane and equals 0.012 GPa. This value is certainly greater than the compression along the graphene layer (0.001 GPa).[19]

Tables S9–S16 present the calculated elastic constants for PETN, ETN, RDX, TNT, HMX, FOX7, and nitromethane. Both experimental and theoretical data are reported for some of these crystals. There is a reasonable agreement between the theoretical and experimental data. It should be noted that due to experimental complexity in determining the elastic constants, experimental data is scattered. Also, the calculations were performed at T = 0 K. Besides, some of the terms are not used in DFT-D3. It is important to take into account several physical effects: zero-point vibration, anharmonicity, Pauli repulsion, and the many-body nature of dispersion.[55] Therefore, the absence of these terms also leads to the discrepancy between experimental and theoretical data.

As in the case of TATB, linear compressibility along different directions was studied using the obtained results. Projections of linear compressibility on various planes are shown in Figure .

Figure 4

Projections of the direction-dependent linear compressibility (GPa–1) of energetic materials in different planes.

On the whole, both methods (DFT-D3(BJ) and rev-vdW-DF2) predict accurate values that are in good agreement with the literature data. Thus, modern density functionals with van der Waals corrections predict accurate values of elastic constants. Therefore, these schemes can be used for determining the mechanical characteristics of energetic materials. It seems quite relevant as it is difficult to study molecular crystals with low symmetry of the unit cell. As can be seen from the study of RDX,[12−14,16] the use of modern precision experimental methods has improved the early results. The result demonstrates the complexity of determining elastic constants for structures where van der Waals forces play the key role. It should be noted that RDX has an orthorhombic unit cell, although molecular crystals more often have a lower symmetry. As has already been said, there are certain difficulties. A larger crystal should be used to obtain more accurate results; however, it increases the number of defects. Besides, the low symmetry of the molecular crystals contributes to the complexity of the task. Therefore, even in the case of the RDX crystal, there are certain difficulties. There is a noticeable difference between the data obtained via Brillouin scattering and the data received using resonant ultrasonic. At the same time, the data obtained via resonant ultrasonic also varies. So, the results obtained by different groups for the RDX crystal have significant differences. From this point of view, it is possible to compare the data for different crystals using ab initio calculations with the same parameters.

Energetic materials have different forms and symmetry of molecules, which is also manifested in the symmetry of unit cells. Some basic mechanical characteristics were calculated using the Voight–Reuss–Hill approximation for ease of comparison (Table ). As the calculations are performed within the framework of a unified theory, it enables a better understanding of the differences between different crystals.

Table 2

Isotropic Aggregate Elastic Properties Based on the Voigt–Reuss–Hill Averagesa

	B, GPa	E, GPa	G, GPa	μ
TATB	16.96	23.07	9.06	0.27
PETN	9.73	15.40	6.23	0.24
ETN	7.73	13.02	5.34	0.22
TNT	10.65	14.16	5.54	0.28
RDX	12.47	16.95	6.66	0.27
HMX	13.39	17.78	6.97	0.28
NM	10.00	12.26	4.73	0.30
FOX7	17.23	19.20	7.30	0.31

The bulk modulus (B), Young’s modulus (E), shear modulus (G) and Poisson ratio (μ) for energetic materials are computed using DFT-D3(BJ).

Table presents the values of linear compressibility along different crystallographic axes. The minimal and maximal values of linear compressibility are quite suitable for the comparison of the properties of different crystals (Table ). It is clearly seen that crystals of energetic materials have different responses to compression. The PETN crystal is almost isotropic. The PETN molecule has a tetragonal symmetry; therefore, they exhibit the same resistance to external forces along all directions. Thus, the isosurface has a shape very similar to a sphere.

Table 3

Values of Linear Compressibility (GPa–1) along Crystallographic Axes are Computed Using DFT-D3(BJ)a

crystal	symmetry	βa	βb	βc	βmin	βmax
TATB	triclinic	0.0118	0.0118	0.0423	0.0102	0.0546	5.3529
PETN	orthorhombic	0.0339	0.0339	0.0350	0.0339	0.0350	1.0324
ETN	monoclinic	0.0410	0.0496	0.0307	0.0275	0.0549	1.9964
TNT	monoclinic	0.0261	0.0207	0.0437	0.0207	0.0491	2.3720
RDX	orthorhombic	0.0240	0.0191	0.0388	0.0191	0.0388	2.0314
HMX	monoclinic	0.0231	0.0344	0.0143	0.0098	0.0344	3.5102
NM	orthorhombic	0.0235	0.0351	0.0423	0.0235	0.0423	1.800
FOX7	monoclinic	0.0162	0.0513	0.0039	0.0034	0.0513	15.0882

The minimum and maximum values of linear compressibility in the crystals.

The properties of the ETN crystal are similar to those of PETN,[25] but its reaction to compression is different. A preferential direction is observed in the ac-plane. The approximate angle formed with the a-axis is 49°. Since the ETN molecule has an inversion center and its molecular form is different, it results in the formation of a monoclinic angle, which is manifested in the response of the crystal to compression. The TNT molecule has an aromatic ring, while the RDX molecule contains a ring consisting of three carbon atoms and three nitrogen atoms. Unlike the TATB molecule, these molecules have atoms located in various planes.

The RDX and nitromethane crystals have an orthorhombic unit cell. The maximal compression for both is implemented along the c-axis. The minimal compression for a crystalline RDX is observed along the b-axis and for nitromethane along the a-axis. The maximum compression for FOX7 is βb = 0.0513 GPa–1. The molecule layers approach each other when compressed along this axis. Van der Waals forces play the main role in holding the layers together in this direction. It is similar to the way molecular layers approach each other in TATB crystal. It should be noted that the minimum value is very small βmin = 0.0034 GPa–1, it is 15 times smaller than βmax. The direction with the minimum compression is very close to the c-axis. The compression along this direction results in molecule deformation, which causes high resistance. It explains the high C11 and C33 values. Thus, FOX7 crystal has the highest anisotropy among all of the crystals of energetic materials under investigation.

Conclusions

Dispersion forces play a key role in molecular packing in molecular crystals. High compressibility is observed along directions where dispersion forces hold the molecules together. There are voids between molecules along these directions. Thus, even slight pressure causes molecules to approach each other. The macroscale manifestation of it is anisotropy of the mechanical properties of the crystals, especially the properties of the TATB and FOX7 crystals.

We used two different schemes (ab initio and semiempirical) to account for the interaction of van der Waals forces. Full sets of elastic constants were determined within the framework of a single approach, which allowed us to compare them with each other. The mechanical properties of typical energetic molecular crystals were also calculated. Since the experimental data is not available for TATB, ETN TNT, FOX7, and NM, this allows a better understanding of the properties of these crystals.

The volume of the unit cells computed in the framework of rev-vdW-DF2 has lower values than that computed with DFT-D3 (BJ). This leads to the fact that the elastic constants calculated with rev-vdW-DF2 generally have a higher value than those calculated with DFT-D3 (BJ).

Thus, it has been established that using density functionals with van der Waals corrections makes it possible to correctly determine the full set of elastic constants for molecular crystals with different symmetries of unit cells. This allows us to understand the impact of mechanical deformations on the properties of molecular crystals of energetic materials. Since the synthesis of new energetic materials is continuing at present, the use of these schemes will make it possible to determine the necessary mechanical characteristics for energetic materials. This is especially useful when studying the properties of new energetic materials, which are often unstable under the impact of mechanical deformations.