Classical and Quantum Mechanical Studies of Crystalline Ammonium

We have developed an intermolecular potential that describes the structure of the ADN crystal in the approximation of rigid ions. This potential is co...
0 downloads 0 Views 135KB Size
Download PDF
6774

J. Phys. Chem. B 1999, 103, 6774-6782

Classical and Quantum Mechanical Studies of Crystalline Ammonium Dinitramide Dan C. Sorescu† and Donald L. Thompson* Department of Chemistry, Oklahoma State UniVersity, Stillwater, Oklahoma 74078 ReceiVed: April 5, 1999; In Final Form: June 8, 1999

Plane-wave ab initio calculations based on density functional theory and the pseudopotential method have been used to investigate the structural properties of crystalline ammonium dinitramide (ADN). The optimization of the crystal structure has been done with full relaxation of atomic positions and lattice parameters under P21/c symmetry. The calculations were performed using periodic boundary conditions in all three directions. The predicted crystal structure is in good agreement with X-ray data and indicates no internal symmetry for the dinitramide ion, in contradistinction to the gas phase theoretical results. The two halves of the dinitramine ion are twisted with respect to each other while the NO2 groups are rotated out of the N-N-N plane. We have developed an intermolecular potential that describes the structure of the ADN crystal in the approximation of rigid ions. This potential is composed of pairwise Lennard-Jones, hydrogen-bonding terms, and Coulombic interactions. Crystal-packing calculations without symmetry constraints performed with the potential reproduce accurately the main crystallographic features and yield very good agreement with the estimated lattice energy. This potential was further tested in isothermal-isobaric molecular dynamics simulations at atmospheric pressure and in the temperature range 4.2-350 K. The main temperature effect is an increase in rotational disorder of the ammonium ions, without any significant change of the translational order of the ions. The thermal expansion coefficients calculated for the model indicate anisotropic behavior.

I. Introduction The behavior of crystalline ammonium dinitramide (ADN), NH4N(NO2)2, under thermal conditions has been the subject of a number of experimental studies in recent years.1-3 The main interest for this new type of stabilized energetic inorganic salt compound4 is mainly due to its potential use as an oxidizing component of propellant formulations. The use of ADN would eliminate the environmental problems resulting from the formation of HCl in the combustion of ammonium perchlorate, which is currently used in solid rocket propellants. Moreover, the new oxy anion of nitrogen, the dinitramide anion N3O4-, can be prepared in combination with a variety of different counterions, including ammonium, lithium, potassium, and cesium,5 leading to materials with high oxygen density. The structure of the dinitramide ion (DN-) has been the subject of several theoretical studies.6-8 It was found that in the gas phase the equilibrium structure of N(NO2)2- has C2 symmetry. However, in solution and solid phases the oxy anion of nitrogen was found generally to have no internal symmetry due to the steric and electronic properties of its environment, except in the cases of the lithium salt, in which DN- has C2 symmetry, and the dinitroazetidinium salt, in which DN- has mirror symmetry.5 Another difference between ADN and the lithium, potassium, and cesium salts is the type of bonding between the cations and anions. It was found5 that extensive hydrogen bonding involving all four ammonium hydrogen atoms are present in ADN, while in the alkali salts the nature of the cation-anion interactions are predominantly electrostatic. Thus, detailed investigations of the influence of different counterions on the properties of DN- are of primary fundamental importance. † Current-mailing address: Department of Chemistry, University of Pittsburgh, Pittsburgh, PA 15260.

In the present study we focus on the ADN crystal. The type of intermolecular bonding in the ADN crystal is a determining factor of its density and consequently of its possible uses as an energetic material. Two polymorphs of crystalline ADN, R and β, have been observed, both of which undergo irreversible thermal decomposition at the melting point.9 The R-ADN phase corresponding to low-pressure regime was resolved by X-ray diffraction measurements.5 The unit cell is monoclinic with four molecules ([NH4]+[N3O4]-) per unit cell (see Figure 1). This phase is stable up to 2.0 GPa over a large range of temperatures. This phase melts above 92 °C. At about 2.0 GPa the R phase transforms to β, which also has monoclinic symmetry. This phase is stable above 2.0 GPa in the temperature range from -75 to 140 °C. Due to the fact that only the R phase has been characterized experimentally, we have considered this phase in the present study. Our first objective was to investigate the structure and electronic properties of R-ADN using first-principles calculations. In previous studies6-8 only the gas-phase structure of the DN ion was considered. In the present study we have performed calculations for DN ions in the presence of ammonium ions by simulating the bulk ADN crystal using density functional theory (DFT) and the pseudopotential method with periodic boundary conditions in all three directions. The equilibrium structure of the crystal is obtained by full relaxation of the unit cell parameters allowed by the P21/c crystal symmetry as well as of the ionic positions inside the unit cell. Ideally, the equilibrium and dynamical properties of molecular crystals would be obtainable from ab initio molecular dynamics simulations. However, they are computationally extremely demanding. The alternative is to develop accurate intermolecular potentials that can be used in classical molecular dynamics (MD) or Monte Carlo (MC) simulations. We have successfully developed such intermolecular potentials for a series of nitramine

10.1021/jp9911447 CCC: $18.00 © 1999 American Chemical Society Published on Web 07/20/1999

Studies of Crystalline Ammonium Dinitramide

J. Phys. Chem. B, Vol. 103, No. 32, 1999 6775 ensemble are presented in section III. In section IV we give a summary of the main conclusions. II. Computational Methods

Figure 1. Representation of the R-ADN crystal unit cell with monoclinic space group P21/c and Z ) 4 molecules per unit cell.

