Electronic structure of solid nitromethane: Effects of high pressure and molecular vacancies

The combined effect of pressure and molecular vacancies on the atomic structure and electronic properties of solid nitromethane, a prototypical energetic material, is studied at zero temperature. The self-consistent charge density-functional tight-binding method is applied in order to investigate changes induced in the band gap of this system by uniform and uniaxial strain of up to 70%, corresponding to static pressure in the range of up to 200 GPa. The effects of molecular vacancies with densities ranging from 3% to 25% have also been considered. A surprising ﬁnding is that uniaxial compression of about 25–40 GPa along the b lattice vector causes the C–H bond to be highly stretched and leads to proton dissociation. This event also occurs under isotropic compression but at much higher pressure, being indicative of a detonation chemistry which is preferential to the pressure anisotropy. We also ﬁnd that the band gap, although evidently dependent on the applied strain, crystal anisotropy and vacancy density, is not reduced considerably for electronic excitations to be dominant, in agreement with other recent ﬁrst-principles studies. © 2002 American Institute of Physics. @ DOI: 10.1063/1.1466830 #


I. INTRODUCTION
A series of experiments have investigated the effects of stress and shocks on energetic materials, [1][2][3][4][5][6] demonstrating strong correlations among the applied mechanical stress, the intrinsic solid-state properties and the onset of chemical reactions that may lead to detonation. 7These observations have in turn stimulated theoretical interest into the microscopic nature of the detonation initiation, in which a combination of shock parameters such as duration and pressure cause widespread chemical reactions.Despite concentrated effort, [8][9][10] the mechanisms by which the propagating shocks affect atomistic conditions that favor bond-breaking and subsequent detonation are not as yet fully understood.A longterm objective is to construct reliable models that encompass nonlinear dynamics, dispersion and chemical ͑or phase͒ transformations, 11,12 in order to predict aging and other properties of energetic materials pertaining to safety.
Experiments 1,4,13 suggest that local heterogeneity and defects in solids, such as dislocations, vacancies, microvoids, impurities, and cracks, play an important role in the sensitivity and performance of energetic materials to shock or impact. 146][17] The nature of these hot spots is an unresolved issue.Dlott and Fayer 16 have proposed that the anharmonic coupling between phonons and low-frequency molecular vibrations is strong at such hot spots, and hence causes defect molecules to transiently attain high vibrational temperatures; these temperatures favor chemical reactions that normally do not occur in the bulk.However, no widely agreed upon microscopic model exists to describe how the size, morphology, and density of the defects in the crystalline solid affect the overall rate of the reaction energy release.
Under shock propagation or impact, the energetic solid experiences a sudden compression of up to 30% of its original volume.At least by intuition, one expects that the atomic and electronic properties of a heterogeneous molecular solid could change dramatically and alter its overall mechanical response at high pressures.Several proposals currently exist to place this speculation on firm physical grounds.According to Williams, 18 initiation of detonation in solids can occur through changes induced in the Fermi level because of a double electron and positive hole injection followed by chemical reactions.Coffey 19  process is due to tunneling of dislocations in a solid; when the shock velocities are sufficiently high, the dislocations can have energy adequate to directly pump the internal vibrational modes of the constituent molecules.More recently, Dremin 20 made an attempt to explain why molecules lying in the shock front dissociate; his proposed dissociation mechanism directly involves electronic excitations.In the same spirit, Gilman [21][22][23][24] suggested that the compression due to the pressure wave front causes local metallization when the bending of certain covalent bonds reduces the energy gap between the highest occupied ͑HOMO͒ and the lowest unoccupied molecular orbitals ͑LUMO͒.A candidate bond angle for bending is that between the C-N bond and the plane of the almost omnipresent nitro (NO 2 ) group in explosives.This effect in turn favors chemical decomposition, in compliance with a criterion given in the 1920s by Herzfeld 25,26 for crystals with dense structures.According to this criterion, a material becomes metallic if its molar refractivity, which is proportional to polarizability and hence approximately proportional to the inverse of the optical gap, exceeds the molar volume. 248][29][30][31] These authors make use of the first-principles Hartree-Fock approach with an approximate treatment of electron correlations via second-order many-body perturbation theory ͑MBPT͒. 32Their results have indicated that compression of the RDX crystal in the presence of single 29,30 and dimer 29,31 vacancies reduces the optical gap of this material appreciably, thus decreasing the excitation energy needed for the insulator-metal phase transition.Specifically, these calculations showed that the edge dislocations cause a dramatic reduction of the optical gap due to the splitting of the local electronic states from both the valence and conduction bands.These findings have in turn prompted a mechanism based on electronic excitations induced by the impact wave propagating through the crystal. 33,34According to this hypothesis, the pressure exerted by the impact wave front causes the dramatic reduction of the band gap to nearly zero values and results in the breakage of the N-NO 2 chemical bond in RDX, thus initiating detonation and chemical chain reactions. 34ecently, Reed et al. 35 examined the nature of electronic excitations in crystalline solid nitromethane (CH 3 NO 2 ) under static and dynamic compression by using densityfunctional theory ͑DFT͒ calculations in the context of the generalized-gradient approximation ͑GGA͒.These authors conclude that the band gap is not lowered sufficiently to produce considerable thermal population of the electronic excited states in the crystal, in apparent contradistinction to the results of Kuklja and Kunz, 29,30 who dealt with different systems and larger primitive cells.Reed et al. 35 also referred to earlier calculations by Manaa and Fried, 36 who considered the role of nonradiative transitions in isolated nitromethane molecules and indicated that metallization in the gas phase can be achieved by bending of the nitro group; this group originally lies in the CNOO plane of the singlet ground state.][39][40][41][42][43] In this paper, we study the combined effect of static pressure and molecular vacancies on the optical gap of the nitromethane crystal for volume strains up to 70% by following the behavior of the band gap at the single kϭ0 point ͑⌫͒.This gap is known to be related to the HOMO and the LUMO electron states that are localized near the nitro group.There are at least two reasons to choose nitromethane as a prototypical energetic material: the first is that the nitromethane molecule provides the simplest model of a C-(NO 2 ) bond in energetic materials; the second is that, as mentioned previously, recent first-principles calculations 36 involving the single molecule of nitromethane have reported that singlet-triplet minimum energy crossings require a significant bending of the nitro group off the CNOO plane of the ground state.It is therefore expected that the crystalline form of nitromethane is a good candidate to check Gilman's hypothesis, [21][22][23][24] or other theories on the hot spot formation.In particular, a question addressed in this paper is whether an appreciable bending of the NO 2 group inside the crystal can be the primary cause of band-gap reduction within the experimentally accessible pressure range.
Our electronic-structure calculations are based on the self-consistent charge density-functional tight binding 44 ͑SCC-DFTB͒ scheme.This is an extension of the standard tight binding ͑TB͒ approach 45 in the context of DFT, 46,47 and self-consistently describes total energies, atomic forces, and charge transfer.This method has been tested for systems such as small organic molecules and biomolecules, 44,48 and semiconductors. 44,49An essential element of our analysis is that, for each fixed set of lattice parameters, all interatomic positions are determined by total energy minimization, i.e., the atoms are fully relaxed.
The tight-binding character of the SCC-DFTB makes it feasible to study large supercells, containing up to 32 molecules, with concentrations of molecular vacancies ranging from 3.125% to 25%.The present study therefore differs with respect to computational tools and system sizes from recent first-principles calculations, 35 which considered smaller supercells with more demanding, first-principles computations.The two methods give very similar results for the band gap of systems accessible to both, as discussed in Sec.III.In particular, our calculations give a small reduction of the HOMO-LUMO gap and no sign of significant distortion of bond angles involving the nitro group for uniform volume strain up to 70%.For the present study, the tendency of the SCC-DFTB to overestimate band gaps is not a handicap, since only the band gap change is of interest, as a function of the applied strain.To corroborate this claim we showed that, by imposing the bending of the nitro group in the nitromethane molecule, the SCC-DFTB method indeed predicts a significant reduction of the HOMO-LUMO gap in remarkable agreement with high-level first-principles calculations. 36 noteworthy finding of our static calculations is that the C-H bond is highly stretched at a relatively early stage of uniaxial compression along the direction of the b lattice vector, a prediction which is made for the first time by simulations to our knowledge.As the pressure increases, this bond is stretched further, while its hydrogen atom becomes acidic and tends to dissociate.Accordingly, a serious reduction of the band gap at high pressures is favored by the formation of nearly free electrons and protons.Our simulations under isotropic compression indicate a similar tendency for proton dissociation at a much higher hydrostatic pressure, pointing to the preliminary conclusion that the chemistry onset is strongly affected by pressure anisotropy.
Molecular-dynamics simulations can reveal changes in the atomic and electronic structure that may not be captured by static calculations. 35However, our purpose here is to probe the atomic and electronic properties in an adiabatic sense, as a first step toward more accurate, first-principles simulations for large systems.A prerequisite at this stage is to understand how the electronic structure depends on macroscopic parameters, such as strain and stress, in the simplest possible context.In this connection, we also need to clarify how reliable our method actually is.
The remainder of this paper is organized as follows: In Sec.II the key ideas of the SCC-DFTB method are outlined, along with a brief description of its computational details and of the first-principles adaptive-coordinate real-space electronic-structure [50][51][52] ͑HARES͒ method, used here on a limited scale for comparison purposes.In Sec.III, the SCC-DFTB method is applied in order to calculate the total energies and the HOMO-LUMO energy gap of the single nitromethane molecule, the band gap at ⌫ of the perfect periodic crystal, as well as of crystals with molecular vacancy densities of 3.125%, 6.25%, 12.5%, and 25%.As further discussed in Sec.III, these calculations show a not strictly monotonic behavior of the gap as a function of the vacancy density for fixed compression, with a distinct dependence on the crystal anisotropy.We find that sufficient uniform or uniaxial compression causes the C-H bond to be highly stretched, before any significant bending of the nitro group occurs.Finally, in Sec.IV we conclude the paper with a summary of our findings and comparisons to other calculations and experimental data available in the literature.We alert the reader that the terms ''band gap'' and ''HOMO-LUMO gap'' are used interchangeably in the present work.

