Dielectric constant predictions for energetic materials using quantum calculations

Pierre-Olivier Roitaille , Hakima Aou-Rachid , Jos′ee Brisson

a RDDC-Valcartier 2459, Route de La Bravoure, Quebec City, Quebec, G3J 1X5, Canada

b D′epartement de Chimie and CERMA (Centre de Recherche sur Les Mat′eriaux Avanc′es), Universit′e Laval,1045 Avenue de La M′edecine, Qu′ebec, G1V 0A6,Canada

Keywords:DFT calculations Energetic materials Dielectric constant Permittivity constant

ABSTRACT The dielectric constant (DC) is one of the key properties for detection of threat materials such as Improvised Explosive Devices(IEDs).In the present paper,the density functional theory(DFT)as well as ab-initio approaches are used to explore effective methods to predict dielectric constants of a series of 12 energetic materials (EMs) for which experimental data needed to experimentally determine the dielectric constant (refractive indices) are available. These include military grades energetic materials,nitro and peroxide compounds, and the widely used nitroglycerin. Ab-initio and DFT calculations are conducted. In order to calculate dielectric constant values of materials, potential DFT functional combined with basis sets are considered for testing. Accuracy of the calculations are compared to experimental data listed in the scientific literature, and time required for calculations are both evaluated and discussed. The best functional/basis set combinations among those tested are CAM-B3LYP and AUG-ccpVDZm, which provide great results, with accuracy deviations below 5% when calculated results are compared to experimental data.

1. Introduction

The detection of energetic materials is a worldwide major security issue,due to personal mines found on old battlegrounds,but also to the increasing number of artisanal bombs (or improvised energetic material devices, abbreviated IED) planted by terrorists.Several techniques of detections are available[1-6].However,most are not reliable,leading to a large number of false positives[7,8].For a long time, metal detectors formed an acceptable solution. Unfortunately, explosive devices nowadays contain less and less, and often no metal components. Metal detectors therefore fail more and more to detect and avoid hazardous explosive devices, which affects the safety of military and civilian operators, and results in equipment losses.

One of the most promising new techniques to circumvent the shortcomings of metal detectors is the ground-penetrating radar(GPR). This type of radar detects the voltage of returning echoes,which is a function of electrical conductivity and dielectric permittivity variations in soil[9-13].Wavelengths used by the GPR are in the microwave spectral range (from 1 mm to 10 cm), corresponding to the frequency spectrum region of 100 MHz to 10 GHz[14,15].The wavelength is therefore of the same order of magnitude as the IED size,which is a few centimeters long.The power used by a GPR is of the order of a few milliwatts,which is quite lower when compared to that of conventional radars, allowing to sweep the ground without any risk of activating the energetic materials[16,17].

Several improvements are still needed to ensure fully reliable,efficient and accurate detection results,related to temperature and humidity changes in soil under which IEDs are buried [18,19].Nevertheless,this mature technique,produced on a large industrial scale [20,21], has shown encouraging results for the detection of IEDs and other dielectric materials in soil,and even more so when used with a metal detector [12,22,23].

In order to interpret signals in terms of potential IED devices,permittivity reference values for energetic materials are required.Permittivity (ε) is traditionally defined as a constant related to the potential energy (V) of a given system, related to two charges, q1and q2, separated by a distance r, as reported in equation (1) [15].

The dielectric constant εr(or relative permittivity) is the ratio between the permittivity of a material and that of the vacuum,which serves as a reference[24].

The dielectric constantvalue is high for very polar molecules,for example water, or easily polarizable materials, and low for less polar molecules. Although values of relative permittivity are frequency dependent,the zero-frequency relative permittivity can be used due to its small variation over the frequency range covered by GPR. This greatly simplifies the calculations [25,26].

Most energetic materials have relatively similar structures and,therefore, have dielectric constants fairly close to each other,varying between 2.5 and 3.5 [27]. In addition to the difficulties stated above, energetic materials have quite varied chemical structures,and include nitramines,nitrated esters,cyclic peroxides and aromatic nitro,among others.This can lead to issues in the use of a force-field based method, such as molecular mechanics and dynamics. Further, molecular dynamics requires building a large cell, and can therefore be difficult to implement efficiently on a wide variety of compounds [28]. Using a precise method for energetic material detections will contribute to reduce or eliminate false alarms. Accurate experimental data of DC of materials and a detailed interpretation of measured signals will allow to distinguish dielectric constant values of energetic materials from those of the surrounding environments. This aspect can generate another issue, which is the need to compare the measured DC values of materials assumed being energetic materials to those of known energetic materials. This means that an exhaustive dielectric constant database for energetic materials is required.However,due to the number of new energetic materials that are developed in the world, often commercially unavailable, and because of the hazard associated with experimental measurement data of these compounds, building such a database is not straightforward. This specific issue is the one that will be explored in the present paper.