crystals10a-c including some of the most important explosives, such as hexahydro-1,3,5-trinitro-1,3,5-s-triazine (RDX), 1,3,5,7tetranitro-1,3,5,7-tetraazacyclooctane (HMX), and 2,4,6,8,10hexanitrohexaazaisowurtzitane (HNIW), and have analyzed the transferability of such potentials to a set of 30 nitramine crystals10d and 50 nonnitramine compounds.10e We report here the development of an intermolecular potential that describes the structure of the R-ADN crystal. This potential will allow direct theoretical investigations of the response of ADN crystals to different external stimuli such as heating or pressure, particularly for conditions that are close to normal pressures and temperatures. In addition, the present potential model can be used to theoretically design crystals with different types of ionic packing. The methodology followed to develop the ADN intermolecular potential is similar to that used previously for the RDX10a and 5-nitro-2,4-dihydro-3H-1,2,4-triazol-3-one (NTO) crystals.11 The intermolecular interactions are assumed to be describable by simple isotropic potentials, such as the (6-exp) Buckingham potentials or the (6-m) Lennard-Jones potentials, with explicit consideration of the electrostatic interactions between the charges associated with the various atoms of different molecules. We find that for R-ADN, the crystal structure and its energy can be reasonably predicted by using a superposition of Coulombic terms, 6-12 Lennard-Jones potentials, and 10-12 potentials to describe the hydrogen bonding. The Coulombic terms were determined by fitting partial charges centered on each atom of the DN and NH4 ions to a quantum mechanically derived electrostatic potential. The parametrization of the potential function is done such that the model reproduces, in molecular packing (MP) calculations without symmetry constraints, the experimental structure of the crystal and its lattice energy. In all these calculations, the ions are assumed to be rigid, and the structure is described by the center of mass position and the orientational parameters (the Euler angles or the set of quaternions) for each of the ions in the unit cell. The intermolecular potential determined in MP calculations was further tested in isothermal-isobaric molecular dynamics (NPTMD) simulations as a function of temperature over the range 4.2-350 K at atmospheric pressure. The organization of the paper is as follows: In section II we describe the computational methods. The results of total energy calculations, lattice energy minimization by molecular packing calculations without symmetry constraints, and the results of trajectories calculations in the constant temperature and pressure

A. Total Energy Calculations. The calculations were performed using the commercial version of the software package CASTEP (Cambridge Serial Total Energy Package).12 This program evaluates the total energy of periodically repeating geometries based on DFT and the pseudopotential approximation. In the latter, only the valence electrons are represented explicitly in the calculations, the valence-core interactions being described by nonlocal pseudopotentials. Periodic boundary conditions are used, with the occupied electronic orbitals expanded in a plane-wave basis. The expansion includes all plane waves with kinetic energies, p2k2/2m < Ecut, where k is the wave vector, m the electronic mass, and Ecut is the chosen cutoff energy. The value of the cutoff energy is chosen to ensure convergence with respect to the basis set. A gradient-corrected form of the exchange correlation functional (GGA) was used in the manner suggested by White and Bird.13 It has been shown14 that the GGA approximation gives more accurate results for the energetics of molecular interactions than the local density approximation (LDA). The Brillouin zone was sampled with the Monkhorst-Pack scheme.15 The pseudopotentials used in this study are norm-conserving of the form suggested by Kleinman and Bylander16 and optimized using the scheme of Lin et al.17 The self-consistent ground state of the system was determined using a band-byband conjugate gradient technique to minimize the total energy of the system with respect to the plane-wave coefficients. The optimization of different atomic configurations reported in this study was done using the Broyden-Fletcher-Goldfarb-Shanno (BFGS) minimization technique, with the following thresholds for the converged structure: energy change per atom less than 1.0 × 10-6 eV, residual force less than 0.02 eV/Å, the displacement of atoms during the geometry optimization less than 0.001 Å, and the residual bulk stress less than 0.1 GPa. B. Molecular Packing Calculations. A general procedure for testing empirical or semiempirical intermolecular potential energy functions for organic crystals is the use of molecular packing calculations.18 The basic idea is to minimize the lattice energy with respect to the structural degrees of freedom of the crystal. For a crystal with Z rigid molecules per unit cell, the degrees of freedom are determined by the positions and orientations of the molecules in the unit cell as well as the dimensions and angles of the unit cell. Significant reduction of the computational time necessary to minimize the lattice energy starting with a trial configuration can be obtained by imposing and maintaining the observed space group symmetry of the crystal throughout energy minimization. However, studies done by Hagler and co-workers19 and more recently by Gibson and Scheraga20 show that maintenance of the observed space group symmetry, after energy minimization, should be considered a necessary condition in assessing the accuracy of the empirical potential energy. We have reached a similar conclusion in our studies of RDX crystals10a and of several nonnitramine crystals;10e that is, for an accurate force field for a crystal there are no significant differences when the minimization of the lattice energy is done using space group symmetry constraints and when these constraints are removed. We have chosen to use the method of crystal packing calculations without symmetry constraints. This can be done using the algorithm proposed by Gibson and Scheraga20 for efficient minimization of the energy of a fully variable lattice

6776 J. Phys. Chem. B, Vol. 103, No. 32, 1999

Sorescu and Thompson

composed by rigid molecules and implemented in the LMIN21 program. This algorithm makes use of Gay’s22 secant-type unconstrained minimization solver (SUMSL routine) to minimize the lattice energy. The gradients of the energy with respect to generalized coordinates, i.e., the 6Z rigid-body parameters of the Z molecules in the unit cell and the six lattice parameters, are evaluated analytically. The potential function used to describe the ADN intermolecular interactions was constructed as a sum of pairwise additive Lennard-Jones (LJ), hydrogen bonding (HB), and Coulombic (C) potentials of the forms:

[( ) ( ) ] [( ) ( ) ]

0 VLJ ij (r) ) ij

0 VHB ij (r) ) ij 5

r0ij 12 r

r0ij r

-2

12

-6

r0ij 6 r

r0ij r

(1)

10

(2)

and

Vcij(r) )

qqii 4π0r

(3)

where 0ij is the energy minimum for the pair of atoms i and j, r0ij is the interatomic distance at the energy minimum, qi and qj are the electrostatic charges on the atoms, and 0 is the dielectric permittivity constant for a vacuum. On the basis of eqs 1-3 the lattice energy can be evaluated as the sum of interactions between the ions in the central unit and all the other ions in the image cells. For practical computational purposes, the lattice sums are evaluated over cells within a sphere centered at the origin and of radius limited by the chosen cutoff distance. To ensure the continuity of the energy and its first derivative, both the LJ and HB interactions are cut off at a distance Qr0 using a cubic feather (spline) applied over the distance Pr0-Qr0 where the P and Q parameters specify, respectively, the start and the end of the feather and represent the multiplicative coefficients for the interatomic distance at the energy minimum r0. In all our calculations the parameter P was set equal to Q - 0.5. The influence of value of the cutoff distance on the final lattice parameters is discussed in the following sections. The evaluation of the Coulombic sums over the infinite periodic lattice was done using the standard Ewald’s transformation technique20 to improve the slow rate of convergence of these sums. Further details of the summations procedures of different energy terms over the lattice as well as of the adjustable parameter which determines the relative contribution of the direct and reciprocal-space terms in Ewald sum are given in ref 20. As the purpose of our calculations was to determine the best values of the potential parameters to describe the R-ADN crystal, the experimentally observed geometry5,23 was taken as the starting point for energy minimization. The orientation of the ions with respect to the Cartesian space-fixed system was described by setting the origin of the body-fixed coordinates at the centroids of the ions and with the axes aligned along the principal axes of inertia, with unit weight assigned to each atom.20 During energy minimization, the values of the translational vectors, describing the origin of the centroids and of the Euler angles, representing the orientations of the ions, give the necessary information on the relative positions of the ions in the unit cell.

TABLE 1: ADN Atom-Atom Intermolecular Potential Parameters pair (R-β) Lennard-Jones H-H N-N O-O N-O N-H hydrogen bonding N-H‚‚‚O-N

0 Rβ (kcal/mol)

0 rRβ (Å)

0.424 160 0.071 390 0.227 390 0.105 747 0.122 860

1.81 3.69 3.36 3.805 2.76

0.118 571 6

2.72

To determine the space group symmetry of the crystal predicted by the semiempirical potential, at the beginning and at the end of the lattice-energy minimization, the symmetry operations that transform the ions in the unit cell are computed. The space group is considered to be conserved if the symmetry operations, as defined in the International Tables of Crystallography,24 between a given ion and the remaining ions in the unit cell remain unchanged and if the lattice parameters fixed by the lattice symmetry have not been modified significantly. If the space symmetry is not conserved as a result of energy minimization, a new space group is deduced on the basis of the final symmetry relations. It has been shown20 that for the cases when the space group symmetry remains unchanged, some of the symmetry transformations may be lost during energy minimization, but they are regained before convergence takes place. In the particular case of R-ADN crystal, with space group P21/c, the angles R and γ of the unit cell were found to remain very close to 90°, as required by space group symmetry. As mentioned before, the optimum set of potential parameter values should minimize the difference between the predicted and experimental geometries of the crystal. In addition, we imposed the condition that the calculated energy matches the static lattice energy U0 of the crystal. We have taken the value of -602.5 kcal/mol determined recently by Politzer et al.25 using a density functional procedure for the lattice energy. The assignment of the atom-centered monopole charges has been made using the set that best reproduces the quantum mechanically derived electrostatic potential calculated over grid points surrounding the van der Waals surface of the DN and NH4 ions. This method of fitting to the electrostatic potential was proposed by Breneman and Wiberg26 and is incorporated in the Gaussian 94 package of programs27 under the keyword CHELPG. We have used this method in conjunction with Møller-Plesset perturbation theory at the MP2/6-311+G** level.28 The atomic arrangement of the ions used in these calculations is consistent with the crystallographic configuration.23 The full list of potential parameters used in our calculations is given in Table 1. C. Constant Pressure and Temperature Molecular Dynamics Calculations. A more comprehensive test of the intermolecular potential for ADN crystal was done in a constant pressure and temperature (NPT) molecular dynamics simulation, in which there are no geometric constraints other than the assumption of rigid-body ions. This method yields average equilibrium properties of the lattice as functions of temperature and pressure. We have used the Nose´-Hoover barostat algorithm29 as implemented in the program DL•POLY•2.0,30 to simulate the ADN crystal at various temperatures in the range 4.2-350 K and at atmospheric pressure. The equations of motion for both the translation of rigid ions and the simulation cell are integrated using the Verlet leapfrog scheme.31 The molecular rotational motion is handled using Fincham’s implicit quaternion algorithm.32