II. COMPUTATIONAL METHODS
The key idea underlying the SCC-DFTB method is the formulation of the traditional TB approach 45,53 in terms of a variational principle; this gives the TB approach as a stationary approximation to DFT. 54 -57 Accordingly, the SCC-DFTB method is based on a second-order expansion of the total energy functional from DFT with respect to charge density fluctuations ␦n relative to a reference density n 0 .The first ͑zeroth-order͒ term amounts to the standard TB approach.More precisely, the total energy is where H ˆ0 is a single-particle Hamiltonian operator resulting from the input density n 0 , i are the single-particle electron wave functions corresponding to occupied states of valence electrons, E rep stands for the combined ion-ion repulsion and double-counting Hartree and exchange-correlation terms represented by an effective short-range pair potential, 57 and is the correction due to density fluctuations ␦n.It is worthwhile noting that ⌬E (2) accounts for the charge transfer among atoms and renders the proposed scheme a selfconsistent extension of TB. 44 ⌬E (2) is formally defined as 44 In the above expression, E XC is the exchange-correlation ͑XC͒ energy contribution within DFT. 57 By expanding i ϭ i (r) in a suitable set of nonorthogonal localized atomic orbitals , the total energy becomes where is an approximation to Eq. ͑2͒ that is found by decomposing the charge density fluctuations into atom-centered contributions, treating them within a multipole expansion and retaining monopolar terms only.The first term in Eq. ͑4͒ is the band-structure energy, while ⌬q ␣ ϭq ␣ Ϫq ␣ 0 are the atomcentered charge density fluctuations in the usual Mullikencharge analysis. 58The Hamiltonian matrix elements H 0 are calculated in the ''two-center approximation'' 55 using a minimal set of atomic orbitals.Equation ͑5͒ represents the long-range Coulomb interactions between point charges at different sites under the monopole approximation, and includes the self-interaction contributions of single atoms. 44By minimizing the right-hand side of Eq. ͑4͒, one derives Roothaan-type equations 58 for the coefficients c i which become self-consistent through the modification of the Hamiltonian matrix elements because of the nonzero term ⌬E (2) .E rep is determined by appropriate fitting of data for the difference between the cohesive energy versus distance in the local density approximation ͑LDA͒ and the corresponding TB band-structure energy for a suitable reference structure. 44n order to apply the SCC-DFTB to the nitromethane crystal, the TB Hamiltonian and overlap integrals have been evaluated as functions of distance between the requisite O, N, C, and H atom types.In all calculations of this paper, we use minimal basis sets consisting of s and p atomic orbitals for O, N, C, and only the s orbital for H.
We have also used the adaptive-coordinate real-space electronic-structure ͑HARES͒ code, 50 in order to calculate the HOMO-LUMO energy gap in the case of the single nitromethane molecule for comparison purposes.The band-gap calculation within this approach was supplemented by a lowest-order correction on the basis of a theory originally formulated by Fritsche, 59 who generalized the Kohn-Sham theory by considering changes in the electron density of the ground state due to electron excitations from the HOMO to the LUMO.