Several methods have been proposed in the literature to predict the dielectric constant of materials. In most cases, a molecular dynamics approach is used [29-31]. A few publications mention the use of the density functional theory (DFT) technique [32-34].However, these references do not systematically propose the testing of various DFT methods and basis sets on energetic materials having diversified structures, and hence do not provide an assessment of the precision afforded by this technique. Other approaches have also been used, such as the theory of effective environments,a technique that takes into account each component of a system and then extrapolates the properties to the overall system[35-37]. Although interesting, these approaches are either imprecise, especially the effective method regarding the surrounding environments,or specific to a particular type of materials.

To overcome these issues,various basis sets in conjunction with different DFT functionals are used in this work. A set of selected energetic materials have been studied,this because of their relative interest in detection domain. Calculated dielectric constant values are compared to values reported in the literature,which will allow to verify the accuracy of the theoretical methods used. The time required to obtain these values will also be compared.This,in order to retain the best methods that are less time consuming while affording the best accurate dielectric constant values of materials.The purpose is to apply these approaches to other energetic materials for which no dielectric constant values are available.

2. Methodology

2.1. Choice of basis sets

Sixteen different levels of calculations were considered,in order to determine which one can provide precise results of predicted dielectric constants and related optical properties in the least possible calculation time.Basis sets which were tested are listed in Table 1,and chosen with the aim of having a fairly varied sampling.Several different functionalities have been used with various basis sets and at various levels of the complexity. The calculations were carried out on 12 different energetic molecules.

Firstly, the ab initio methods HF/STO-3G and HF/6-31G were selected as a “benchmark” for calculation time, as these methods are less time consuming. Their use allows to verify how a substantial gain can be made when going towards more complex methods that are more time consuming.The DFT functional B3LYP was selected, combined with two basis sets having different complexities, 6-311G and 6-311+G (3df, 2p). Functionals using the prefix LC, which provide a better description of longer-distance interactions, were included to test the sensitivity of the dielectric constant determination to this specific approximation, and whether this would increase computational time duration disproportionately. The MP2 ab initio method [38,39]. This method was first developed by Christian M?ller and Milton Plesset in the 1930s[40]. Here, MP2 was tested to demonstrate how the use of the Rayleigh-Schr¨odinger perturbation theory would yield better results. Several tests have been made with basis sets using the Dunning model rather than the Pople model, cc-pVNZ rather than the generic X-YZg approach. The functional developed by Perdew,Burke and Ernzerhof named as (PBE), was also included in comparison with B3LYP. Being also a hybrid functional,it is interesting to compare its performance with B3LYP[41].In the LC-w,terms are added to better consider long-distance interactions. Finally, the CAM-B3LYP[42]functional was included because of its capacity to take into account interaction distances between the atoms.Indeed,depending on the distance,interactions are not treated in the same way by using this functional. For short-distance interactions,greater importance will be given to the B3LYP method(81%against only 19% for HF). The situation is reversed for long-range interactions,65%HF for only 35%B3LYP.This adaptability makes this functional particularly interesting [42].

He was to travel this path now, forno goblet of life is all bitterness; no good man would pour out such adraught for his fellow-man, and how should He do it, Who is lovepersonified?

2.2. Calculation details

Calculations were made using the Gaussian 09 package[43]. In order to run several calculations,and obtain shorter computational time durations, all calculations for the same basis set were launched simultaneously on a supercomputer.For calculations,the geometry startpoint was the conformation of the molecule in the crystal phase, as taken from the literature. Molecular geometries were then optimized to converge to a minimum energy state, and optical properties were calculated. It should be noted that calculations were performed on isolated molecules, without using periodic boundary conditions. Although the molecular conformation remains that of the crystal phase, the exact geometry slightly changed during energy optimization, as expected, due to the absence of interchain interactions in the isolated molecules. For comparison purposes,time durations are reported in time per node(each node has eight cores for the supercomputer used). The ’POLAR’keyword was used to determine the electric field polarization.In order to speed up the convergence a self-consistent field wasused with maximum 600 cycles, Once the polarization of the electric field is obtained, it is possible to calculate the polarizabilities from the matrix provided in the Gaussian software, which in turn is linked to the dielectric constant by the Clausius-Mossotti relationship,

Table 1 Level of theory and basis sets.