Studies of Crystalline Ammonium Dinitramide The MD simulation cell consists of a box containing 24 crystallographic unit cells (a 3 × 2 × 4 box of unit cells). In the initial simulation, corresponding to the lowest temperature, the position and orientation of the ions in the unit cell were taken to be identical to those for the experimental structure. In all production runs done in the Nose´-Hoover NPT ensemble, the system was integrated for 12 000 time steps (1 time step ) 2 × 10-15 s), of which 2000 steps were for equilibration. In the equilibration period the velocities were scaled after every five steps so that the internal temperature of the crystal mimicked the imposed external temperature. Then, average properties were calculated over the next 10 000 integration steps. In subsequent runs, performed at successively higher temperatures, the initial configurations of the molecular positions and velocities were taken from the previous simulation at the end of the production run. The lattice sums were calculated subject to the use of minimum-image periodic boundary conditions in all dimensions.31 The interactions were determined between the sites (atoms) in the simulation box and the nearest image sites within the cutoff distance Rcut ) 9.1 Å. The Coulombic long-range interaction were evaluated by using Ewald’s method.31 Several types of calculations were performed to obtain information about structural characteristics of the crystal. First, the mean positions of the ions in the simulation box were calculated. This was done by using a folding procedure similar to that used in MDCSPC4B program,33 in which the ions considered in their average positions in different subcells of the MD box are superimposed over the primitive cell. The content of these subcells is considered identical if these positions coincide within the root-mean-square (rms) deviations of their averaged values. Additional information about the structure of the crystal has been obtained by calculating the center-of mass (COM) and the site-site radial distribution functions (RDF). Finally, the average orientations of the ions in the unit cell were determined by evaluating the averages of the Euler angles of the principal molecular axis. The above quantities, that is, the mean positions, RDFs, and Euler angles, were calculated from recordings done at every 10th step during the trajectory integrations. Results and Discussion A. Ab Initio Total Energy Calculations. We investigated the bulk structure of the R-ADN crystal by performing a geometry optimization of both the unit cell parameters and ionic positions under the experimentally determined5,23 space group symmetry P21/c. The relaxation of the unit cell was done with respect to the independent lattice parameters a, b, and c and with respect to the crystallographic angle β while the angles R and γ were fixed at 90.0°. In these calculations, a number of 2k-points were used for Brillouin zone sampling, which were generated using the Monkhorst-Pack scheme15 with mesh parameters 2 × 1 × 2 along the three reciprocal lattice vectors. All calculated results are given in Table 2. An analysis of the results given in Table 2 shows that the agreement between the calculated and the experimental crystallographic lattice parameters is very good. There is also a very small dependence of the calculated values on the cutoff energy Ecut. The errors relative to experimental values for cell dimensions a, b, and c are, respectively, 1.62%, 1.34%, and 0.59% at a cutoff energy of Ecut ) 1000 eV and 1.69%, 1.66%. and 0.64% at Ecut ) 1400 eV, which indicate that the lattice parameters are practically converged for these cutoff energies. The corresponding errors for the β angle of the unit cell are, respectively,

J. Phys. Chem. B, Vol. 103, No. 32, 1999 6777 TABLE 2: Ab Initio Predicted Lattice Parameters for Bulk ADN Crystal at Different Cutoff Energiesa Ecut (eV)

a (Å)

1000.0

7.0264 (1.62%) 7.0298 (1.67%) 7.0308 (1.69%)

1200.0 1400.0

b (Å)

c (Å)

R (deg)

β (deg)

γ (deg)

11.9453 5.6478 90.000 99.694 90.000 (1.34%) (0.59%) (-0.70%) 11.9808 5.6472 90.000 99.702 90.000 (1.65%) (0.58%) (-0.70%) 11.9860 5.6504 90.000 99.788 90.000 (1.66%) (0.64%) (-0.61%)

a The values in parenthesis represent the percentage difference from the experimental values (ref 23).

-0.70% and -0.61% at these cutoff energies, while the other two angles R and γ remain equal to the experimental value of 90.0°, consistent with space group symmetry P21/c. Our analyses show that when the P21/c symmetry constraints are removed, the final optimized structure does not change significantly from the above values. We calculated the self-consistent band structure for the theoretical structure obtained at Ecut ) 1200 eV. The band gap at Γ point is 3.18 eV. No experimental data are available for comparison. The calculated geometries of the dinitramide and ammonium ions are compared with experimentally observed data23 in Table 3. For comparison we have included the results of previous calculations using nonlocal DFT6 and MP2.34 The labeling scheme in Table 3 is depicted in Figure 1. It can be seen by analyzing the data in Table 3 that there is a generally good agreement between the calculated and the experimental structures of N(NO2)2-. In particular there is the predicted nonequivalence of N1-N2 and N1-N3 bond lengths and of the N-N-O angles of each nitro group. In addition both N2-O5 and N3-O6 bond lengths are shorter than the corresponding N2O4 and N3-O7 bonds, in agreement with the crystallographic observations. The analysis of the torsional angles in Table 3 indicates that the two NO2 groups are not in the same plane. These groups are rotated out of the N-N-N plane to minimize the interaction between O5 and O6 by increasing their separation. This rotation might also be due to the increased conjugation between the NO2 groups and the lone pair on the central nitrogen atom. Similar findings were obtained using a nonlocal density functional procedure by Politzer et al.,6 but as can be seen in Table 3 the bond lengths predicted at this level have slightly larger errors than the experimental values. In contradistinction with geometries predicted in this study and by the nonlocal density functional procedure,6 the calculated minimum-energy structure of N(NO2)2- at the MP2/6-31+G* level of theory was found to possess C2 symmetry35 with the oxygen atoms above and below the plane of the three nitrogens. However, it was argued that the potential energy surface for rotation of the nitro groups is very shallow with rotational barriers of the NO2 groups of less than 3 kcal/mol.35 In Table 3 we also give the calculated values of N-H bonds in the ammonium ion as well as the corresponding hydrogenbonding distances. All of the H atoms of the NH4+ are involved in hydrogen bonds with the oxygen atoms of the surrounding dinitramide ions. The corresponding N-H‚‚‚O bond lengths vary between 1.833 and 2.145 Å. These values are smaller than the corresponding experimental values for the range 2.1432.281 Å. The main reason for such differences is the uncertainty in the H atom positions obtained using X-ray diffraction techniques. It can also be seen from the data in Table 3 that N-H bonds lengths in NH4+ are not equivalent. The longer N8-H9 and N8-H10 bonds correspond to the shorter intermolecular hydrogen-bonding distances N8-H9‚‚‚O4 and