A. The nitromethane molecule
The structure of the nitromethane molecule in its singlet ground state is shown in Fig. 1͑a͒.The molecule exhibits C S symmetry, with a dominant single closed-shell configuration character. 60,61The constituent atoms were fully relaxed by SCC-DFTB and the C atom was found to be coplanar with the NO 2 group.Molecular orbitals of importance are the and * of the C-N pair, the ''inplane'' NO 2 bonds, the and * of the NO 2 group, and the and * of the C-H pair ͑for more details see Ref. 36͒.The starting point for the calculations is the equilibrium structure given by Manaa and Fried, 36,62 which was obtained via the multiconfiguration self-consistent field 63 ͑MCSCF͒ method.The bond lengths and angles determined after atomic relaxation by SCC-DFTB are within 2-4% of the values given in their Table 1. 62The bending angle ␥ depicted in Fig. 1͑b͒ is defined as the angle between the plane of the NO 2 group and the original CNOO plane of the fully relaxed structure.
The HOMO-LUMO gaps corresponding to the planar and bent structures, for ␥ϭ53°, are shown in Table I evalu-ated separately by the SCC-DFTB and HARES methods.A lowest-order correction has also been supplied ͑fourth column of Table I͒ to account for the rearrangement of the electron occupied states due to excitations, according to Fritsche's generalization of the Kohn-Sham theory. 59This correction is calculated using the electron density and singleelectron wave functions from HARES within LDA. 52 The value given in the last column of Table I for the undistorted molecule is in very good agreement with the electron-impact spectroscopy experiments by Flicker et al., 64 where a singlet-triplet transition was observed to have an onset at 3.1 eV and an intensity peak at 3.8 eV.Note that the SCC-DFTB overestimates the HOMO-LUMO gap by nearly 1 eV when ␥Ϸ0, as expected. 44he geometry corresponding to the second row of Table I has been chosen to be very close to the geometry of the minimum energy crossing determined recently by Manaa and Fried 36 by use of MCSCF, with all atomic positions kept fixed during the calculations.The values reported here are in qualitative agreement with their results, according to which there is a singlet-triplet minimum energy crossing point for ␥ close to 50°.We emphasize that the precise location of this minimum energy crossing of course depends on the method employed.The essential point is that the SCC-DFTB correctly indicates a tendency toward closure of the HOMO-LUMO gap when ␥ deviates significantly from zero.

