Theoretical Charge density Proof of N–N weak bonds of RDX Energetic Molecule

Received : 23-03-2019 Accepted : 12-05-2019 ABSTRACT: The bond topological analysis of Cyclotrimethylene-trinitramine (RDX) energetic molecule has been carried out for the wave function obtained from the ab initio and DFT methods of quantum chemical calculations. The geometrical parameters of all bonds are compared with that of experimental reports. The inclusion of diffuse function in HF basis set levels makes the significant shift of bond critical point towards carbon atoms of C–N bonds. The heteroatomic bond density character is well understood from unequal C-cp and cp-N distances in all C–N bonds. For all the level of calculations, the maximum bond density was found for all N=O bonds, attributes the maximum potential energy V(r). The N–N bond properties are strongly depends upon the equilibrium bond length which clears from charge concentration in shorter N1–N4 bond and charge depletion found in longer N2–N5 and N3– N6 bonding regions. The bond topological analysis of all bonds in RDX molecule resulted that the N–N bond is the weakest among all the other bonds. The weakness of N2–N5 and N3–N6 bonds than N1–N4 bond of RDX has also been analyzed from energy density calculation from various level of theories as an alternate for Laplacian of electron density. From the analysis of CHELPG charges at the MP2 level, the N–N bonds of RDX appears to have a significant ionic nature which attributes strong hyperconjugation effect. The hyperconjugation effect of RDX, due to polarization of N–N bonds, is the additional proof of weak N–N bonds in RDX explosive. The isosurface electrostatic potential shows the electro positive and negative region in the molecule. A large negative potential found at the vicinity of oxygen atoms.


Introduction
To serve fuels or explosives, there has been an extensive search for new high energy density materials (HEDM) for the last couple of decades [1]. Good HEDMs have high density, a fast velocity of detonation (D), which are energetically unstable with respect to their reaction products. The characteristics of such energetic materials depend on various electronic and chemical properties at the molecular level. Such system includes solids, liquids, and gases. Here our studies on RDX are limited to gas phase systems, and primarily consider the electronic properties. For the design of high performance explosives, the analysis of geometry and electronic charge distribution of molecules via theoretical modeling is prerequisite. Several reports [2] includes the characterization of molecules has been carried out using semi-empirical continnum models and further a DFT model also used for these analysis. Therefore, it is important to make a note about the application of electronic and structural parameters in designing the high energetic molecules. The energetic character based on charge density was investigated by Coffey [3] and Kunz and Beck [4]. The Coffey model is fair less specific in providing the precise physical effects of such regions on the molecules constituting the solid. In addition to the incompleteness of the Coffey model, this theory does not provide a convincing model for the presence of the charges found in the fractoemmision from energetics. Some studies performed by Kunz and Beck provide strong evidence that the presence of charges inside an energetic solid can provide significant diminution of the strength of molecular bonds, and in some instances even cause their dissociation. Scheme 1. hexahydro-1,3,5-trinitro-1,3,5-triazine Compounds that contain polynitro groups are always highly energetic. They can provide a large amount of energy and heat after they are burnt and are employed extensively as the main ingredient in explosives. Nitramines have long been used for technological as well as military purposes [5][6][7][8][9]. It is able to withstand much larger mechanical and thermal shocks without igniting. These characteristics make the material particularly well-suited to a variety of defence and civilian applications. [10,11] To this end, it is very important to investigate one of the common very sensitive explosives, RDX (hexahydro-1,3,5-trinitro-1,3,5-triazine) (scheme 1), which is the energetic polynitro compound. In this article, we present a theoretical investigation of RDX. In the attempt to have a better characterization of the geometric parameters, bond topology of RDX, in the present study, we report the results of ab initio calculations [12] at Hartree-Fock(HF) and secondorder Moller-Plesset (MP2) levels using different types of basis sets. For the purpose of comparison and as an alternative to the computationally demanding MP2 methods we have also used density functional theory (DFT), in the Kohn-Sham formulation. [13,14] The major aim of the this work is to provide some insight to understanding the strength of various types of bonds and the energy density distribution of the RDX molecule using AIM theory [15]. Hence, the quantum chemical approaches at various level of sophistication coupled with AIM theory allows to characterize the weak and strong bonds, and helpful to predict the bond break in the molecule while it undergoes decomposition. RDX is a nitramine type highly important explosive [6][7][8]. In practice, usually new energetic materials are designed by modifying known substances by addition and/or modification of energetic group(s) in the molecules.
The theory of atoms in molecules (AIM) [15] calculations were also performed, which permits the investigation of chemical systems on a common basis, as the theory uses only information contained in the electron density ρ(r). The critical points (CP) in the electron density and the Laplacian are the points where ρ(r) vanish. The topology of the Lapalcian field allows one to recover the chemical model of localized bonded and non-bonded pairs and to characterize local concentrations  2 ρ(r) < 0, and depletions,  2 ρ(r) > 0, of the electronic charge density distribution. Bader's theory of Atoms 16 in molecules is the theory of chemical structure and reactivity based on the topological properties of the electronic charge density ρ, the formation of the chemical bond is the result of a competition between the perpendicular contractions of ρ towards the bond path, which lead to a concentration of the charge density along this line and parallel expansion of ρ away from the interatomic surface, which leads to separate concentration in each atomic basin. This behavior results in the formation of a critical point in the charge density at which the Hessian of ρ has two negative eigen values (λ1 and λ2) and one positive eigen value (λ3). This means that ρ exhibits two negative curvatures (λ1 and λ2) perpendicular to the interatomic line and one positive curvature (λ3) along the interaction line . In this theory, the  atomic interactions are classified between two limiting   behaviors: the open and closed-shell interactions. The openshell interactions are characteristics of covalent and polar bonds. In this limiting situation, the charge distribution at the BCP is dominated by the perpendicular negative curvature of the electron density. These open-shell interactions are characterized [17] by large values of ρ,  2 ρ(r) < 0, and |λ1|/λ3 > 1 at the BCP. In contrast, for the closed shell interactions, characteristics of ionic bonds, hydrogen bonds, and van der Waals molecules, the value of ρ is small,  2 ρ(r) > 0 and |λ1|/λ3 << 1. These behaviors can be better understood, if we recall that the local form of the virial theorem [18] can be written as where G(r) > 0 is the electronic kinetic energy density and V(r) < 0 is the electronic potential energy density, defined as the virial of the forces exerted on the electrons. Thus the sign of the  2 ρ serves to summarize the essential physical characteristics of the interactions, which create the BCP. For the closed-shell interactions the kinetic energy density is dominant contribution at the BCP, while for the open-shell interactions the potential energy density makes the prevailing contribution.