6778 J. Phys. Chem. B, Vol. 103, No. 32, 1999

Sorescu and Thompson

TABLE 3: Ab Initio Calculated and Observed Geometries for N(NO2)2- and NH4+ Ions in ADN Crystala parameter

this work

N(1)-N(2) N(1)-N(3) N(2)-O(4) N(2)-O(5) N(3)-O(6) N(3)-O(7) N(2)-N(1)-N(3) N(1)-N(2)-O(4) N(1)-N(2)-O(5) N(1)-N(3)-O(6) N(1)-N(3)-O(7) O(4)-N(2)-O(5) O(5)-N(3)-O(6) N(8)-H(9) N(8)-H(10) N(8)-H(11) N(8)-H(12)

1.363 (-0.95) 1.359 (-0.07) 1.215 (-1.93) 1.200 (-2.12) 1.191 (-2.69) 1.225 (-2.07) 115.0 (1.59) 113.2 (0.35) 123.1 (-0.24) 124.1 (-0.79) 112.8 (0.17) 123.6 (0.24) 122.9 (0.65) 1.039 1.028 1.026 1.027

N(3)-N(1)-N(2)-O(4) N(3)-N(1)-N(2)-O(5) N(2)-N(1)-N(3)-O(6) N(2)-N(1)-N(3)-O(7)

158.6 (0.81) -26.0 (-7.80) -24.6 (2.92) 160.8 (-0.80)

NLDF/GGA/b

MP2/6-31+G* c

observedd

1.407 (2.25) 1.402 (3.08) 1.270 (2.58) 1.269 (3.50) 1.259 (2.85) 1.262 (1.37) 112.4 (-0.70) 113.4 (-0.35) 123.3 (-0.08) 123.3 (-1.43) 113.5 (0.80) 123.0 (-0.25) 123.0 (0.73)

1.3778 (0.13) 1.3778 (0.13) 1.2523 (1.15) 1.2575 (2.56) 1.2575 (2.73) 1.2523 (0.10) 112.79 (-0.36) 113.88 (0.78) 122.67 (-0.59) 122.67 (-1.94) 113.88 (1.14) 123.14 (-0.21) 123.14 (0.85)

1.376 1.360 1.238 1.226 1.224 1.251 113.2 113.0 123.4 125.1 112.6 123.4 122.1 0.834 0.852 0.834 0.878

149.93 (-4.68) -23.82 (-15.5) -23.82 (-0.33) 149.93 (-7.5)

157.3 -28.2 -23.9 162.1

Torsion Angles

N(2) nitro group N(3) nitro group N(8)-H(9)‚‚‚O(4) N(8)-H(10)‚‚‚O(7) N(8)-H(11)‚‚‚O(7) N(8)-H(12)‚‚‚O(5) O(4)‚‚‚O(6)

-26.0 (-7.80) -27.8 (16.31)

Twist and Bend Anglesd 23.6 (-6.71) 4.3 (-15.6) 21.8 (4.80) 4.9 (-7.54) 1.833 2.029 2.149 2.100 2.588 (-0.19)

25.3 5.1 20.8 5.3 2.143 2.172 2.235 2.281 2.593

a The values in parenthesis indicate the percentage difference from experimental values. b Data from ref 6. c Data from ref 34. d Data from refs 5 and 23. The twist angles were calculated as the angles between the plane NO2 group and the plane of N(2)-N(1)-N(3) atoms. The bend angle was determined as the angle between the O-N(i)-O plane and N(i)-N(1) axis, with i ) 2, 3.

N8-H10‚‚‚O7. We can conclude that the geometrical parameters of both DN and ammonium ions are strongly influenced by their environment and by the presence of hydrogen-bonding interactions between these ions. Two descriptions have been previously proposed by Politzer et al.6 for the bonding mechanism in the dinitramide anion. They suggested for the first mechanism that if the dinitramide is twisted at an N-N-N angle of 112.4° (their ab initio calculated value), the electron pairs around the central nitrogen is tetrahedrally distributed and the conjugation optimized by rotation of the NO2 planes such that they are approximately perpendicular to the lone-pair orbitals. The second mechanism is based on planar dinitramide. The lone pairs are in different orbitals, an sp2 orbital in the N-N-N plane and a p-orbital perpendicular to it. Since both the experimental N-N bonds and the N-N-N angle are intermediate between the ideal tetrahedral and trigonal planar geometries it was also concluded5 that the best description of the bridging nitrogen is between these two extremes. This situation is characterized by the delocalization of all lone pairs on the bridging nitrogen and by allowing conjugation with both NO2 groups. This last interpretation of the bonding mechanism is supported by our calculations. First, the calculated N-N and N-N-N geometrical parameters are very close to the corresponding experimental results. Second, as can be seen in Figure 2a, the HOMO orbital for dinitramide ion is highly asymmetrical and with significant contributions along the N-N bonds. In Figure 2b we present a contour plot of the electron density along a plane perpendicular to the b crystal axis containing a section through the H10 site of the NH4+ ion. As can be seen, the valence electron density is highly concentrated around individual N and O sites of N(NO2)- and the N atom of NH4+.