B. Perfect crystal 1. Crystal structure
Nitromethane is known to be a liquid at room temperature.Its crystalline form is employed here for the sake of simplicity; it also makes meaningful future comparisons to properties of other energetic materials, such as RDX and PETN, which are known to be solids at room temperature.The crystal structure of nitromethane at nearly zero pressure is orthorhobic, with 4 molecules per unit cell. 65The values used here for the initial lattice constants are aϭ5.198Å, b ϭ6.266 Å, cϭ8.649Å ͑space group P2 1 2 1 2 1 ͒, which have been determined via lattice relaxation by starting with the experimental lattice parameters and using the CASTEP plane-wave code.The initial volume of the primitive cell is approximately V 0 ϭ281.7 Å 3 .These numerical values are within 2% of the experimental ones, determined by x-ray single crystal diffraction and neutron diffraction data. 65,66he bond lengths and angles found after full atomic relaxation by SCC-DFTB are within 2-4% of those determined by MCSCF for the nitromethane molecule in the gas phase. 62niform or uniaxial compression was applied to the unit cell by successively decreasing the lattice parameters and relaxing the atomic positions at each step.Use was made of the conjugate gradient technique for the atomic relaxations, with a threshold of the order of 10 Ϫ4 a.u.for the maximum magnitude of the atomic forces.The resulting change in the total energy per unit cell is shown in Fig. 2 for compression up to 50% of the original volume.Because the initial lattice parameters of the primitive cell have not been relaxed by the SCC-DFTB method, the total-energy curves of Fig. 2 are found to be monotonically decreasing in the volume strain (V 0 ϪV)/V 0 for sufficiently small values of the strain, where V is the volume of the unit cell at each step of compression.The pressure, estimated by invoking the formula P ϭϪ‫ץ‬E tot /‫ץ‬V where E tot is the total energy per unit cell, is relatively small and negative for this range of strain.P becomes zero when the total energy attains its minimum at V/V 0 Ϸ0.76 for uniform compression, and V/V 0 Ϸ0.83, 0.72, and 0.91 for uniaxial compression in x, y, and z, respectively.For higher values of the strain ͑smaller values of V/V 0 ͒, the total energy is of course found to increase with compression and P is positive and increasing, expressing the strength of the intermolecular interactions which progressively develop as the molecules approach each other.The unit cell of the nitromethane crystal is shown in Fig. 3͑a͒.Before compression, the HOMO-LUMO gap is approximately equal to 4.6 eV, as seen in Fig. 3͑b͒.This value is close to that given in Table I for the single nitromethane molecule, because of the weak intermolecular forces in the initial primitive cell.

Uniform compression
To simulate the effect of uniform strain on solid nitromethane, the initial lattice parameters were successively decreased by keeping their ratio fixed.Throughout this paper, the pressure P is estimated by using the low-temperature formula PϭϪ‫ץ‬E tot /‫ץ‬V.The initial pressure is found to be of magnitude about 1.6 GPa, close to the value given by recent first-principles calculations. 35For hydrostatic compression up to 50% of the original volume V 0 , the pressure rises to about 50 GPa, while the HOMO-LUMO gap dropped by 0.6 eV, i.e., by 13% of its original value.This compression results in a simultaneous increase of the HOMO and the LUMO energies ͑Kohn-Sham eigenvalues͒, while the decrease of the HOMO-LUMO gap is almost monotonic in the volume strain, as depicted in Fig. 3͑b͒.The change induced in the band gap by very high hydrostatic compression is shown in Fig. 4 for strain equal up to 70%.
It is of interest to note that the mutual orientations of the C-N axes vary smoothly with the strain, (V 0 ϪV)/V 0 , when this is less than 40%, with the corresponding pressure not exceeding 20 GPa.When V/V 0 is between 58% and 60% and P is estimated to lie between 15 and 25 GPa, the atomic configuration undergoes an abrupt change accompanied by rotations of the methyl (CH 3 ) groups, in agreement with first-principles calculations; 35 a similar transition was observed in that work for V/V 0 between 59% and 77% and P in the range 10-30 GPa.This transition is expected because the barrier for the rotation of the methyl group is known to be very low. 66he following changes of the bond lengths and bond angles were observed as the strain varied from 0 to 50%.FIG. 2. Total energy per unit cell of the nitromethane perfect crystal under uniform ͑hydrostatic͒ and uniaxial compression of the unit cell along x ͑for the lattice parameter a͒, y ͑for b͒, or z ͑for c͒.V is the volume of the compressed unit cell.The geometry of the initial unit cell of volume V 0 is very close to the experimental one.Since no lattice relaxation is applied, each energy curve is monotonically increasing in V/V 0 when this is sufficiently close to 1.The corresponding pressure is negative but small.While practically no bending of the nitro group in the molecules took place, the C-N bond lengths were invariably shortened at most by 6%.In contrast, the eight N-O bond lengths in the unit cell exhibited different variations, with a maximum deviation of Ϯ1% from their equilibrium values.The bond angles involving the nitro group in each molecule remained practically intact.A few of the C-H bonds were stretched up to 2% of their initial values.No significant change in the bond angles involving the nitro group was noticed even as the uniform strain was increased from 50% to 70%.However, two of the C-H bonds were highly stretched when the strain was extended to values above 60%.The same effect on the C-H bond in all molecules of the unit cell was also observed under uniaxial stress, yet the estimated pressure is appreciably lower.This case is discussed in more detail below.