where M represents the molar mass, ρ is the density, NAis the Avogadro number and α′represents the polarizability of the material. Therefore, the polarizability estimated by DFT calculations can predict the dielectric constant if the density of the material is known.

Dielectric constant values rarely appear in the literature for energetic materials. However, Maxwell has shown that the refractive index n is, to a very good approximation, equal to the square root of the dielectric constant for a given material and at a specific orientation of the electric field to the material [44]:

It is hence easy to calculate the dielectric constant value from the refractive index of the material.

Refractive indices taken from the literature to calculate the dielectric constant, which will be called ‘experimental dielectric constant values’ throughout this article, resulting experimental dielectric constant values and experimental density values used to calculate the predicted dielectric constants from the DFT-calculated polarizability values are reported in Table 2.

2.3. Selection of energetic material

In order to provide an interesting and application-related comparison between the various calculation methods, it is important to carefully choose the energetic material to be studied.Indeed,a distribution too small or not diversified enough in terms of chemical structures could lead to a false trend. On the other hand, as this is a preliminary work to screen as many datasets as possible,using a larger number of molecules would have taken too much time.We therefore decided to limit the number of molecules,which choosing, at random, molecules which did have varying chemical functionality. The selected molecules also had some interest in connection with their use as an explosive, whether as a military or artisanal one. Finally, experimental data (n and d) was also be available in the open literature. Chosen molecules for this investigation are reported in Fig.1.Abbreviations for the molecules,along with their full names are given in Table 2.

RDX, HMX, TATB, FOX 7, PETN, as well as TNT, Tetryl and Nitroguanidine are all military-grade energetic material.They have been included because, in addition to being commonly used, they are among the best documented energetic materials.HMTD,MEKP and TATP,three peroxides commonly used in improvised explosive devices, are also interesting because of their structures. Unlike other energetic materials in the selection, these are not nitro compounds. Nitroglycerin was added due to its widespread use in many applications. It is also the only liquid in the selection, along with MEKP.

3. Results and discussions

Dielectric constants calculated from DFT results and from experimental density,thereafter designated as calculated dielectric constants, are reported, for each basis set used, in Table 3. These data were subsequently compared to experimental values. Fig. 2 reports the average difference between calculated and experimental values of dielectric constants for all energetic materials studied, along with the time required for calculations. In all cases,when more than one crystal form exists for the explosive, in this table, the crystal form for which experimental data is reported in indicated by the appropriate greek letter in front of the explosive abbreviation.

In Fig.2,several methods show disappointing results.Firstly,as was expected,predictions using the simplest ab initio Hartree-Fock methods are very fast, but show the largest deviations from experimental values (over 20% deviation for HF/6-31G) and, up to 35%deviation for HF/STO-3G,and can only be considered as semiquantitative. The use of the cc-pVDZ basis set in HF/cc-pVDZ provides better results, with approximately 16% average deviations,without significantly increasing the calculation time.

As expected, methods using the hybrid functional B3LYP perform quite well with an average deviation of only 5% when a robust basis set is used such as 6-311+g(3df,2p).The ab initio MP2 method,although used with a fairly simple basis set,shows a very acceptable calculation time, but relatively low accuracy (12.5% deviation). It is also interesting to note that the methods that take better into account long-range interactions (LC and CAM prefix)show significantly better results than others, with deviations of 4-10%). Another important aspect is that the AUG-cc-pVTZ basis set,while being very complete and accurate,is not interesting in a practical point of view,as the calculation time is too long(from 200 to 500 h,which is up to five times longer than most other methods,which average a 100-h calculation time), and accuracy is not significantly better, sometimes worst, than that obtained with the double ζ version (AUG-cc-pVDZ). It is quite interesting to see that the most complex basis sets are not necessarily the best options to consider.

Table 2 Experimental data (refractive index n and density ρ), taken from the literature,and resulting experimental dielectric coefficients of selected energetic materials.

Fig.1. List of modelled energetic material [45].

Table 3 Experimental and calculated values of dielectric constant at different levels of calculations.The symbol?-?represents no DC values for calculations that have not converged.

Four methods show especially good mix of results for the calculation time and accuracy. These are CAM-B3LYP/AUG-ccpVDZ, B3LYP/6-311 + g (3df, 2p), LC-wPBE/AUG-cc-pVDZ and B3LYP/cc-pVTZ. These basis sets present the best results for all molecules,and also the smallest standard deviations.Three of these four combinations are supplemented with a term to better account for long-range interactions, three have the cc-pVNZ basis set, and three are based on the hybrid functional B3LYP.Unsurprisingly,the best option in terms of precision(lowest deviation,with an average value slightly below 5%),CAM-B3LYP/AUG-cc-pVDZ,has all of these three characteristics while having a reasonable calculation time of 50 h, on average.