Figure 2. (a) Configuration of the HOMO as determined from ab initio calculations. (b) Charge density distribution map along a plane perpendicular to the b axis which passes through the H10 atom of the ammonium ion.

In addition, the polarization of charge density is along the N8H10‚‚‚O7 hydrogen bond. B. Molecular Packing Calculations. The results of the molecular packing calculations obtained using the potential parameters given in Table 1 are summarized in Table 4 as functions of the cutoff distance parameter Q. Since evaluation

Studies of Crystalline Ammonium Dinitramide

J. Phys. Chem. B, Vol. 103, No. 32, 1999 6779

TABLE 4: Lattice Parameters and Energies Obtained in Crystal Packing Calculations without Symmetry Constraints lattice energya Q

c

EXPd LMIN 5.5 10.5 15.5 20.5 25.5

NB

lattice parametersb ES

-602.496e -47.1255 -47.7271 -47.7933 -47.8103 -47.8150

25.5

-561.1697 -561.2220 -561.2220 -561.2220 -561.2220

b

c

R

β

6.9138

11.7867

5.6143

90.000

100.405

90.000

6.8274 6.8267 6.8262 6.8264 6.8264 (-1.26)h

12.0078 12.0057 12.0061 12.0056 12.0057 (1.86)

5.6473 5.6469 5.6471 5.6470 5.6470 (0.58)

90.007 90.008 90.004 90.007 90.009 (0.01)

101.180 101.183 101.187 101.185 101.185 (0.78)

89.990 89.982 89.991 89.987 89.981 (-0.02)

a

γ

∆xf

∆y

∆z

∆Φ(deg)g

∆Θ (deg)

∆Ψ (deg)

-0.012 -0.002

-0.004 -0.005

-0.019 -0.002

-1.03 5.29

0.70 5.39

3.57 -1.57

a Nonbonded (NB) and electrostatic (ES) lattice energies per molecule in kJ/mol. b Lattice dimensions a, b, c in angstroms and angles R, β, and γ in degrees. c Cutoff parameter as described in text. d Experimental values from ref 23. e Calculated lattice energy in kJ/mol from ref 25. f Change in fractional coordinates of molecular centroids for DN and ammonium ions. g Change in Euler angles of molecular centroids for DN and ammonium ions. h The percentage change in lattice and molecular parameters after energy minimization is determined as a function of the experimental geometry (ref 23) and is given in parentheses.

TABLE 5: Predicted Lattice Parameters as Functions of Temperature and the Thermal Expansion Coefficients (χ) at T ) 300 K T (K) 4.2 100.0 200.0 273.1 300.0 350.0 χa a

a (Å)

b (Å)

c (Å)

R (deg)

β (deg)

γ (deg)

volume (Å3)