Computational Details
The geometry optimization of RDX molecule leading to energy minima were achieved by using quantum chemical calculations including ab initio and Density functional methods (DFT) and these calculations were performed using the GAUSSIAN03 program. [19] Ab initio calculations includes Hartree-Fock (HF) [20] with the basis sets 6-311G** and 6-311++G**. To overcome the electron correlation effects second order Moller-Plesset (MP2) [21] with 6-311G** were employed. Further, DFT calculations include B3LYP, [22] BLYP [23] and BP86 [24] with 6-311G** were also performed. The wave functions obtained from the various optimization procedures were used to calculate the topological properties such as electron density, Laplacian of electron density, and ellipticity at the bond critical points by using the Bader's theory of atoms in Molecules (AIM) implemented in AIMPAC software. [25] The atomic charges derived from NPA, MPA and electrostatic potential derived population (CHELPG) different schemes were calculated for each atom. The deformation densities for each bond of the molecule were plotted by using the software wfn2plots. [26] 3. Results and Discussion But the nitro group at 1-position is essentially coplanar [1.7, 3.7 and 0.8˚] with less deviation from the carbon atom plane.

Structural Aspects
The same trend appears in the molecule optimized at all the levels of theories including the electron correlation effect incorporated MP2 and DFT level of theory. The unique configuration of the 1-nitro groups appears to come from the effect of repulsive non-bonded interaction between adjacent nitro groups. As one may notice in Table 1, a controversial geometric feature among the calculation levels is the N1-N4 bond length. Among the N-N bond length, N1-N4 distance is significantly shorter than N2-N5 and N3-N6 distances attributes the repulsive interaction between the nitro groups attached at N2 and N3 atoms. The N1-N4 bond distance  [27] investigation also reflects the same trend. The N-N bond distances predicted from HF/6-311++G** are found larger, when the electron correlation effect was introduced in the calculation, this distance is further lengthened in DFT level of optimization. The N2-N5 and N3-N6 bond distances predicted for MP2/6-311G** and BP86/6-311G** methods are much longer than experimental [27] crystallographic bond distances [1.392 and 1.398 Å] and are highly overestimated. The differences between the levels are due to the electron correlation and basis set effects. Besides the calculations on N-N bonds, further geometrical calculations also have been made on C-N bonds of RDX. The interesting feature observed in C-N distances for all level of calculations gave the adjacent bond distances in the sixmembered ring are equal, and their distances calculated for both ab initio and DFT theory are C2-N2=C2-N3: [ Due to the crystal field effect, this trend is not found in the experimental [27] reported structures as it differ by 0.01 , 0.014 and 0.003 Å respectively. The C-N-C angles of six membered ring are very close except the angle C3-N1-C1 calculated from MP2, which is too narrowed. The average C-N-C bond angles predicted by the methods HF/6-311++G**, MP2/6-311G** and BP86/6-311G** are 114.8, 112.1 and 115.3˚ respectively. Among these angles, the angle from HF well agree with experimental value. The average N-C-N bond angles for the above levels of calculations are 109.6, 110.6 and 110.6˚ respectively, which are found slightly larger than the experimental value 109.3˚. Among these calculations, HF/6-311++G** calculated C-N-C and N-C-N bond angles are almost agrees with the reported experimental average bond angles are 144.8 and 109.3˚ respectively. [27] From, C-N bonds, the maximum bond twist noticed in C1-N1 and C3-N1 bonds. The torsion angles predicted by HF/6-311++G**, MP2/6-311G** and BP86/6-311G** methods for the bonds C1-N1: 167.5, 169.9, 166.3˚ and C3-N1: -167.5, -169.9, -166.3˚ respectively. Also among the N-N bonds, the maximum bond twist was noticed in N2-N5 and N3-N6 bonds. For N2-N5 bond, it is 168.4, 166.5 and 166.3˚ and for N3-N6 bond 158.1, 151.0 and 159.0˚ respectively. The MP2 and DFT values are expected to be same as in HF but here they differ by ~10˚, this unequal twist may be due to the folding of six membered ring as it is differentiated among the different level of calculations. This inequality of conformation is important. The C-N and N-O bonds are trans oriented with respect to N-N. Table I. Important geometric parameters a of RDX calculated by various levels of theories.   The critical points in the C-N bonds are off from the middle points, but they shifted towards the C-atoms, which indicates that C-N bond densities are highly polarized to C-atom. This effect is more pronounced on the introduction of polarization function in HF calculation. The C-N bond densities obtained from MP2 and DFT methods are slightly smaller except C3-N3 bond (consistently all methods found higher density) on comparing with HF, the corresponding 2 bcp(r) values, which ranges from -13.8 to -18.1 eÅ -5 . Among the DFT calculations, the maximum bond density bcp for all C-N bonds was found in B3LYP/6-311G** level, which was randomly decreased in BLYP and BP86 methods. The unequal C-cp and cp-N distances in all C-N bonds prove the location of the heteroatomic bond density which never lie at the middle of internuclear axis as it is mostly found at the centre in the homo atomic bonds. The charge accumulation in the N=O bonds of three NO2 groups in the molecule are found to be almost equal, despite, the different basis set levels. However, the charge accumulation in these bonds are found to large in HF model calculation. And the MP2 calculation gave the moderate -values as 3.377 -3.384 eÅ -3 and the average is 3.380 eÅ -3 . This value is almost close to the values calculated from DFT level [B3LYP/6-311G**, BP86/6-311G**]. The  2 bcp(r) values for these bonds from MP2 level are -23.5 to -23.9 eÅ -5 . These values are slightly lower than B3LYP but higher than BLYP and BP86 level DFT calculation. The CP's in bonds are found shifted towards N-atoms of N-O bonds, indicates the polarization . Fig 3 (a-c), shows the deformation density of N-NO2 groups in the molecule, and its corresponding Laplacian of density [Fig 3(d-f)] at the bond critical points. Depending upon the method of calculation, the bond properties of N-N bonds varies significantly. The N-N bond densities in the molecule are expected to equal, but, surprisingly the N1-N4 bond density is larger than N2-N5 and N3-N6 bond densities as they are almost equal. At HF/6-311++G** level, the electron density bcp(r) of N1-N4, N2-N5 and N3-N6 bonds are 2.554, 2.472 and 2.472 eÅ -3 respectively. The Laplacian of electron density  2 bcp(r) at the bond critical points of N-N bonds are -22.5, -20.9 and -20.9 eÅ -5 respectively. The bcp(r) and  2 bcp(r) of N-N bonds are found decreased significantly in the electron correlation method MP2 level and the values for the bonds N1-N4, N2-N5 and N3-N6 are 2.188 eÅ -3 , -12.2 eÅ -5 ; 2.08 eÅ -3 , -10.5 eÅ -5 and 2.08 eÅ -3 , -10.5 eÅ -5 respectively. These values are found to be small in DF level of calculation. The decrease in density at CP can be clearly confirmed from the less negative value of Laplacian obtained from DFT calculations indicates, the large charge depletion in these bonds and the values are -8.31(BLYP), -9.5 eÅ -5 (BP86) for N1-N4 and -5.8 (BLYP), -7.0 eÅ -5 (BP86) for N2-N5 and N3-N6. The less negative value of Laplacian, indicates the N-N bond charges are largely depleted on comparing all other bonds in the molecule. This confirms that the bonds are very weak. This prediction supports the structural investigation on RDX molecule. [27] As expected the CP's of the N-N bonds are not at the middle of these bonds, and they are shifted to either sides, the minimum and maximum shift in this bond in HF, MP2 and DFT levels are ranges for N1-N4: 0.072 to 0.024, N2-N5: 0.067 to 0.024 and N3-N6: 0.067 to 0.021 Å. The CP positions of N-N bonds were shifted to middle of the internuclear axis as the density decreases and the bond path length also increased. Based on these results, we conclude that N-N bonds are the weakest bonds among all the bonds in RDX molecule. In the N-N bonds, specifically, the N1-N4 having shorter bond length and exhibit high electron density and high charge concentration. The low electron density and depleted charges were found for N2-N5 and N3-N6 bond having longer bond length. This reveals that, the bond properties strongly depend upon the equilibrium bond length .  Fig 4(a-c) depicts the exact variation of bcp(r) and  2 bcp(r) for various bonds of molecule.

CHARGE DENSITY AND LAPLACIAN OF 
The bond ellipticities is the measure of anisotropy of electron density distribution at CP. It can be calculated from the ratio of the negative values of  = (1/2)-1. The Table 2 reveals the whole spectrum of the shape of electron density distribution at the critical points of bonds in the molecule. The ellipticity values of C3-N1 and C1-N1 bonds are same (0.092 ) and these values are much greater than the values found in the similar type of bonds C1-N2 and C2-N2; C3-N3 and C2-N3 (ranges 0.052-0.056), these low values of ellipticity in C-N bonds indicate that the bond densities in the bond are slightly distorted. Further, on compared to above values with the N-N bonds (average 0.227) and N=O bonds (average 0.117) are much smaller, attributes different bonding nature and shows the anisotropy in the bond densities. The order of ellipticities of the bonds are C-N < N-N < N=O. The trend in DFT calculation is C-N < N=O < N-N and in HF also found the same order. On the whole, the ellipticity obtained from correlation functions incorporated in MP2 and DFT calculations are consistently gave smaller values and specifically, the higher level DFT calculations gave still smaller values. This may be attributed to the less screening of bonding electrons in the molecule, hence the bonding densities preserve the isotropicity even though they are depleted.

Atomic Charges
The NPA MPA, and CHELPG charges have been calculated (table 4) for MP2/6-311G** level and further the group charges of all NO2 groups in the molecule are also calculated and the values are shown in Fig 6(a)-(c). From the NPA and MPA charge layouts, it was found that, the -NO2 fragments bearing slightly negative charges (-0.08, -0.06, -0.06 e and -0.09, -0.07, -0.07 e) and the one for N1, N2 and N3 atoms having highly negative charges, lead no hyperconjugation of the molecule. Going to CHELPG charge layouts, the charge distribution becomes totally different. The CHELPG charges for -NO2 fragments attached at N1, N2 and N3 atoms are -0.04, +0.08 and +0.08 e respectively. The corresponding charge for N1, N2 and N3 atoms are -0.24, -0.47 and -0.47 e respectively. This leads to strong hyperconjugation effect in the molecule and make N-N bonds be highly polarized, especially N2-N5 and N3-N6 bonds, as N2 δ--NO2 δ+ and N3 δ--NO2 δ+ respectively, which confirms the weakness of N-N bonds in the molecule.

Electrostatic Potential
The molecular electrostatic potential (MEP) calculations have been done to predict the polarization, electron correlation and charge transfer effects within the molecule. Fig 7 shows the theoretical MEP obtained from MP2 level of calculation as negative (red) and positive (blue) regions of the property at the +0.5 eÅ-1 and -0.05 eÅ-1 isosurface value. The polarization nature of N-NO2 fragments in RDX and the hypercojugation effect of N-N bonds in the molecule was quite clear (Fig 7), the negative regions are concentrated around the oxygen atoms, while the rest of the molecule has positive ESP. In this iso-surface we noticed one surface of the molecule highly electronegative region and the other surface mounted with positive. This dominant electronegative region may be important for the RDX's extra-ordinary less impact sensitivity compare with other explosive materials.

Conclusion
Depending upon the various ab initio and DFT levels of calculations, the N-N bond length varies significantly, and obtain higher values, when the electron correlation effect is included. N1-N4 bond length is shorter than N2-N5 and N3-N6 bonds for all level of calculations. The shorter N1-N4 bond, have high electron density and significant charge concentration at the bond critical point. On the other hand, the bcp(r) in N2-N5 and N3-N6 bonds are found lesser and the charges are more depleted. This ensures that, predicting the geometric features accurately is of considerable importance in determining the bond properties. This was obvious from energy density calculation, as the N-N bonds have the minimum total energy density H(r) obtained from DFT level of theory. Finally, weakness of N-N bonds, again pointed out from hyperconjugation effect, as N2δ--NO2δ+ and N3δ--NO2δ+ respectively, recommended from CHELPG charge and MEP analysis. On the whole, we conclude that the N-N bonds, especially, N2-N5 and N3-N6 are the weakest bonds of RDX energetic molecule.