Fig. 2. Comparison of accuracy and calculation time of DFT using different methods.

An important aspect of interest is the accuracy of the prediction as a function of each of the studied molecules. So far, only the average on all molecules has been considered, but it is worth verifying how accurate results are for each type of molecule studied. Data for all studied methods are reported in Table 3. Upon visual inspection, four methods clearly afforded better results. For these four methods, comparison between calculated and experimental dielectric constant values for each energetic molecule are reported, in Fig. 3 and in Table 3. In Fig. 3 are presented, as single points, data from every method studied in the present work. For comparison purposes, a straight line with slope one was added,consistent to positions where experimental data would be equal to DFT-calculated values. Therefore, cases in which data points are close to this line represent the best matches between experimental and DFT calculated values.In all four cases,the fit is good for most compounds studied.Deviations from the experimental value seem larger for molecules having intermediate values, more specifically RDX,TNT,Tetryl and FOX-7.No clear link between the structures of these different energetic materials can be found.It is therefore not easy to draw a trend for these compounds. It must, however, be kept in mind that experimental data on energetic materials often show large experimental errors, which could partly explain observed differences. This is all the more likely that no basis set seems to be able to adequately predict the dielectric constant values for these compounds.

Fig. 3. Comparison of DFT-calculated and experimental dielectric constant values.

Linear regression was used to assess more precisely which methods afforded the best fit. Regression results for the four best methods are reported in Table 4, where x stands for the experimental dielectric constants and y for the calculated values. The x(experimental dielectric constants) and y (calculated dielectric constants using a given theoretical approach) values are therefore reported in Table 3.The CAM-B3LYP clearly shows the best linearity,with a R2value of 0.94.The slope is also closest to one,which is the value expected for a perfect fit.Nevertheless,other B3LYP methods as well as LC-wPBE/AUG-cc-pVDZ methods show a very good agreement. For purpose of clarification for the readers, Table 5 shows an example of the use that can be made of the linear regression equations reported in Table 4, where y is the theoretically determined value, and x is the dielectric constant calculated using this y value and the linear regression equation. Values are reported for the four best DFT methods considered and for each explosive studied.

In light of the obtained results,the best calculation option is the functional combination of CAM-B3LYP and AUG-cc-pVDZ. In addition to being the most accurate, it remains one of the fastest.Further,as shown in Fig.3,this method provides good results for all molecules considered in this study.Only the LCBLIP/AUG(cc-pVDZ)method provides satisfactory results, particularly for molecules having high dielectric constant values. However, as mentioned earlier,this method is the most time-consuming,and the difference in accuracy of data is not large enough to warrant its use, except perhaps to fine-tune the database for molecules with high dielectric constants.

Fig. 4 further shows a comparison of experimental and calculated values,for each molecule studied,using CAM-B3LYP/AUG-ccpVDZ. Contrary to what one might expect, this method shows no clear tendency to under or overestimate the dielectric constant values.Further,variations are in the same order of magnitude as the error expected for experimental values.

4. Conclusion

The purpose of this explorative study was to develop an approach to determine the dielectric constant values of energetic material of interest, to do so accurately and with a reasonable calculation time.This method is the starting point for developing a future database to identify energetic materials,on the basis of their dielectric constant parameters,as one means among others for the detection of threat materials.Future work should be carried out in order to better understand the influence of medium (soil) temperature on the density, and therefore on the dielectric constant values, of the energetic material of interest.

Table 4 Linear regression results for the best four methods of DFT calculations, where x stands for the experimental dielectric constants and y for the calculated dielectric constant values.

Table 5 Values of x parameters (Experimental values) for each explosive and method.

Fig. 4. Comparison between CAM-B3LYP/AUG-cc-pVDZ values and experimental values.

In this work, explorative theoretical calculations were chosen.Their use circumvents two problems: The hazards associated with experimental measurements of energetic material and the limited accessibility of some energetic material, which becomes a nonlimiting factor. Among the available techniques, the DFT approach has emerged as the best option because of its accuracy, versatility,relative rapidity,and ease of execution.Sixteen basic and functional combinations were compared with a set of thirteen energetic materials. Comparison was made in terms of calculation accuracy, as compared to experimental data from the scientific literature, and calculation time. Following this comparison, it became clear that the best combinations among those tested was CAM-B3LYP/AUGcc-pVDZ. In addition to being the most accurate with an average error under the five percent mark, it is also one of the fastest.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.