6.8263 (0.0166 6.8141 (0.0273 6.8222 (0.0204 6.8305 (0.0303 6.8338 (0.0357 6.8400 (0.0334 1.52

12.0094 (0.0385 12.0857 (0.0447 12.1970 (0.0529 12.2752 (0.0558 12.3086 (0.0602 12.3654 (0.0669 8.69

5.6493 (0.0122 5.7032 (0.0219 5.7413 (0.0226 5.7661 (0.0293 5.7760 (0.0351 5.7919 (0.0365 7.15

89.998 (0.074 90.024 (0.226 90.030 (0.263 89.995 (0.335 90.011 (0.346 90.025 (0.406

101.197 (0.100 100.964 (0.182 100.954 (0.168 100.980 (0.335 101.002 (0.303 101.025 (0.373

90.00 (0.05 90.031 (0.317 89.977 (0.188 90.019 (0.209 89.960 (0.337 90.021 (0.303

454.333 (0.409 461.083 (0.989 469.041 (1.323 474.583 (1.822 476.916 (1.681 480.833 (2.090 16.96

Units: linear expansion coefficients are in 10-5 K-1; volume expansion coefficient 10-5 K-1.

of the dispersive terms of the form r-6 in the nonbonded potential was done without an accelerated convergence method,10a a slight dependence of the lattice energy on the cutoff distance is expected. The results in Table 4 show that this dependence is more significant for Q values smaller than 10, after which only small variations of the total energy per molecule and of the geometrical lattice and molecular parameters take place. For all the cases in Table 4, the symmetry space group was maintained as verified by the existent symmetry relations of the P21/c space group between the Z - 1 ions in the unit cell and the corresponding central ions with Z ) 1. The change in lattice and molecular parameters after energy minimization relative to the experimental structure5 are also given in parentheses in Table 4. The maximum deviations of the lattice dimensions and angles to experimental values are 1.86% and 0.78%, respectively, with small translations of the rigid molecules. There is essentially no change in the values of the angles R and γ, which remain close to 90°, consistent with space group symmetry. The Euler angles that described the molecular orientation in the unit cell also have small variations for the DN ion, but somewhat larger rotations takes place for ammonium ion. For the minimized lattice, the total energy is U ) -609.04 kJ/mol, with the major contribution given by the electrostatic interactions. The predicted lattice energy is in good agreement with the value -602.50 kJ/mol, which was determined by Politzer et al.25

These results suggest that accurate descriptions of the equilibrium properties of the R-ADN crystal can be predicted by this intermolecular potential. C. NPT Molecular Dynamics Calculations. A more realistic prediction of the structural lattice parameters can be obtained by considering the molecular motion as a function of temperature and pressure. For this purpose, we have determined the structural crystal parameters at atmospheric pressure and in the temperature range 4.2-350 K. These results are summarized in Table 5. As can be seen by comparing results in Tables 4 and 5, the lattice dimensions obtained at T ) 4.2 K are in very close agreement with those determined in the molecular packing calculations without symmetry constraints. This is expected since the thermal effects at 4.2 K should be minimal and the thermal averages at this temperature should be close to the values corresponding to the potential energy minimum. At 200 K the average lattice dimensions of these crystals agree satisfactorily with the experimental values measured at 223 K,5 the corresponding differences for the a, b, and c lattice dimensions being 1.3%, 3.4%, and 2.2%, respectively. Moreover, in both low and ambient temperature simulations, the unit cell angles R and γ remain approximately equal to 90°, while the β angle varies slightly around the value of 101° over the temperature range 4.2-350 K. In Figure 3 we compare the results of the average fractional coordinates of the four molecules in the unit cell with the

6780 J. Phys. Chem. B, Vol. 103, No. 32, 1999

Sorescu and Thompson

Figure 3. Comparisons of the time-averaged center-of-mass fractional positions with experimental results at 350 K. The values for each molecule in the unit cell were averaged over the 24 unit cells in the simulation box.

corresponding experimental data. It is evident that increasing the temperature from 4.2 to 300 K does not produce any significant displacement of the molecular center of mass (COM). Additional support for the small degree of translation of the molecules inside the unit cell with the temperature increase can be obtained from the COM-COM radial distribution functions (RDF) given in Figure 4. Indeed, the RDFs at these temperatures correspond to well-ordered structures with correlation at long distances. The positions of the major peaks do not change significantly and the main temperature effect is the broadening of the peaks with partial overlapping of some of them. Further insight into the crystallographic structure can be obtained by analyzing the site-site RDFs for the case of shorter contacts, i.e., the N-O‚‚‚H intermolecular bonds. These RDFs were calculated as averages over all molecules in the MD box but limited by the cutoff distance used in potential calculation. The results are represented in Figure 5. It can be seen that at the highest temperature of our simulations, T ) 350 K, the RDFs of different O‚‚‚H pairs are very similar with loss of peak individuality. This trend is significantly different from the behavior at 4.2 K, where different pairs O-Hi, (i ) 1, 2, 3, 4) have specific variations. This indicates that by increasing the temperature the individuality of the H atoms in N-O‚‚‚H bonds is lost. This suggests a large degree of rotational disorder of the ammonium ions with the increase of temperature. We have also determined the linear and volume expansions coefficients at 300 K for the ADN crystal using the corresponding temperature dependence of these parameters. These results are summarized in Table 5. It can be seen that there is anisotropic thermal behavior of the crystal with the largest expansion along the b axis. At present, no experimental data for the thermal coefficients are available for comparison.

Figure 4. Radial distribution functions for center-of-mass-center-of mass pairs as functions of temperature: (a) 4.2 K and (b) 350 K.

Summary and Conclusions The results of this study indicate that reliable geometrical parameters for the R-ADN crystal can be obtained by using plane-wave total energy calculations based on DFT theory and the pseudopotential method. In particular we were able to describe the main geometrical features of the DN and NH4 ions in the ADN crystal. It was found that the DN ion has inequivalent N-N and N-O bond distances and that nitro groups are twisted from the N1-N2-N3 plane. This is mainly due to minimization of the steric effects between the O5 and O6 atoms and maximization of the conjugation between the NO2 groups and the lone pairs on the central nitrogen atom. A strong influence of environment was observed for the H atoms of NH4 ions, which presents nonequivalent N-H bond lengths. We have developed an intermolecular potential for the R-form of the ADN crystal based on 6-12 Lennard Jones potentials, 10-12 hydrogen-bond terms, and Coulombic interactions. Electrostatic charges of the atoms in the DN and ammonium ions were determined from fits to ab initio electrostatic potentials calculated at the MP2/6-311+G** level. The values of intermolecular potential were optimized to minimize the difference between the theoretical and experimental5,23 structures and their corresponding energies. A prime test of the proposed potential was done in symmetryunconstrained molecular packing calculations. It was shown that the space group symmetry is maintained throughout energy minimization and that the crystallographic parameters are

Studies of Crystalline Ammonium Dinitramide

J. Phys. Chem. B, Vol. 103, No. 32, 1999 6781

Figure 5. Radial distribution functions for the N-O‚‚‚H pairs as functions of temperature at (a-d) 4.2 K and (e-h) 350 K.

reproduced with an acceptable degree of accuracy. Very good agreement was obtained for the predicted lattice energy and previous theoretical estimations.25 The temperature dependencies of the physical parameters of the lattice have been investigated by performing isothermalisobaric molecular dynamics simulations at zero pressure over

the temperature range 4.2-350 K. The results of these calculations show that at 200 K the model reproduces the experimental unit cell dimensions within 3.4%. Additionally, little translational disorder occurs in thermal, unconstrained trajectories. However, the increase of temperature leads to the increase of rotational disorder of the ammonium ions.

6782 J. Phys. Chem. B, Vol. 103, No. 32, 1999 The linear and volumic thermal expansion coefficients of the crystal were determined from the averages of lattice dimensions extracted from trajectory calculations. The results obtained indicate the anisotropic behavior of the crystal with expansions preferentially along the b and c axes. This model is useful for predicting nonreactive processes in the R-ADN crystal. Refinement of this model can be made by including the effects of intramolecular motions, particularly of low-frequency torsional motions of the nitro group and of the ring. In addition, the model can be extended to reproduce not only geometry and energy parameters but also spectroscopic data for the ADN lattice. Further developments of the model to include the intramolecular degrees of freedom are expected to facilitate full atomistic investigations of the dynamics of this energetic material. Acknowledgment. This work was supported by the Air Force Office of Scientific Research (Grant No. F49620-95-l0310). We would like to thank Dr. Richard D. Gilardi of the Naval Research Laboratory, Washington, DC, for making available to us the crystallographic parameters of ADN crystal. References and Notes (1) Brill, T. B.; Brush, P. J.; Patil, D. G. Combust. Flame 1993, 92, 178. (2) Wight, C. A.; Vyazovkin, S. J. Phys. Chem. A 1997, 101, 5653. (3) Oxley, J. C.; Smith, J. L.; Zheng, W.; Roger, E.; Coburn, M. D. J. Phys. Chem. A 1997, 101, 5646. (4) Bottaro, J. C.; Schmidt, R. J.; Renwell, P. E.; Ross, D. S. World Intellectual Property Organization, International Application Number PCT/ US91/04268, December 26, 1991. (5) Gilardi, R.; Flippen-Anderson, J.; George, C.; Butcher, R. J. J. Am. Chem. Soc. 1997, 119, 9411 and references herein. (6) Politzer, P.; Seminario, J. M.; Concha, M. C.; Redgern, P. C. J. Mol. Struct. (THEOCHEM) 1993, 287, 235. (7) Michels, H. H.; Montgomery, J. A., Jr. J. Phys. Chem. 1993, 97, 6602. (8) Mebel, A. M.; Lin, M. C.; Morokuma, K.; Melius, C. F. J. Phys. Chem. 1995, 99, 6842. (9) Russsell, T. P.; Piermarini, G. J.; Block, S.; Miller, P. J. J. Phys. Chem. 1996, 100, 3248. (10) (a) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. 1997, 101, 798. (b) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. B 1998, 102, 948. (c) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. B 1998, 102, 6692. (d) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. A 1998, 102, 8386. (e) Sorescu, D. C.; Rice, B. M.; Thompson, D. L. J. Phys. Chem. A 1999, 103, 989. (11) Sorescu, D. C.; Thompson, D. L.; J. Phys. Chem. B 1997, 101, 3605.

Sorescu and Thompson (12) Payne, M. C.; Allan, D. C.; Arias, T. A.; Johannopoulus, J. D. ReV. Mod. Phys. 1992, 64, 1045. (13) White, J. A.; Bird, D. M. Phys. ReV. 1994, B 50, 4954. (14) White, J. A.; Bird, D. M.; Payne, M. C.; Stich, I., Phys. ReV. Lett. 1994, 73, 1404. (15) Monkhorst, H. J.; Pack, J. D. Phys. ReV. 1976, B 13, 5188. (16) Kleinman, L.; Bylander, D. M. Phys. ReV. Lett. 1980, 45, 566. (17) Lin, J. S.; Qteish, A.; Payne, M. C.; Heine, V. Phys. ReV. 1993, B 47, 4174. (18) Pertsin, A. J.; Kitaigorodsky, A. I. The Atom-Atom Potential Method, Applications to Organic Molecular Solids; Springer-Verlag: Berlin, 1987. (19) Hagler, A. T.; Lifson, S.; Dauber, P. J. Am. Chem. Soc. 1979, 101, 5122. (b) Hagler, A. T.; Dauber, P.; Lifson, S. J. Am. Chem. Soc. 1979, 101, 5131. (20) Gibson, K. D.; Scheraga, H. A.J. Phys. Chem. 1995, 99, 3752. (21) Gibson, K. D.; Scheraga, H. A. LMIN: A Program for Crystal Packing, QCPE No. 664. (22) Gay, D. M. ACM Trans. Math. Software 1983, 9, 503. (23) Gilardi, R. Private communication, 1998. (24) International Tables for Crystallography; Hahn, T., Ed.; Reidel: Dordrecht, Boston, 1983. (25) Politzer, P.; Seminario, J. M.; Concha, M. C. J. Mol. Struct. (THEOCHEM) 1998, 427, 123. (26) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (27) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Patersson, G. A.; Montgomery, J. A.; Raghavachari, K.; AL-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanyakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andres, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; HeadGordon, M.; Gonzales, C.; Pople, J. A. Gaussian 94, Revision C.3; Gaussian, Inc.; Pittsuburgh PA, 1995. (28) Mo¨ller, C. M. S. Phys. ReV. 1934, 46, 618. Hehre, W. J.; Ditchfield, R.; Pople, J. A. J. Chem. Phys. 1972, 56, 2257. Hariharan, P. C.; Pople, J. A. Theor. Chim. Acta 1973, 28, 213. Gordon, M. S. Chem. Phys. Lett. 1980, 76, 163. (29) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol. Phys. 1993, 78, 533. (30) DL•POLY is a package of molecular simulation routines written by W. Smith and T. R. Forester, copyright The Council for the Central Laboratory of the Research Councils, Daresbury Laboratory at Daresbury, Nr. Warrington, 1996. (31) Allen, M. P.; Tindesley, D. J. Computer Simulation of Liquids; Oxford University Press: New York, 1989. (32) Fincham, D. Mol. Simul. 1992, 8, 165. (33) Smith, W. MDSCPC4B, A Program for Molecular Dynamics Simulations of Phase Changes, CCP5 Program Library (SERC), May, 1991. (34) Christe, K. O.; Wilson, W. W.; Petrie, M. A.; Michels, H. H.; Bottaro, J. C.; Gilardi, R. Inorg. Chem. 1996, 35, 5068. (35) Michels, H. H.; Montgomery, J. A., Jr. J. Phys. Chem. 1993, 97, 6602.