Uniaxial compression
Compared to uniform compression, uniaxial strain along one of the lattice vectors is more likely to lead to detonation.Dick has proposed 67 that detonation initiation in nitromethane is favored by shock-wave propagation in specific directions related to the orientation-dependent steric hindrance to the shear flow.This proposal is based on a model according to which the sterically hindered shear process causes preferential excitation of optical phonons strongly coupled with vibrons.
We applied uniaxial compression along each of the x, y, and z axes, and tracked the HOMO-LUMO gap as a function of the volume change, for strains not exceeding 50% and corresponding estimated pressures as high as 100 GPa ͓see Fig. 3͑b͔͒.Among the three types of compression, the ones in x and y yielded the highest estimated pressure.The gap was decreased at most by 1.4 eV ͑a 30% reduction of the original value͒ under stress in y.This drop is higher than in the case with uniform strain, while the gap behavior here is not strictly monotonic, i.e., the gap does not continuously de-crease as the uniaxial strain increases.The strains in x and y seem to cause a similar overall behavior of the gap.The compression in y, however, causes a steeper reduction of the gap as V/V 0 approaches 50%.In contrast, the compression in z reduces the energy gap more drastically for intermediate values of V/V 0 , roughly between 65% and 80%.This last effect can be visualized as follows.The molecules approach each other in z by bringing the methyl group of one close to the nitro group of the next, while the cell lattice parameter c approaches the intermolecular distance.Consequently, the electron densities of the HOMO and the LUMO, which are localized near the nitro group, are distorted more than in any other type of compression of similar magnitude and the band gap is reduced more rapidly for this type of uniaxial strain.
A closer look at the atomic configurations reveals structural transitions which are triggered preferentially by pressure anisotropy.The most frequent transition in the x compression involves of course molecule translations, rotations of the methyl groups ͑which occur even for V/V 0 close to 1͒ and reorientations of the molecules, primarily via relative translations and rotations around their C-N axes.Two such abrupt transitions are observed when V/V 0 ϭ67-70% and the estimated pressure is 25-30 GPa, and V/V 0 ϭ59-62% with an estimated pressure of 45-55 GPa.In the z compression the methyl groups seem to undergo smoother transitions as each molecule's center-of-mass is translated in z.In contrast, the y compression causes abrupt transitions when V/V 0 ϭ76-81% and 63-66% and the pressure is estimated to lie in the ranges 2-5 and 13-16 GPa, respectively.
Notably, the dramatic stretching of four C-H bonds, one in each nitromethane molecule of the unit cell, occurs under stress in y when V/V 0 is between 59% and 62%, with an estimated pressure of 25-40 GPa; more precisely, these bond lengths are stretched by more than 10-12% of their original values in this case.The increase of strain in y above 50% causes these bonds to be stretched further and leads to the abstraction of their protons.This indication of proton dissociation corresponds to the steep part of the curve in Fig. 3͑b͒, and renders the y and x compressions qualitatively different.
The C-N bond lengths are again invariably shortened, at most by 6%, 9%, and 2-3% of their initial values under x, y, and z compression, respectively.In contrast, the N-O bond lengths exhibit positive or negative deviations from their equilibrium values that depend both on the location of the corresponding molecules and the direction of compression.It is worthwhile noting that only the compression in y causes the N-O bond length to be continuously stretched in all molecules, at most by 3% of the initial value.The bond angles do not exhibit any change worthy of reporting for this range of compression.

C-H bond high stretching under uniaxial compression
As mentioned previously, the SCC-DFTB calculations predict a high stretching of the C-H bond that leads to proton abstraction under y compression when the strain becomes of the order of 40% or higher; for uniform compression, the corresponding strain takes the significantly higher values, in the range 60-62%, at estimated pressure of 150-180 GPa.One of the intermediate atomic configurations leading to the proton dissociation in the former case of uniaxial strain is shown in Fig. 5.Note that the high stretching of the C-H bond is visualized graphically by enlarging the symbols that correspond to each abstracted acidic hydrogen atom ͑the atomic sizes in these figures are proportional to the van der Waals radius of the ionized form͒.
A physical picture of the C-H bond stretching is the following: upon contraction of the initially single C-N covalent bond, electronic charge is transferred to the region of space between the C and N atoms.Consequently, the C-N bond becomes stronger and its character approaches that of a double covalent bond.The adjacent C-H bonds act as sources of the electronic charge; at least one of these H atoms is deprived of some electronic charge and its bond is weakened.On the other hand, the orbitals localized near the NO 2 group are not significantly distorted by the compression.The Mulliken charge changes given in Table II are in agreement with this picture.The total Mulliken charge loss of the H atoms is found to increase monotonically with uniaxial strain in each molecule, being dominated for high strain by the charge pertaining to the markedly stretched C-H bonds ͑shown in Fig. 5͒, as expected.While these values follow closely the Mulliken charge gain of the C atoms for low uniaxial strain, they start to deviate as the strain approaches 50%.
These findings imply that chemical reactions and transitions from covalent to ionic C-H bonds prior to detonation in nitromethane are more likely to start before the closure of the optical gap.Roughly speaking, an actual precursor to detonation appears to be the proton dissociation and not the hypothesized metallization induced by drastic changes in the bond angles. 22Interestingly, energetic materials can therefore be considered as more sensitive to their chemistry under uniaxial rather than uniform compression, in agreement with Dick's 67 elaborate results that invoke a steric hindrance model for nitromethane.
There have been numerous experimental studies of the reaction kinetics that may lead to detonation of liquid nitromethane. 68 -70Most of these studies are concerned with the nature of chemical products that may speed up detonation.Shaw et al., 68 for example, measured the effect of the concentration of protonated compounds within the energetic material on the time needed until explosion.Their experi-FIG. 5. Atomic configurations of the unit cell for uniaxial compression along bϭby ˆ, with the nitromethane molecules numbered from 1 to 4 in the clockwise sense ͑top left molecule is numbered 1͒.In each molecule, the H atoms are designated by the letters a, b and c ͑not to be confused with the lattice parameters͒.The x axis is perpendicular to the plane of the page.͑a͒ Original unit cell of volume V 0 ϭ281.7 Å 3 .͑b͒ Configuration for V/V 0 ϭ0.79; the cell undergoes phase transformations ͑molecule reorientations and methyl-group rotations͒.͑c͒ Configuration for V/V 0 ϭ0.60; the C-H bond lengths have been stretched by more that 12% of their original values.͑d͒ Configuration for V/V 0 ϭ0.51.In ͑c͒ and ͑d͒ the dissociated acidic H atoms ͑protons͒ are shown as larger light gray circles.TABLE II.Mulliken charge changes ͑in units of e͒ of the C and H atoms of the nitromethane molecules 1 and 2 of the unit cell under uniaxial compression along b, in close correspondence to Fig. 5.All values are relative to the Mulliken charges in the primitive cell of volume V 0 ϭ281.7 Å 3 .Each H ic atom (iϭ1,2) is found to become acidic and start dissociating from the corresponding C atom when V/V 0 ϭ0.60.The Mulliken charges for the other two molecules ͑3 and 4͒ of the unit cell show an essentially identical behavior and their values are therefore omitted from this table.ments indicated that a reaction involving hydrogen atoms or protons is involved in the first steps of this fast-reaction process.This suggestion was later corroborated by Blais et al. 69 by more direct measurements, who also referred to a sequence of previous investigations by Engelke and others. 70ur zero-temperature, static calculations point to similar conclusions for the nitromethane molecular crystal.A similar prediction for this material was also reached concurrently 71 via ongoing first-principles, molecular-dynamics simulations at high densities ͑1.5-2.5 gr/cm 3 ͒ and high temperatures (T ϭ2000-4000 K).We plan to supplement our results with first-principles static calculations in the near future. 72

