Charge Distributions of Nitro Groups Within Organic Explosive Crystals: Effects on Sensitivity and Modeling

The charge distribution of NO2 groups within the crystalline polymorphs of energetic materials strongly affects their explosive properties. We use the recently introduced basis-space iterated stockholder atom partitioning of high-quality charge distributions to examine the approximations that can be made in modeling polymorphs and their physical properties, using 1,3,5-trinitroperhydro-1,3,5-triazine, trinitrotoluene, 1-3-5-trinitrobenzene, and hexanitrobenzene as exemplars. The NO2 charge distribution is strongly affected by the neighboring atoms, the rest of the molecules, and also significantly by the NO2 torsion angle within the possible variations found in observed crystal structures. Thus, the proposed correlations between the molecular electrostatic properties, such as trigger-bond potential or maxima in the electrostatic potential, and impact sensitivity will be affected by the changes in conformation that occur on crystallization. We establish the relationship between the NO2 torsion angle and the likelihood of occurrence in observed crystal structures, the conformational energy, and the charge and dipole magnitude on each atom, and how this varies with the neighboring groups. We examine the effect of analytically rotating the atomic multipole moments to model changes in torsion angle and establish that this is a viable approach for crystal structures but is not accurate enough to model the relative lattice energies. This establishes the basis of transferability of the NO2 charge distribution for realistic nonempirical model intermolecular potentials for simulating energetic materials.


I. ROTATING MULTIPOLES FOR BOTH ORIENT AND DMACRYS CALCULATIONS
The ORIENT 1 program can transform the iterated stockholder atoms 2, 3 (ISA) distributed multipoles of a molecule (originally calculated in the global axis frame of the molecule) into the atomic local axes of the molecule as defined in Figure S1. Figure S1. An example of the global (molecule fixed) axis (red) and local axis (blue) definitions used in this study. In the global axis frame, z is perpendicular to the ring plane, thus the π orbitals give a significant Q 20 moment (c.f. the dz 2 orbital). When the quadrupole tensor is rotated into the local axis, the z axis is now along the bond, hence the Q 20 moment is in the ring plane and has a very different value, with the π orbitals being represented by Q22s. This is done by analytically rotating the global multipoles into the local axis frame. The molecular geometry and the user defined local axes must first be defined. The multipole moments calculated in the global axis are then analytically rotated into this user defined local axis.

II. DIFFERENCES IN OPTIMIZED AND EXPERIMENTAL CONFORMATIONS OF RDX CONTRASTED WITH TNT
Overlays of experimental conformations, the PBE0/aug-cc-pVTZ 4-7 optimized conformations, and the conformations with only the NO2 angles rotated (optexptNO2) are shown in Figure S2, to illustrate the changes that accompany nitro-group rotation. In the AAE conformation of RDX, the equatorial NO2 moves more into the plane of the ring, as the aliphatic ring is more open for the isolated molecule, possibly due to repulsion from lone pairs on the ring sp 3 nitrogens, that are not being counteracted by crystal packing forces. A similar change in the aliphatic ring conformation and nitrogens in the axial NO2 groups is seen for the AAA conformation, therefore, it is not possible to consistently analytically rotate the nitro-group multipole moments. While these changes may not look severe, we can see by a comparison of differences in cell parameters, RMSD15 values and lattice energies between using the experimental and optimized molecular conformations in Table S4 that these differences result in very different minima. In contrast, the aromatic ring in the optimized and experimental conformations of the nitro-aromatic TNT is unchanged, with only small methyl group rotation ( Figure S2), meaning analytical rotation is feasible. Table S3 shows that although there is a significant difference in the nitro groups, which does have a major effect on the crystal packing, changing the nitro-torsion angles to the experimental values can result in a good reproduction of the structure.

III. SURVEY OF CAMBRIDGE STRUCTURAL DATABASE
Distribution curves illustrating the occurrence of each ---torsion in a specific chemical environment (for example, an NO2 adjacent to two hydrogens or a methyl group in an aromatic ring) were generated using an in-house python program developed by Luca Iuzzolino for determining the possible torsion angles that could be adopted in crystal structures based on the torsion angle statistics of all 870,000+ molecules in the Cambridge Structural Database (CSD). 11 The specific NO2 fragments using in this study are as follows The occurrence of each torsion angle (rounded to the nearest degree) was plotted in the histograms. The optimized isolated structures and the most stable experimental conformations within the polymorphs were used. In RDX this is the AAE conformation in CTMTNA03, and the AAA conformation in CTMTNA06. The range of (X is either C or N, the atom connected to the nitro-group) and is the min and max values from the 3NO2 groups (6 for HNB). The impact sensitivity used is given with other impact sensitivity determinations in brackets; for RDX the polymorph used is not specified, thus it is probable that the impact sensitivity of the form is very different from literature values, if the experiments used the most stable form. 1 of the longest bond. 2 The trigger bond with the most positive .

V. EFFECTS OF ROTATING NITRO-GROUPS IN TNB, HNB AND RDX (AAE)
The relative changes of the charge and dipole moment magnitudes with NO2 torsion angle are compared with Cambridge Structure Database (CSD) distributions and the intramolecular energy relative to the optimized isolated molecular structure for the molecules where this data is not given in the manuscript.  A comparison of the expt (experimental structure), xminexpt (DMACRYS 17 minimized crystal structure, with all molecules held rigid in their experimentally observed conformations), xminopt (DMACRYS minimized crystal structure, with all molecules held rigid in their gas-phase optimized conformation), xminoptexptNO2 (DMACRYS minimized optimized structure with experimental NO2 torsions using computed after rotation) and the xminanarot (DMACRYS minimized optimized structure with experimental NO2 torsions using analytically rotated multipole moments) crystal structures and intermolecular lattice energies of TNT, TNB and HNB. The intermolecular lattice energy ( ) estimates above are calculated using a distributed multipole electrostatic force-field derived from the molecular (PBE0/aug-cc-pVTZ) wave-function and empirical repulsion-dispersion model; the empirical FIT potential 18  (b) calculated for the PBE0/aug-cc-pVTZ optimized structure and then analytically rotated (anarot) into the experimental observed NO2 torsion angles (grey). 14 = 0.214Å Figure S4 shows that even for TNT, which has the greatest differences in optimized and observed NO2 torsion angles, analytical rotation of the atomic multipoles into the experimental NO2 torsion angles and recalculation of the wave-function for each NO2 torsion angle results in very small structural differences. The only significant differences are seen in the intermolecular energies (Manuscript Table 1 & Table S3).

Table S4
and Figure S5 highlight the differences between the experimental and optimized conformations of RDX and how it would be inappropriate to use molecular charge distribution models derived from the optimized structure to predict experimental crystal properties such as lattice energy or cell geometries. A comparison of the expt (experimental structure), xminexpt (DMACRYS 17 minimized crystal structure, with all molecules held rigid in their experimentally observed conformations) and xminopt (DMACRYS minimized crystal structure, with all molecules held rigid in their gas-phase optimized conformation) crystal structures and intermolecular lattice energies of RDX. The intermolecular lattice energy ( ) estimates above are calculated using a distributed multipole electrostatic force-field derived from the molecular (PBE0/aug-cc-pVTZ) wave-function and empirical repulsion-dispersion model; the empirical FIT potential 18 & ISA 2 (PBE0/aug-cc-pVTZ). This highlights how slight differences in the experimentally observed and optimized conformations of both polymorphs of RDX can result in large changes in intermolecular energy.