C. Effect of molecular vacancies
In view of Dick's 67 proposal and the persistence of uniaxial strain following shock propagation, the cases of uniform and uniaxial compression are herein presented separately.In order to simulate vacancies of desired densities in the nitromethane crystal, we first superimposed two, four, and eight unit cells in order to create 2ϫ1ϫ1, 2ϫ2ϫ1, and 2ϫ2ϫ2 supercells where the numbers express multiples of the primitive unit cell dimensions.More precisely, a vacancy concentration of 25% was formed in two different ways: first, by removal of one entire molecule from the primitive unit cell and, second, by removal of two molecules from the 2ϫ1ϫ1 supercell.Vacancy concentrations of 12.5%, 6.25%, and 3.125% were also modeled by removal of one molecule from the 2ϫ1ϫ1, 2ϫ2ϫ1, and 2ϫ2ϫ2 supercells, respectively.In the following, unless it is stated otherwise, there is no lattice relaxation and the initial supercell structure merely stems from a multiple of the starting primitive cell along the x, y, and z axis.
The unit cell with a molecular vacancy is shown in Fig. 6͑a͒ and the corresponding HOMO-LUMO gaps for uniform and uniaxial compression are given in Fig. 6͑b͒.The 2ϫ1 ϫ1 supercell models are shown in Figs.7͑a͒ and 8͑a͒ with their HOMO-LUMO gaps depicted in Figs.7͑b͒ and 8͑b͒.The larger supercells are not shown.A pair of vacancies in the 2ϫ2ϫ1 supercell is created by removing the molecules nearest to the endpoints of the cell diagonal.To realize some of the underlying complexity, it suffices to mention that a 2ϫ2ϫ2 supercell with one molecular vacancy contains 31 nitromethane molecules, which amounts to 217 atoms with 744 valence electrons.Each supercell is strained up to 50% of its original volume.FIG. 6. ͑a͒ Unit cell of the nitromethane molecular crystal with 1 molecular vacancy, identified by a large circle ͑25% vacancy concentra-tion͒.͑b͒ The HOMO-LUMO gap ͑in eV͒ for uniform and uniaxial compression of the unit cell with 1 molecular vacancy.FIG. 7. ͑a͒ 2ϫ1ϫ1 supercell of the nitromethane crystal with 1 molecular vacancy, identified by a large circle, and relaxed atomic positions, corresponding to a 12.5% vacancy concentration.͑b͒ HOMO-LUMO gap ͑in eV͒ for uniform and uniaxial compression of this supercell.

Uniform compression
The HOMO-LUMO gap with one molecular vacancy ͑a concentration of 25%͒ in the primitive unit cell is shown in Fig. 6͑b͒.As V/V 0 is changed from 100% to 50%, the pressure is estimated to vary from 1 GPa to 15 GPa.These values are lower than in the case of uniform compression of the perfect crystal, as expected.A comparison with Fig. 3͑b͒ shows that the band gap exhibits a tendency to remain flat for a wider range of the volume strain in the presence of a vacancy, but the overall drop is about the same as in the case with no vacancy, i.e., 0.6 eV.This resistance of the gap to changes induced by stress stems from the increase in the size of the effective free space that each molecule experiences because of the vacancy.The gap reduction is so small that it cannot allow for considerable electron excitations from the HOMO to the LUMO.The molecules reorient themselves so that the structure relieves the applied stress; roughly speaking, the molecules tend to take positions that are closest to those of minimum pressure for the given lattice parameters.The most frequent structural change observed is again the rotation of the methyl group, yet the variation in the atomic positions appears to be smoother in the presence of a vacancy.
The effect of vacancy concentrations of 12.5% and 25% on the HOMO-LUMO gap are shown in Figs.7͑b͒ and 8͑b͒ on the basis of a 2ϫ1ϫ1 supercell.With V/V 0 varying from 100% to 50%, the estimated pressure reaches almost 30 GPa for concentration of 12.5% and 10 GPa for concentration of 25%.Note the differences from the values in Fig. 6͑b͒, which may reflect the anisotropy of the crystal since the relative positions of the vacancy are different in the two cases.Again, no appreciable narrowing of the band gap is found.The small jumps appearing in the simulation data for the band gap correspond to abrupt structural transitions involving relative reorientations of the C-N axes and rotations of the methyl groups, especially of the molecules that lie in close proximity to the vacancies.
In Fig. 9͑a͒ the HOMO-LUMO gap of the 2ϫ1ϫ1 supercell is depicted for vacancy concentrations of 12.5% and 25%, with the inclusion of partial optimization of the initial lattice parameters.Accordingly, V 0 is now the volume that corresponds to the relaxed supercell in each case.Notably, any two of these curves intersect.In other words, for any fixed V/V 0 , the gap is not a strictly monotonic function of the vacancy concentration.The pressures at the intersection points for the perfect crystal are estimated to be 5 GPa and 18 GPa ͑from right to left͒.In Fig. 9͑b͒ the gaps are shown for vacancy concentrations of 12.5% and 6.25% from the 2ϫ2ϫ1 supercell, and 3.125% from the 2ϫ2ϫ2 supercell, where the starting point (V/V 0 ϭ1) is that of the minimum pressure for the given vacancy concentration and hydrostatic compression.The curves are well-separated now.For any fixed volume strain, the gap is roughly a monotonically decreasing function of the density.In all these cases, the gap is reduced at most by 1.3 eV, dropping to 3.25 eV.The smallest FIG. 8. ͑a͒ 2ϫ1ϫ1 supercell of the nitromethane crystal with 2 molecular vacancies, identified by large circles, corresponding to a 25% vacancy concentration.Note that the arrangement of vacancies in the periodic extension of the supercell is different from the one of Fig. 6. ͑b͒ HOMO-LUMO gap ͑in eV͒ for uniform and uniaxial compression of this supercell.FIG. 9. HOMO-LUMO gap ͑in eV͒ for uniform compression and different vacancy concentrations: ͑a͒ Perfect crystal, and vacancy concentrations of 25% and 12.5% from an approximately optimized 2ϫ1ϫ1 supercell with defects.͑b͒ Vacancy concentrations of 3%, 6.25% and 12.5% from the 2ϫ2ϫ1 and 2ϫ2ϫ2 supercells with defects.decrease in the gap occurs when the vacancy concentration is the highest, i.e., equal to 25%.

Uniaxial compression
Application of uniaxial strain does not alter the order of magnitude of the HOMO-LUMO gap reduction, as shown in Figs.6͑b͒, 7͑b͒, 8͑b͒, and 10 where only strain in the x direction is considered.In fact, for this uniaxial compression the gap is reduced at most by 1 eV.The unstrained structure of Fig. 10 corresponds to the minimum pressure for the given compression in x by starting with the nonoptimized lattice parameters of the unit cell.In each set of curves in Figs.10͑a͒ and 10͑b͒ the smallest change in the gap is again observed when the vacancy density is the highest, 25%.In Fig. 10͑a͒ one sees that the gap decreases with the vacancy density, and thus the curves appear to be more distinct, compared to the case with uniform strain.

IV. SUMMARY AND CONCLUSIONS
We studied the combined effect of static pressure and molecular vacancies on the crystal geometry, atomic configuration and electronic structure of the prototypical solid nitromethane.We found that exerting pressure on the crystal does not cause a dramatic band-gap reduction associated with the bending of the nitro group.
As a step toward simulating the effect of shocks, we applied uniform and uniaxial compression to the nitromethane crystal, with and without molecular vacancies.Our calculations demonstrated only a slight reduction of the HOMO-LUMO gap for uniform compression down to 30% of the original volume and hydrostatic pressures up to 150 GPa or higher.This small reduction essentially prohibits any consideration of electron excitations, and hence of nonradiative transitions due to the crossing of adjacent energy surfaces of molecules. 36The band-gap drop is most likely attributed to the relatively weak intermolecular interactions, as the molecules tend to maintain positions that yield the lowest possible pressure in the crystal.Accordingly, no appreciable bending of the nitro group is observed for this range of volume strains.
In a physical picture that can explain our results, the constituent molecules undergo changes in their bond lengths, relative orientations of the C-N axes, and rotations of their methyl groups in order to relieve the applied stress.What configuration prevails of course depends on how the energy landscape evolves under adiabatic compression, and therefore on how high the barriers are for transitions to neighboring local minima of other configurations.
The most common structural transition observed in our simulations is the rotation of the methyl group.We emphasize that it is crucial to be able to describe accurately the atomic forces and to include full atomic relaxation, an essential feature of the SCC-DFTB method, which makes it possible to capture this effect.This aspect of the calculations renders our analysis distinctly different from that by Kuklja et al., 34 who applied a ''rigid molecule approximation'' and thus calculated a large decrease of the band gap in the RDX crystal.Direct comparisons of our results to theirs are of course not compelling because of the different systems involved.However, we expect that allowing for full atomic relaxation necessarily leads to a band-gap behavior that precludes dominant electron excitations in either system.
Because of the crystal anisotropy, on the other hand, the sequence and the frequency of phase transformations depend on the character of compression.For example, the compression in z forces the molecules to develop more rapidly electrostatic interactions, which in turn affect the electronic charge of the nitro group and therefore result in a steeper decrease of the band gap with applied pressure.Since charge transfer is important in the description of atomic forces, we expect that its realistic description within the SCC-DFTB method captures the general features of the molecular response to external stress.
It follows from the physical picture outlined previously that the stress relief is in general facilitated by the presence of vacancies, since the effective size of free space in the unit cell increases.Indeed, our results show a higher overall persistence of the band gap of the crystal with vacancy defects with applied pressure.It should be pointed out, however, that the present model of periodically distributed vacancies is restricted in its applicability; in real materials the vacancies follow an inhomogeneous, and even random, distribution.
We also paid attention to the distortion of bonds under compression, even above 100 GPa, of the unit cell or supercell.The C-H bond in particular was found to be highly stretched, and its proton most likely dissociates, in the per-FIG.10.HOMO-LUMO gap ͑in eV͒ for uniaxial compression and different vacancy concentrations, in close correspondence to Fig. 9. ͑a͒ Perfect crystal, and vacancy concentrations of 25% and 12.5% from an approximately optimized 2ϫ1ϫ1 supercell with defects.͑b͒ Vacancy concentrations of 3%, 6.25%, and 12.5% from 2ϫ2ϫ1 and 2ϫ2ϫ2 supercells with defects.
fect crystal under both uniform and uniaxial strain in the y direction, but in the latter case the requisite pressure is estimated to be significantly lower, lying in the accessible regime of 20-40 GPa.We interpreted this effect as being associated with the induced electron-charge transfer from the hydrogen to the region between the carbon and the nitrogen, since the C-N bond is noticeably contracted and its character is progressively modified from a single to a double covalent bond.
These results point to the preliminary conclusion that the chemistry prior to detonation is strongly dependent on the pressure anisotropy, and starts before any band-gap closure, i.e., metallization which would be caused by bending of covalent bonds involving the nitro group.Accordingly, further decrease of the band gap seems to be favored by the formation of protons in the crystal.It remains to check this hypothesis by applying a more accurate, first-principles method.Because the formation of a free-electron gas is favored under very high pressure, it can be argued that LDA within DFT should suffice to produce reliable results for this regime.It is remarkable that our results, being derived under adiabatic compression and for zero temperature, are in qualitative agreement with ongoing first-principles, molecular-dynamics simulations at both high density and temperature. 71Specifically, these last simulations also predict that the C-H bond is highly stretched at an early stage of the compression.
It is of course premature to make any firm conjectures based on the present results alone.Yet, it is tempting to place our discussion in the context of detonation experiments for liquid nitromethane. 68 -70It has not escaped our attention that the prediction of a C-H bond high-stretching here suggests that a proton abstraction may be indeed involved in the first activation step, thus reducing the observed time to selfexplosion under high pressure. 68This type of proton dissociation has a long history in high pressure studies.While our prediction is not rigorous, the reasoning that implies it warrants some attention; we believe that more theoretical work is needed in this direction.
It is expected that our results for the band-gap behavior will be modified with inclusion of dynamical effects.Molecular-dynamics simulations are only suggestive at the moment. 35,71Although restricted to the ground state, the present analysis has indicated properties in agreement with more expensive, first-principles calculations that are typically limited to smaller systems. 35It appears therefore promising to apply the SCC-DFTB method to more complex energetic materials such as RDX and PETN.
postulated that the initiation a͒ On leave from the Department of Mathematics, Massachusetts Institute of Technology, Cambridge, MA 02139.b͒ Also at: Department of Physics, Harvard University, Cambridge, MA 02138.c͒ Also at: Universita ¨t-GH, Paderborn, Fachbereich Physik, Theoretische Physik, D-33098 Paderborn, Germany, and Department of Molecular Biophysics, German Cancer Research Center, D-69120 Heidelberg, Germany.

FIG. 1 .
FIG. 1. Two different configurations of a nitromethane molecule.͑a͒ The molecule with the atoms fully relaxed; C, N, and O all lie in the same plane (␥ϭ0).͑b͒ The molecule with the NO 2 group bent out of the original CNOO plane of the relaxed structure (␥Ϸ53°); the structure is taken to be very close to the one given in Ref. 36, which corresponds to a singlet-triplet minimum energy crossing point found by using a first-principles method.

FIG. 3 .
FIG. 3. ͑a͒ Unit cell of the nitromethane molecular crystal.Atoms are relaxed and the nitro group lies in the CNOO plane in each molecule.͑b͒ HOMO-LUMO energy gap ͑in eV͒ of the perfect crystal under uniform and uniaxial compression of the primitive unit cell.

FIG. 4 .
FIG. 4. HOMO-LUMO gap change for the nitromethane crystal under uniform compression down to 30% of the original volume.The gap values are shown relative to the initial configuration of the ͑unstrained͒ unit cell.The sudden drop of the gap to lower values corresponds to the reported C-H bond high stretching.