ARMY RESEARCH LABORATORY

Theoretical Studies of

Solid Nitromethane

by Dan C. Sorescu, Betsy M. Rice,

and Donald L. Thompson

^^ff^^^^^^ß^^^''^^^'^

ARL-RP-14 February 2001

A reprint from the Journal ofPhysical Chemistry B, vol. 104, no. 35, pp. 8406-8419.

20010316 069

Approved for public release; distribution is unlimited.
The findings in this report are not to be construed as an official

Department of the Army position unless so designated by other

authorized documents.

Citation ofmanufacturer's or trade names does not constitute an

official endorsement or approval of the use thereof.

Destroy this report when it is no longer needed. Do not return

it to the originator.
Army Research Laboratory

Aberdeen Proving Ground, MD 21005-5066

ARL-RP-14 February 2001

Theoretical Studies of Solid Nitromethane

Dan C. Sorescu

Oklahoma State University

Betsy M. Rice

Weapons and Materials Research Directorate, ARL

Donald L. Thompson

Oklahoma State University

A reprint from the Journal ofPhysical Chemistry B, vol. 104, no. 35, pp. 8406-8419.

Approved for public release; distribution is unlimited.
Abstract

A classical potential to simulate the dynamics of a nitromethane crystal as a function

of temperature and pressure is described. The intramolecular part of the potential was

taken as superposition of bond stretching, bond bending, and torsional angles terms.

These terms were parametrized on the basis of the geometric and spectroscopic

(vibrational frequencies and eigenvectors) data obtained using ab initio molecular orbital

calculations performed at the B3LYP/6-31G level on an isolated molecule. The

intermolecular potential used is of the Buckingham 6-exp form plus charge-charge

Coulombic interactions and has been previously developed by us (Sorescu, D. C; Rice,

B. M.; Thompson, D. L. J. Phys. Chem. 1997, B101, 798) to simulate crystals containing

nitramine molecules and several other classes of nitro compounds. The analyses

performed using constant pressure and temperature molecular dynamics simulations and

molecular packing calculations indicate that the proposed potential model is able to

reproduce accurately the changes of the structural crystallographic parameters as

functions of temperature or pressure for the entire range of values investigated. In

addition, the calculated bulk modulus of nitromethane was found in excellent agreement

with the corresponding experimental results. Moreover, it was determined that the

present potential predicts correctly an experimentally observed 45° change in methyl

group orientation in the high-pressure regime relative to the low-temperature

configuration. The analysis of the linear expansion coefficients and linear compression

data indicate anisotropic behavior for the unit cell edges.
Theoretical Studies of Solid Nitromethane

Dan C. Sorescu, Betsy M. Rice, and Donald L. Thompson

Department of Chemistry, Oklahoma State University, Stillwater,

Oklahoma 74078, and The U. S. Army Research Laboratory,

Aberdeen Proving Ground, Maryland 21005

The Journal of

Physical Chemistry W

Reprinted from

Volume 104, Number 35, Pages 8406-8419
8406 J. Phys. Chem. B 2000, 104, 8406-8419

Theoretical Studies of Solid Nitromethane

Dan C. Sorescu,^5 Betsy M. Rice,*'* and Donald L. Thompson*

Department of Chemistry, Oklahoma State University, Stillwater, Oklahoma 74078, and

The U. S. Army Research Laboratory, Aberdeen Proving Ground, Maryland 21005

Received: March 13, 2000; In Final Form: June 13, 2000

A classical potential to simulate the dynamics of a nitromethane crystal as a function of temperature and

pressure is described. The intramolecular part of the potential was taken as superposition of bond stretching,

bond bending, and torsional angles terms. These terms were parametrized on the basis of the geometric and

spectroscopic (vibrational frequencies and eigenvectors) data obtained using ab initio molecular orbital

calculations performed at the B3LYP/6-31G* level on an isolated molecule. The intermolecular potential

used is of the Buckingham 6-exp form plus charge—charge Coulombic interactions and has been previously

developed by us (Sorescu, D. C; Rice, B. M.; Thompson, D. L. J. Phys. Chem. 1997, B101, 798) to simulate

crystals containing nitramine molecules and several other classes of nitro compounds. The analyses performed

using constant pressure and temperature molecular dynamics simulations and molecular packing calculations

indicate that the proposed potential model is able to reproduce accurately the changes of the structural

crystallographic parameters as functions of temperature or pressure for the entire range of values investigated.

In addition, the calculated bulk modulus of nitromethane was found in excellent agreement with the

corresponding experimental results. Moreover, it was determined that the present potential predicts correctly

an experimentally observed 45° change in methyl group orientation in the high-pressure regime relative to

the low-temperature configuration. The analysis of the linear expansion coefficients and linear compression

data indicate anisotropic behavior for the unit cell edges.

I. Introduction

This is the seventh in a series of papers describing our

development and assessment of interaction potentials to be used

in the study of dynamic processes in energetic materials. Our

original intent was to develop a model to study nonreactive

processes in the nitramine explosive RDX (1,3,5-hexahydro1,3,5-s-triazine).1 The potential that was developed consisted

of atom—atom (6-exp) Buckingham potential terms plus electrostatic interactions. The Coulombic interactions were obtained

through fitting atom-centered partial charges to a quantummechanically-determined electrostatic potential for a single RDX

molecule whose structure corresponded to that in the crystal at

ambient conditions. The remaining Buckingham parameters

were adjusted to reproduce the experimental structure of the

RDX crystal at ambient conditions. We found that this interaction potential could also describe the geometric parameters and

lattice energies of different polymorphic phases of two other

nitramine crystals: the polycyclic nitramine 2,4,6,8,10,12-

hexanitrohexaazaisowurtzitane (HNIW, or CL-20)2and the

monocyclic nitramine octahydro-l,3,5,7-tetranitro-l,3,5,7-tetraazacyclooctane (HMX).3

Isothermal—isobaric molecular dynamics (NPT-MD) simulations for these crystals predicted cell

parameters within a few percent of experimental values, and

little translational or rotational disorder ofthe molecules. Further

investigations to explore the limits of transferability of this

interaction potential to other energetic molecular crystals were

undertaken through performing molecular packing (MP)

* Author to whom correspondence should be addressed.

t Oklahoma State University.

* The U.S. Army Research Laboratory, Aberdeen Proving Ground.

§ Current mailing address: Department of Chemistry, University of

Pittsburgh, Pittsburgh, PA 15260.

calculations for 30 nitramine crystals. These included several

types of mono- and polycyclic nitramines, particularly crystals

of importance in energetic materials.4For most of the crystals,

the predicted structural lattice parameters obtained from molecular packing (MP) calculations differ by less than 2% from

the experimental structures, with small rotations and practically

no translations of the molecules in the asymmetric unit cell. A

further assessment on the limits of transferability of the

interaction potential was accomplished through molecular

packing calculations of 51 crystals containing non-nitramine

molecules with functional groups common to energetic materials.5 MP calculations using this interaction potential reproduced

the crystal structures to within 5% of experiment for these 51

non-nitramine systems, including the explosives pentaerythritol

tetranitrate (PETN), nitromethane, 2,4,6-trinitrotoluene (TNT),

and several nitrocubane derivatives.

More recently we have analyzed the dynamics of the

important energetic crystals RDX, HMX, HNIW, and PETN

under hydrostatic compression conditions using NPT-MD

simulations and this intermolecular potential.6In that study we

found that the predicted lattice parameters for the RDX, HMX,

and HNTW crystals are in good agreement with experimental

values over the entire range of pressures investigated experimentally. For the PETN crystal, the calculated crystallographic

parameters were in acceptable agreement with experimental data

for pressures up to 5 GPa. For higher pressures, the disagreements of predictions and experiment were attributed to the

inadequacy of the rigid-body approximation used to simulate

floppy molecules such as PETN.

The validity of the rigid molecular approximation was

assumed in our previous studies.1-6 As indicated by our results,

this model has been very successful in its ability to describe

10.1021/jp000942q CCC: $19.00 © 2000 American Chemical Society

Published on Web 08/15/2000
Solid Nitromethane J. Phys. Chem-'B, Volfl04, No. 35, 2000 8407

HC6!

Figure 1. Representation of the nitromethane crystal unit cell with

orthorhombic space group P2A21 and Z = 4 molecules per unit cell.

Atom labels are consistent with the indices given in Table 1.

the equilibrium structures of a variety of organic molecular

crystals under ambient conditions and with moderate increases

in pressure and temperature. However, the physical and chemical

processes of energetic materials that are of most interest occur

within the regime of high pressures and temperatures, a regime

in which conformatiönal molecular changes become important.

Consequently, further developments of the interaction potential

are necessary to describe more realistically the intramolecular

motion, molecular deformations, and the energy flow inside

these crystals. For this purpose, in the present paper we eliminate

the previously used rigid-molecule approximation1""6and extend

the current intermolecular potential to include a full intramolecular potential for use in simulations of energetic materials.

Particularly, in this work we consider the prototypical explosive,

nitromethane.

Nitromethane was selected for our first attempt at developing

a fully flexible model of an energetic molecular crystal since

numerous experimental investigations of its properties under a

wide range of conditions have been performed, thus providing

significant data for use in assessing the model potential.7-14

Single-crystal X-ray diffraction and neutron powder diffraction

measurements indicate that the low-temperature structure of

crystalline nitromethane belongs to the orthorhombic space

group P2i2i2i with Z = 4 molecules in the unit cell (see Figure

l).7In this configuration the methyl group is in the staggered

symmetric position with respect to the C—NO2 plane. At

ambient pressure, the crystal symmetry was observed to remain

unchanged when the temperature was increased from 4.2 to 228

K.7'

9

In the temperature range 243-244 K an intermediate

crystal phase (mesophase) was recently observed,10 followed

by transition to the liquid state at 244.73 K.11 No crystallographic details have been provided about the intermediate

phase at 243 K.10

The internal rotation of the methyl group in the crystal has

received significant attention, since experimental evidence

indicates that this motion is almost completely governed by

intermolecular interactions. Analyses of the neutron powder

diffraction patterns7for deuterated nitromethane at temperatures

ranging from 25 to 125 K provided temperature-dependent

amplitudes of librations of the methyl group about the C-N

bond. Quasielastic neutron scattering spectra8recorded at

temperatures ranging from 50 to 150 K are consistent with

rotations of the methyl groups through 120° jumps rather than

the 60° jumps observed in the gas phase. Furthermore, it was

found8that in the crystalline phase the activation energy for

the internal rotation of the methyl group is about 234 cal/mol,

significantly higher than the corresponding gas-phase barrier

of 6 cal/mol.15 Inelastic neutron-scattering measurements14 of

tunneling splittings and of the lowest rotational states indicated

that the methyl group activation barrier is significantly smaller

than the energy of the second excited rotational state by ~170

cal/mol. These findings have been confirmed by infrared16 and

Raman10 measurements, but the corresponding activation energies reported in these studies are slightly higher, namely 401

and 576 cal/mol, respectively. Several mechanisms have been

proposed to contribute to the dynamics of methyl groups

including the enhanced hopping due to interactions with

phonons, band tunneling through thermally broadened excited

states, or changes of the barrier with temperature.14

The influence of pressure on the internal rotational motion

of the methyl group has also been examined. The crystal

structure for isotropically compressed nitromethane has been

determined through X-ray diffraction' at pressures ranging from

0.3 to 15.3 GPa, at room temperature.912 It was found that the

space group and the packing arrangement remain unchanged

relative to the low-temperature structure.7 However, the measurements indicate a continuous change in the orientation of

the methyl group with the increase of pressure. For pressures

below 0.6 GPa the methyl group is freely rotating. At intermediate pressures (between 0.6 and 3.5 GPa) the rotation of the

methyl group is hindered, while for pressures above 3.5 GPa

the orientation of the methyl group becomes fixed. Moreover,

it was determined that at pressures above 3.5 GPa, the methyl

group is rotated by 45° relative to the low-temperature configuration.7

The ensemble of the above-presented experimental data

provides many metrics against which to assess the performance

of our interaction potential by molecular dynamics simulations.

In this work, we describe such an assessment. In Section n, we

provide a description of our potential model and our choice and

parametrization of the intramolecular and intermolecular interaction terms. In Section ED, we provide a brief description of

computational details. Section IV contains the results of our

MP and NPT-MD simulations and a comparison of predictions

of crystal structural parameters with experimental values over

a wide range of temperatures and pressures. Summary and

conclusions are given in Section V.

II. Potential Energy Functions

We have assumed that the potential energy for a system of

Nnitromethane molecules can be described as the sum of interand intramolecular interaction terms:

N I 1 N \

ytotal ^j yintramolecular V

1

yintennolecular j , »-.

The intermolecular potential is the same as described in the

earlier studies,1-6 and consists of the superposition of a pairwise

sum of Buckingham (6-exp) (repulsion and dispersion) and

coulombic (C) potentials of the form:

V

a^r) =

A

aßexP(-B<xßr) ~

C

a^r

and

v c(

,=ML

(2)

(3)

where r is the interatomic distance between atoms a and ß, qa

and qß are the electrostatic charges on the atoms, and eo is the

dielectric permittivity constant of free space. The parameters
8408 J. Phys. Chem. B, Vol. 104, No. 35, 2000 Sorescu et al.

TABLE 1: Force Field Parameters for Crystalline

Nitromethane'1

atom* charge/e

C, -0.305342

N2 0.820603

03 -0.470445

04 -0.484277

Hs 0.143365

Hs 0.155443

H7 0.140652

Morse potential parameters

bond A(kJ/mol) ß(k~') M)

C-N

N-O

C-H

251.0419414 2.005501

390.3702976 2.459942

426.7713101 1.892486

bending potential parameters

1.499574

1.226747

1.090000

angle fc«(kJ/molrad 2) <?(deg)

C-NO-NN-CH-C--O 294.5211162

-O 657.8485114

-H 224.8140504

-H 149.9402928

torsional potential parameters

117.039960

125.890000

107.560708

111.312289

angle V# (kJ/mol) <5(deg) m

H,-Ci-N2-03 (i = 5,6,7) 0.27000 (0.49)< -90.0 3.0

H,-C,-N2-O40" = 5,6,7) 0.27000 (0A9Y 90.0 3.0

N2-O4-O3-C1 240.37076 180.0 2.0

" The electrostatic charges have been determined by the CHELPG

procedure as implemented in Gaussian 94vat the HF/6-31G** level.

'The atom designation numbers are defined in Figure 1. c The choice

of a torsional barrier of 0.49 kJ/mol is also discussed in text.

Aaß, Baß, and Caß for different types of atomic pairs have been

previously published and have been used in the present study

without change.1

The set of partial charges used in these calculations were

determined through fitting these to the quantum-mechanicallyderived electrostatic interaction potential for an isolated molecule

whose atoms are arranged in the experimental crystallographic

arrangement. These calculations have been done using the

CHELPG procedure as implemented in the Gaussian 94

package.17 We have previously shown45 that for the majority

of the crystals studied, the best agreement between simulations

and experiment occurs when the set of partial charges is

determined using methods that employ electron correlation

effects such as second-order Möller-Plesset (MP2) perturbation

theory.18 However, this was not true for nitromethane; the best

agreement between the measured and the predicted crystallographic lattice parameters was obtained when unsealed partial

charges derived from the Hartree—Fock (HF) wave function

were used. The lack of improvement in the predicted lattice

parameters when electron correlation effects are considered

might be due to omission of other important interactions. These

could include the omission of the polarization effects of

neighboring molecules in the crystal when evaluating the

electrostatic charges for the isolated molecule, or the assumption

of a simple charge-charge interaction model when a multipolar

interaction description would be more accurate. Nevertheless,

our choice for the set of charges determined at the HF level

(see Table 1), together with the rest of the potential parameters,

give a very good representation of the crystallographic parameters.

In a previous study we have developed a complex intramolecular interaction potential to describe the unimolecular

decomposition reactions of gas-phase nitromethane.19 However,

we do not require such a complex function in order to

accomplish the goals described in this work. Thus, we have

assumed a simpler form of the intramolecular interaction

potential:

yiniitramolecular V 4- V + V streich bend torsion (4)

to describe the bond stretching, angle bending, and torsional

motions that occur within an isolated molecule. The bond

stretches are represented by Morse potentials,

V«««* = X Oei{exp[-2A(ri - r?)] - 2 expf-ft^ - /•/>)]}

(5)

1=1

where r, are the bond distances, Dcl are the bond dissociation

energies, /J, are the curvature parameters, and rp are the

equilibrium bond lengths. The six bond stretches that are

described by this potential correspond to the six covalent bonds

in nitromethane.

The bending potentials are represented by harmonic functions

of the form

v

bend = Z;W-ö,0)

2

1=1 2

(6)

where k$ is the force constant and 6° the equilibrium value of

the angle. The nine angles used to describe these interactions

correspond to the HCH, HCN, CNO, and ONO angles in

nitromethane.

Cosine-type torsional potentials have been used to describe

the relative positions of the C—NO2 atoms and the orientations

of the H atoms relative to the C—N—O planes. The potentials

have taken the form

torsion = X^[1+COS(^1.- 6)] (7)

i=l

where V* is half of the intramolecular torsional barrier, <J> is

the torsional angle, and m = 2 or 3. Six of the torsional angles

used to describe these interactions correspond to HCNO angles,

the seventh corresponds to the N—O—O—C torsional angle. The

complete list of potential parameters from eqs 5—7 is given in

Table 1.

The bond dissociation energies in eq 5 have been taken from

our previous study of nitromethane in the gas phase,19 while

the remaining set of geometrical equilibrium values and force

constants in eqs 5—7 were parametrized based on quantum

mechanical information generated using density functional

theory. These calculations have been done for the isolated

molecule using Becke's three-parameter hybrid method,203 in

combination with the Lee, Yang, and Parr correlation functional201'

(B3LYP) and the basis set 6-3 IG*21 using the Gaussian 94

package of programs.17 Our earlier molecular dynamics studies22

on reactions of gas-phase energetic molecules suggested that a

potential energy function fitted to only the structural parameters

and normal-mode frequencies could lead to a model that

produced anomalous and unexpected results, such as nonstatistical behavior in the unimolecular decomposition of the molecule. However, when the function was fitted to the ab initio

Cartesian second derivatives of the energy, the potential model

produced the expected statistical behavior in the unimolecular

decomposition of the molecule. In this study, both the geometrical parameters ofthe optimized nitromethane molecule and
Solid Nitromethane /. Phys: Chem.B;Völ. 104, No. 35, 2000 8409

TABLE 2: Comparison between the Experimental, the Scaled Harmonic Vibrational Frequencies Calculated at B3LYP/

6-31G*Level,aand the Predicted Frequencies by Our Force Fields; the Corresponding Projections (Scalar Product) of the

Predicted Eigenvectors on the Ab Initio Eigenvectors Are Also Indicated*

band

type

ref

24a

ref

24b

ref

24c

ref

24d

ref

24e

B3LYP/

6-31G*0 Mgas Pgas Ml PI M2 P2

ft 51 50 0.94 117 0.94 157 0.94

v% K 416 477 480 475.2 461 459 0.99 461 0.99 461 0.99

V6 A" 599 605 607 602.5 584 596 0.94 601 0.94 605 0.94

V5 A'2 647 658 656.0 655 657.4 636 602 0.92 605 0.92 605 0.92

n *2 921 918 917.1 917 917.9 892 828 0.93 821 0.92 821 0.92

Vll ET 1097 1087 1091.0 1103 1099.0 1076 1052 0.99 1053 0.99 1054 0.92

"u Ef 1153 1100 1146.0 1103 1119.0 1101 1098 0.99 1101 0.99 1103 0.99

"3 A2 1384 1377 1378.8 1379 1380.4 1364 1509 0.95 1508 0.95 1508 0.99

V2 A2 1397 1397.0 1402 1397.4 1388 1372 0.96 1368 0.96 1368 0.95

Vio E 1449 1443 1426 1432 1432 1.00 1434 1.00 1437 0.99

"10 E 1488 1482 1438.0 1426 1444 1434 0.99 1436 0.99 1439 0.99

V7 A'{ 1582 1586 1583.3 1561 1583.8 1615 1647 1.00 1618 1.00 1618 1.00

"1 A2 2965 2972 2964.3 2967 2973.9 2981 2961 1.00 2961 0.98 2961 0.98

V9 ET 2984.0 3045 3044.0 3070 3093 0.98 3093 1.00 3093 0.98

V9 E 3048 3065 3080.4 3065 3080.0 3100 3095 0.98 3095 1.00 3095 1.00

a B3LYP/6-31G* values are scaled by 0.9613, as recommended in ref 23. *The columns denoted with Mgas, Ml, and M2 correspond to the

reference force fields for gas phase, for solid phase (with torsional HCNO barrier of 0.27 kj/mol), and for solid phase with increased torsional

barrier (torsional HCNO barrier of 0.49 kJ/mol), respectively. The corresponding projections of the predicted eigenvectors are indicated in columns

Pgas, PI, and P2, respectively. The frequency units are cm-1.

the ab initio eigenvalues (scaled by a factor of 0.961323) and

the corresponding eigenvectors of the ab initio Hessian matrix

have been used in the fitting procedure.

Special attention was paid in this work to the choice of the

torsional barrier V* corresponding to the H-C—N—O torsional

angles. As indicated in previous experimental studies,8,101416

an exact value for this barrier is not yet well-known and several

values have been proposed. The major difficulties in detenriining

a precise barrier for the internal rotation of the methyl group

were due to the existence of a complex mechanism that

incorporates enhanced hopping due to interactions with phonons,

tunneling effects, and changes of the barrier with temperature.

However, it is important to point out that independent of

experimental technique used, the numerical values of the

activation energy for the methyl rotation are quite small81016

and represent only fractions of 1 kcal/mol. In the present study

we have investigated several values for V* potential parameters

of H-C-N—O torsions and have analyzed their role on the

crystallographic parameters. The vibrational frequencies predicted by our classical potential for an isolated nitromethane

molecule are given in Table 2 along with the corresponding ab

initio values and experimental data. In this table we also indicate

the values of the projections ofthe eigenvectors obtained from

the normal-mode analysis using the proposed potential onto their

quantum mechanical counterpart. A projection whose magnitude

is one would indicate perfect agreement between the model and

the B3LYP/6-31G* eigenvectors. A projection that has a

magnitude near zero indicates that the atomic motion of the

vibration that is predicted by the model is extremely different

from that predicted by the quantum mechanical calculations.

The projections ofthe normal-mode eigenvectors generated by

different potential models (see below) onto the quantum

mechanical eigenvectors are also indicated in Table 2.

In an initial attempt we described the internal methyl group

rotational motion using a 6-fold torsional potential with a barrier

of 6 cal/mol, in agreement with the experimental gas-phase

barrier.15 The vibrational frequencies corresponding to this

model are indicated in Table 2 where they are denoted as Mgas.

As can be seen there is an overall very good agreement between

the entire set of predicted vibrational frequencies and the

corresponding ab initio results. In particular, the frequency of

the torsional mode vris practically identical to the corresponding

quantum mechanical value. The corresponding projection of the

eigenvector predicted by the model onto the quantum mechanical eigenvector is 0.94, indicating also a very good prediction

of the vibrational motion associated to this mode. However,

when this potential was tested in NPT-MD simulations in the

regime oflow temperatures and pressures it was found that while

the unit cell size and shape of nitromethane were reasonably

predicted, the methyl groups of the molecules in the unit cell

were rotated by approximately 30° relative to the experimental

orientation. This fact indicated the need to increase the barrier

height of the torsional H-C-N-O potentials in order to

reproduce accurately the crystallographic structure.

The second set of potentials we have considered is similar

to the Mgas model but we have changed the HCNO torsional

potential to a 3-fold form and increased the barrier parameter

Vd,(H-C-N-0) to 0.27 kJ/mol, which corresponds to an

overall barrier for methyl rotation similar to that given in ref

14. The complete list of intramolecular parameters for this

potential set is that indicated in Table 1. The modifications of

the corresponding vibrational frequencies for this model are

shown in columns denoted as Ml and PI in Table 2. It can be

seen that the most important variations of the vibrational

frequencies take place for the mode vz, which corresponds to

the internal rotation of the methyl group. Overall, excepting

modes vT and v3 (which describes the inversion of the methyl

group), there is good agreement between the predicted, the ab

initio calculated, and the experimental24 sets of vibrational

frequencies with rms deviations of 6.9 and 10.8 cm-1, respectively. In addition, the results indicate thatthe model reproduces

the quantum mechanical eigenvectors very well, with the

smallest projection values being 0.92 for modes v4 and vs. As

will be described more fully hereafter, MP calculations and

NPT-MD simulations using this potential model reasonably

reproduce the experimental crystal structure and molecular

conformations for the entire range oftemperatures and pressures

investigated. Consequently, all the results presented in the next

sections correspond to this choice of potential parameters

denoted as model Ml.

A final choice oftorsional potential parameters we have tested

corresponds to an even larger value for V<j>(H—C—N—O) of

0.49 kJ/mol. Our results indicate that this potential change has

only a very small influence on the predicted crystallographic
8410 J. Phys. Chem. B, Vol. 104, No. 35, 2000 Sorescu et al.

TABLE 3: Lattice Parameters and Energies Obtained in Crystal Packing without Symmetry Constraints (the percentage

change in lattice and molecular parameters after energy minimization is determined as a function of the experimental geometry

and is given in parentheses)

lattice energy" lattice parameters*

Q

c NB ES a b c a ß 7

EXP*

LMIN

20.5

-52.3'

-25.0401 -35.0560

5.1832

5.2369

(1.04)

6.2357

6.2653

(0.48)

8.5181

8.6214

(1.21)

90.000

90.001

(0.00)

90.000

89.989

(-0.01)

90.000

89.997

(0.00)

A/ Ay Az A<J> (deg)* A© (deg) A«/ (deg)

-0.037 -0.105 -0.014 3.44 3.02 6.41

" Nonbonded (NB) and electrostatic (ES) lattice energies per molecule in kJ/mol. * Lattice dimensions a. b, and c in angstroms, and angles a, ß,

and y in degrees.cCutoff parameter as described in text. "* Experimental values from ref 7.' Estimated using the theoretical heat of sublimation of

—47.3 kJ/mol as provided by Politzer.34 sChange in fractional coordinates of molecular centroids for nitromethane molecule. s Change in Euler

angles of molecular centroids for nitromethane.

parameters relative to the Ml model and that this influence is

particularly limited to the region of low temperatures. Consequently we will make only limited references to this potential

choice in the present work. For completeness, we indicate in

Table 2 columns M2 and P2, the predicted vibrational frequencies and the corresponding projections for this potential model.

As can be seen only the value of the torsional frequency vris

increased by about 40 cm-1 relative to the Ml value while the

projection of its eigenvector on the corresponding ab initio

vector remains 0.94. The large value of this projection, close

to unity, indicates that even for this high barrier value the model

is representing correctly the corresponding vibrational atomic

motions.

m. Computational Details

Molecular Packing Calculations. A first set of calculations

used to test the semiempirical intermolecular potential energy

functions proposed here is based on the use ofmolecular packing

calculations,25,26 in which the lattice energy of the crystal is

rninimized with respect to the structural degrees of freedom of

the crystal. These calculations have been done using the

algorithm proposed by Gibson and Scheraga27 for efficient

miriimization ofthe energy of a fully variable lattice composed

of rigid molecules and implemented in the program LMIN.28

The nonbonded interactions were cut off with a cubic spline

function from Po to Qo, to ensure the continuity of the function

and its first derivative. Here a is the value ofrin eq 2 at which

Vafi(r) = 0 and dVaß(r)/dr < 0. The parameters P and Q, which

specify the start and the end of the cubic feather (see refs 1 and

27 for details), were set to 20.5 and 20.0, respectively. The

coulombic potential terms of the form given in eq 3 are summed

over the lattice using the Ewald technique as previously

described.1Finally, the effect of pressure on the crystallographic

parameters has been simulated by adding a potential term of

the form P(V - Vb),27 where Vo is the volume of a suitably

chosen unit cell at zero pressure.

Constant-Pressure and -Temperature Molecular Dynamics Calculations. The dynamics of nitromethane as a function

oftemperature and pressure have been investigated using NPTMD simulations based on the total potential described in eqs

1—7. For comparison we have also performed calculations using

the rigid-body approximation1-6 ofthe molecules in the system.

In all simulations we have used the Nose-Hoover thermostat—barostat algorithm29 as implemented in the program

DL_POLY_2.0,30 to simulate the crystals at various temperatures and pressures. In this case, the equations of motion for

molecules and the simulation cell are integrated using the Verlet

leapfrog scheme.31 In the case of rigid-molecules simulations

the molecular rotational motion is handled using Fincham's

implicit quaternion algorithm.32

The MD simulation cells consist of boxes containing 5x4

x 3 crystallographic unit cells. This choice of the simulation

box ensures the use of a cutoff distance for the intermolecular

potentials of about 10 A. The initial configuration corresponding

to the lowest temperature was chosen to be identical to that for

the low-temperature experimental structure. The system was then

equilibrated at that temperature and atmospheric pressure. In

all production runs done using the Nos6-Hoover implementation

for the NPT ensemble, the system was integrated for 34000

time steps (1 time step = .75 x 10-15 s), of which 4000 steps

were equilibration. In the equilibration period, the velocities

were scaled after every 5 steps so that the internal temperature

of the crystal mimicked the imposed external temperature. Then,

properties were calculated and accumulated for averaging over

the next 30000 integration steps in the simulation. In subsequent

runs, performed at successively higher temperatures or pressures,

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. In these calculations, the coulombic longrange interactions were handled using Ewald's method.31

The main quantities obtained from these simulations were

the average lattice dimensions and the corresponding volume

of the unit cell. Additional information about the structure of

the crystal has been obtained by calculating the radial distribution functions (RDF) between different atomic sites. Such

quantities have been calculated from recordings done at every

10th step during the trajectory integrations.

IV. Results and Discussions

A. Molecular Packing Calculations. The results of MP

calculations without symmetry constraints are presented in Table

3. As can be seen the predicted structural lattice parameters for

nitromethane differ by less than 1.21% from the experimental

values reported by Trevino et al.7 Also, there are small

translations but slightly larger rotations (with a maximum

deviation of 6.4°) of the molecule in the asymmetric unit cell.

Using the HF set of charges, the total lattice energy is predicted

to be —60.1 kJ/mol. This value compares acceptably well with

the lattice energy of -52.3 kJ/mol that has been estimated using
Solid Nitromethane J. Phys. Cherrt. B,Vol. 104, No. 35, 2000 8411

TABLE 4: Lattice Parameters Obtained in NPT-MD Calculations for Nitromethane as a Function of Temperature; the

Calculated Thermal Expansion Coefficients (jf) at 250 K Are Also Indicated"

lattice dimensions

T(K) a(Ä) Kk) c(Ä) a(deg) £(deg) y(deg) volume(Ä3)

4.2* 5.1832 6.2357 8.5181 275.3125

78.0" 5.1983 6.2457 8.5640 278.0476

228.0* 5.2440 6.3200 8.7300 289.3303

4.2 5.2116(0.54) 6.3416 (1.69) 8.6464 (1.50) 89.998 89.999 90.000 285.7666 (3.79)

25.0 5.2134 6.3608 8.6935 90.012 90.018 89.996 288.3000

50.0 5.2184 6.3794 8.7266 90.001 90.034 90.000 290.5000

78.0 5.2195 (0.40) 6.4032 (2.52) 8.7543 (2.22) 89.991 89.996 89.991 292.5833 (5.22)

100.0 5.2222 6.4173 8.7734 90.001 89.962 90.007 294.0166

125.0 5.2267 6.4329 8.7979 89.987 89.994 90.002 295.8166

150.0 5.2298 6.4605 8.8204 89.974 89.991 90.010 298.0166

175.0 5.2373 6.4889 8.8606 89.995 90.005 89.987 301.1166

200.0 5.2452 6.5167 8.8808 90.016 90.012 90.020 305.3166

228.0 5.2534(0.18) 6.5481 (3.60) 8.9195 (2.17) 90.015 90.020 - 90.043 303.5667 (4.92)

250.0 5.2620 6.5691 8.9561 90.001 89.993 89.946 306.8167

t 68.9 x 10"6 181.7 x 10"6 131.9 x lO"6 380.5 x 10-6

" The values in parentheses represent the percentage differences relative to the available experimental results. * Experimental data from ref 7.

c The units for the linear and volume expansion coefficients are K_I

.

the relation U\m *w —Affsobi -2Ä7".33 In this analysis we have

used as the heat of sublimation for nitromethane the value of

47.2 kJ/mol estimated theoretically by Politzer.34

B. NPT-MD Calculations. NPT-MD simulations of nitromethane were used to predict the crystal structure of nitromethane over a large range of temperatures and pressures.

We will first examine thermal effects on the crystallographic

cell and molecular parameters at 1 arm, and compare them

against experimental values.

Bl. Temperature Effects. Average unit cell edge lengths and

volumes determined from NPT-MD simulations for the flexible

model for temperatures ranging from 4.2 to 250 K are given in

Table 4. These data are represented in Figure 2, along with the

corresponding experimental values and the NPT-MD results

obtained using the rigid-body approximation. For the flexible

model, the NPT-MD lattice dimensions obtained at T = 4.2 K

are in very close agreement with those determined in the MP

calculations and with the experimental values.7 At this temperature, the percentage differences between the predicted and the

experimental lattice dimensions are 0.54%, 1.69%, and 1.50%

for the a, b, and c axes, respectively. The deviation of the

predicted unit cell volume from the experimental value at this

temperature is 3.79%. This good agreement is decreased only

slightly with an increase in temperature. At T = 78 K, the cell

edge lengths and the unit cell volume agree with the experimental values to within 2.52% and 5.22%, respectively. At the

highest temperature where experimental data are available (228

K), the deviations from experiment for the cell edges are less

than 3.60% while the difference is 4.92% for the unit cell

volume. Also, the unit cell angles remain close to 90.0° for the

entire temperature range analyzed, as expected for the orthorhombic symmetry of the crystal. For comparison, calculations

using the M2 potential model indicate slightly larger deviations

for the predicted lattice parameters from experimental values

of 1.31%, 1.03%, and 1.72% at T= 4.2 K and 0.60%, 3.61%,

and 2.13% at T = 228 K for the a, b, and c axes, respectively.

The linear and volume thermal expansion coefficients extracted from the data are also given in Table 4. The expansion

ofthe lattice is highly anisotropic with the largest length changes

along the b and c axes. These effects can be partially attributed

to the internal rotation of the methyl group. To illustrate our

argument, consider the limiting case in which a nitromethane

molecule is arranged in a unit cell such that the C-N bond is

parallel to the a axis. Changes in the positions of the hydrogen

9.5

^9.0

S 8.5 H

8.0 /

6.5-

6.0 -

5.5

5.0

4.5

QOO QOOOOO 00

• a(EXP)

o a (MD, R)

0 a (MD. F)

b(EXP)

A b (MD, R)

A b <MD, F)

c(EXP)

V c(MD. R)

V C(MD. F)

0 50 100 150 200

Temperature (K)

250 300

330

320-

«£.310-

o

E 300

290

280

270

260

O

>

250

A

D

A

I*

• EXP

A MD, R

D MD, REL

50 100 150 200 250 300

Temperature (K)

Figure 2. Variation ofthe lattice dimensions (a) and unit cell volume

(b) with temperature from NPT-MD simulations using the rigid-body

approximation (MD, R) and the flexible model (MD, F). The available

experimental data from ref 7 (EXP) are also indicated.

atoms due to the rotational motion of the methyl group about

the C—N bond would occur in the b—c plane only. Thus, the

increasingly large librational motion occurring with increases

in temperature would affect only the cell dimensions associated

with the b and c axes. In the actual crystal, the arrangement of

the molecules in the unit cell is such that the C—N bonds of all

four are nearly perpendicular to the c axis (as can be seen in

Figure 1). The angles that each C—N bond makes with the axes

of the crystal cell are 40, 51, and 86° for the a, b, and c axes,

respectively. This alignment of the molecules relative to these

axes indicates that rotations of the methyl group about the C—N

bond would more strongly affect the b and c axes. Indeed, by

comparing the data in Figure 2a, it is clear that the major

differences between NPT-MD data from the rigid and flexible
8412 J. Phys. Chem. B, Vol. 104, No. 35, 2000 Sorescu et al.

TABLE 5: Comparison of the Average Orientational Parameters (fractional coordinates sx, sy, sz and Euler angles 0, tf», V)

for the Four Molecules in the Unit Cell and of the Internal Geometrical Parameters of Nitromethane Molecule as Obtained

from Trajectory Calculations at Different Temperatures; Corresponding Experimental Values (where available) Are Also

Indicated"

temperature(K)

4.2 4.2 4.2 78 78 78 150 150

parameter expt*' NPT-MD ideal expt*' NPT-MD ideal NPT-MD ideal

sxl 0.3493 0.3513 0.3475 0.3483 0.3484

syl 0.4196 0.4284 0.4188 0.4288 0.4290

szl 0.3750 0.3682 0.3803 0.3722 0.3740

sx2 0.1507 0.1487 0.1487 0.1525 0.1509 0.1517 0.1512 0.1515

sy2 -0.41% -0.4284 -0.4284 -0.4188 -0.4291 -0.4288 -0.4290 -0.4290

sz2 -0.1250 -0.1318 -0.1318 -0.1197 -0.1279 -0.1278 -0.1262 -0.1260

sx3 -0.3493 -0.3513 -0.3513 -0.3475 -0.3494 -0.3483 -0.3488 -0.3485

sy3 -0.0804 -0.0717 -0.0716 -0.0812 -0.0711 -0.0712 -0.0705 -0.0710

sz3 0.1250 0.1317 0.1318 0.1197 0.1277 0.1278 0.1260 0.1260

sx4 -0.1507 -0.1488 -0.1487 -0.1525 -0.1513 -0.1517 -0.1516 -0.1515

sy4 0.0804 0.0716 0.0716 0.0812 0.0708 0.0712 0.0703 0.0710

sz4 -0.3750 -0.3683 -0.3682 -0.3803 -0.3719 -0.3722 -0.3739 -0.3740

01 90.8 94.4 96.6 92.9 92.2

4)1 53.5 51.9 51.6 53.2 53.7

VI 238.9 245.3 240.1 245.4 245.4

02 -90.8 -94.4 -94.4 -96.6 -92.9 -92.9 -92.4 -92.2

<t>2 53.5 52.0 51.9 51.6 53.1 53.2 53.7 53.7

W2 238.9 245.3 245.3 240.1 245.5 245.4 245.4 245.4

03 89.2 85.5 85.6 83.4 87.2 87.1 87.7 87.8

4>3 -53.5 -52.0 -51.9 -51.6 -53.1 -53.2 -53.8 -53.7

V3 58.9 65.3 65.3 60.1 65.4 65.4 65.4 65.4

04 -89.2 -85.6 -85.6 -83.4 -87.1 -87.1 -87.7 -87.8

4>4 -53.5 -52.0 -51.9 -51.6 -53.2 -53.2 -53.7 -53.7

tp4 58.9 65.3 65.3 60.1 65.4 65.4 65.4 65.4

temperature(K)

4.2 4.2 78 78 150

parameter expt* NPT-MD expt.' NPT-MD NPT-MD

C-N 1.481 1.497 1.488 1.499 1.499

C-H 1.098 1.088 1.077 1.089 1.090

N-03 1.209 1.230 1.198 1.230 1.231

N-O4 1.223 1.229 1.226 1.230 1.231

HCH 111.6 110.7 111.2 110.6 110.5

HCN 107.3 108.2 107.7 108.2 108.2

CNO3 118.9 118.1 119.1 117.8 117.8

CNO4 117.8 117.7 116.8 117.8 117.9

O3NO4 123.3 124.3 124.1 124.2 124.2

" In the case of molecules 2—4 the ideal orientational parameters obtained based on P2]2i2i space group symmetry operators and the orientation

parameters of molecule 1, positioned in the asymmetric unit, are given. * Experimental values from ref 7. c Unit cell fractional values have been

shifted to fall within the range from —1/2 to 1/2.

potential models are in the expansions along the b and c axes

while the variations for the case of the a axis are practically

absent. The additional expansion due to the methyl group

rotation can be also seen in Figure 2b where the volume values

obtained using the flexible model are consistently larger than

the corresponding data obtained for the rigid model.

In Table 5 we compare first the results of the average

fractional coordinates and orientational Euler parameters of the

four molecules in the unit cell with the corresponding experimental data (where available). These values are averaged over

time and all unit cells in the simulation box. Additionally, in

the same table we provide the "ideal" orientational Euler

parameters for three of the four molecules ofthe unit cell relative

to the orientation of the reference molecule in the asymmetric

unit cell. These parameters have been calculated assuming the

symmetry operations of the P1{1\1\ space group. In this way

we can analyze the degree of deviation of the predicted crystal

structure from the P2\2\2\ symmetry.

It is evident that increasing the temperature from 4.2 to 150

K does not produce any significant displacement of the

molecular centers of mass and that the degree of rotational

disorder is small. The largest deviations of the system from the

orientational parameters corresponding to a perfect crystal with

P2i2]2i symmetry are 0.0011 and 0.2° for the fractional

coordinates and Euler angles, respectively. A similar conclusion

is obtained for temperatures between 150 and 250 K (not

shown). At T = 4.2 K, the largest deviation between the

experimental and predicted molecular orientations is about 6.4°

for the Euler angle W (see Table 5), in agreement with previous

MP findings. At T = 78 K this maximum deviation decreases

to about 5.4°.

Additional support for the small degree of translation of the

molecules inside the unit cell with the temperature increase can

be obtained from the C—C RDFs given in Figure 3a. As can

be seen, the RDFs at these temperatures correspond to wellordered 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.

The good agreement of the structural molecular parameters

predicted by the model with experiment can be also observed

by analyzing the internal coordinates (bond lengths and angles)

indicated in Table 5. At 4.2 K the maximum deviations for bond

lengths is about 1.7% and about 1.0% for bond angles. This
Solid Nitromethane J. Phys. Chem. B, Vol. 104, No. 35, 2000 8413

a:

8

6

4

2

OH

C...C A A ^ a,

T=250K /|\ '"^' ^

T=J50_K_

T=78K

T=4.2K

/ LL

8-

6

4-

C...C

P=7.0 GPa

P=5.4 GPa

Jk*H

?

P=3.5 GPa

••*/\A""

0

P=0.3 GPa _h^y^

10

r(A)

4 6

r(Ä)

10

Q

T= 4.2 K

"- 4 ^

2 -

O...H .

P=7.0GPa/ /^

d)

\

P=5.4 GPa/

P=3.5 GPa/

P=0.3 GPa,

\.

4 6

r(Ä)

8 10 0 2 4 6

r(A)

8 10

Figure 3. Radial distribution functions for O-C and O—H as functions oftemperature (a-b) at atmospheric pressure and as function of pressure

(c-d) at r = 293 K.

level of agreement is particularly good taking into account that

in our potential parametrization, we have considered the internal

geometrical parameters ofthe nitromethane molecule determined

from gas-phase ab initio calculated values rather than parameters

corresponding to nitromethane in the crystalline phase. These

results suggest also a small influence of the lattice field on the

internal parameters of the nitromethane molecule. Changes in

predicted molecular structural parameters with temperature are

negligible, again in agreement with experiment. The most

significant change appears due to the internal rotation of the

methyl group as a function of temperature. This is evident by

inspecting the behavior of the radial distribution functions

(RDFs) for the O—H intermolecular bonds given in Figure 3b.

There is a rapid destruction of the correlation at long distances

with the increase of temperature and consequently the individuality of the H atoms in the C-H—O bonds is lost. These

effects can be understood as being due to the increase in

rotational disorder of the methyl group with the increase of

temperature.

The increase in the librational motion of the methyl group

with temperature is illustrated in Figure 4. This figure shows

the cumulative distributions of the H,—Q-N2-O3 (/ = 5-7)

dihedral angles at temperatures ranging from 4.2 to 150 K. These

distributions have been determined from the molecular configurations recorded at every tenth step during trajectories of

16000 steps (12 ps). In this figure, the distributions are peaked

at the angles corresponding to the equilibrium orientation of

the methyl group in the low-temperature structure, namely,

-78°, 42°, and 161° respectively. These angles are within 10°

from those corresponding to the experimental structure determined by the neutron diffraction technique.7 As can be seen in

Figure 4 the increase oftemperature does not change the position

of these peaks. However, the distributions broaden with temperature, indicating an increase in librational motion with

temperature. For the case of the potential model with the

increased HCNO torsional barrier (M2) we have found that the

agreement between the predicted (-85°, 34°, 154°) and

experimental (-88°, 31°, 151°) torsional angles is even better,

with deviations within 3°. Similarly, in this case the increase

of temperature did not modify the position of the peaks in the

cumulative distributions.

The magnitudes of the methyl group libration have been

previously determined experimentally for deuterated nitromethane in the temperature range 4.2-125 K using neutron

diffraction powder experiments.7 As can be seen in Figure 5

where the experimental values are compared with our results

the agreement is quite good for the entire temperature range

25-125 K. For example, the experimental results indicate that

at 25 and 125 K the rms amplitudes of libration of the CD3

group are ~13° and 25°, respectively. Our NPT-MD calculations

of deuterated nitromethane predict rms amplitudes of libration

of 12° at 25 K and 23° at 125 K (see Figure 5). Comparison of

the amplitude of libration at 4.2 K with experiment is not

appropriate, since quantum mechanical effects would be significant at this temperature. The experiments7

suggest that the

amplitude of libration at this temperature is as large as that at
8414 J. Phys. Chem. B, Vol. 104, No. 35, 2000 Sorescu et al.

c

O 0.04

'•-*

3

JD

is 0.03

W

Q

_g>

c

<

to

c

o

0.02

0.01

0.00

a) I

T=4.2 K

T=78 K

V ^T=150K U

„ ,yf V ,=4TCT=250Kyfl

-200 -150 -100 -50 0 50 100 150 200

Torsional Angle (deg.)

C 0.007

O

Z5 0.006

.Oi

000S

Q1-1 0.004

j? 0.003

<

— 0.002

CO

c

O 0.001

"co

O 0.000

H- -200 -150 -100 -50 0 50 100 150 200

Torsional Angle (deg.)

Figure 4. Distribution of the H,—Q-N2-O3 (i = 5,6.7) dihedral

angles for all nitromethane molecules in the simulation box as a function

of temperature (a) and pressure (b).

30

25

20

P=0 3GPa Gpa

b) jj P=7.0 GPa a 1 P=3.5 GPa

- >

i\ -^ A A /

s*>

—i 1—1

/¥- k

r\ './

mw / Af^

1 *

v

*V il

\ **""-^W/ \

W7V/ \ '

v

"*'V/ w vy v >-•'

a

«

3 10

en

E

5

• Expt

V NPT-MD

•

• * V

• 9

V

V

• V

V

20 40 60 80

Temperature (K)

100 120 140

Figure 5. Root-mean-square amplitude oflibration ofthe methyl group

in CD3NO2 as a function of temperature.

25 K; our classical simulations predict a value that is smaller

by ~7°. Previous quasielastic and inelastic neutron scattering

studies8performed on protonated nitromethane have provided

additional information about the temperature dependence of the

rotational reorientation of the methyl groups. By assuming a

3-folded reorientation model and an Arrhenius dependence

(~e£a/<:7) of the mean residence times of methyl groups on the

temperature T, an activation barrier of 234 cal/mol to the methyl

group rotation has been determined. However, more recent

experimental studies based on infrared and Raman measurements have indicated higher activation energies with values of

401 cal/mol16 and 576 cal/mol,10 respectively. We monitored

0.002 0.006 0.010 0.014 0.018 0.022

l/T(Kr')

Figure 6. The mean residence times r of the methyl group versus

inverse temperature. The activation energy of 387 cal/mol is obtained

from the slope ofthe least-squares linear fit to the data in the region of

higher temperatures (125-250 K) indicated by a solid line.

the frequency of 120° methyl group rotations for the duration

of 12 ps in trajectories integrated at temperatures ranging from

4.2 to 250 K. From these trajectories we have observed that

methyl groups undergo large vibrational amplitudes, indicative

of small barrier heights for torsional motions in agreement with

experimental findings.

To obtain a more quantitative description of the internal

rotational motion we have considered the evaluation of temperature effects on the rates of reorientation of the methyl group.

A "rotational jump" was defined in the following manner. At

the beginning of each trajectory, the H—Q—Ni—O3 angles (1

= 5—7) of each molecule were calculated. Each angle was

assigned to one of the three minima in the 3-fold rotational

potential by virtue of its proximity to the minima, and the time

denoted as to. The minima corresponded to the values of the

H,—Q-N2-O3 angles (/' = 5—7) in the low-temperature

structure, i.e., —78°, 42°, and 161°, respectively. The assignment

of these angles to the minima are made at each subsequent step

in the trajectory and compared to the original assignment at f0-

When the assignments for all three ofthe H,—C\—N2-O3 angles

differ from the original assignments at to, it is assumed that a

rotational jump has occurred. At this point, the time interval

between the time of the new assignments and fo was recorded.

The value of fo was reset to the time of the new assignments,

and the process repeated for the duration of the trajectory.

Methyl group reorientation rates at each temperature were

determined by extracting first-order rate coefficients from linear

fits of plots of ln(F) versus time, where P is the fraction of

nitromethane molecules that have not jumped at time t. The

lifetimes indicate first-order behavior. An Arrhenius plot of the

residence times, which are the inverse of the reorientation rates,

is shown in Figure 6. Fitting a line to all the points in Figure 6

predicts £a = 241 cal/mol. However, the nonlinearity of the

Arrhenius plot is a result of the statistical errors in the lowtemperature rates, and if we base our prediction on the 7 points

for the highest temperaures (100-250 K) we predict £a = 387

cal/mol. These values are in accord with the reported experimental activation energies (234—576 cal/mol).81016

B2. Pressure Effects. The calculated lattice dimensions at

different pressures are given in Table 6 and a visual comparison

with the experimental data is presented in Figure 7. For the

entire pressure range investigated (0.3—7 GPa) the differences

between the calculated NPT-MD data using the flexible potential
Solid Nitromethane J. Phys; Ghenu B,Vol. 104, No. 35, 2000 8415

TABLE 6: Lattice Parameters Obtained in Crystal Packing and NPT-MD Calculations for Nitromethane as a Function of

Pressure

NPT-MD molecular packing

press (GPa)

0.30

0.33

0.50

1.00

2.00

3.00

3.50

4.00

5.00

5.45

6.00

6.50

7.00

10

lattice lengths (Ä) unit cell

volume (Ä3)

lattice lengths (A) unit cell

a b c a b c volume (A3

)

5.2070 6.4829 8.8597 299.0666 5.2068 6.1989 8.5628 276.3758

5.2006 6.4874 8.8416 298.2833 5.2040 6.1929 8.5576 275.7927

5.1713 6.4236 8.7800 291.6500 5.1861 6.1561 8.5172 271.7649

5.0997 6.3001 8.6595 278.2166 5.1397 6.0786 8.4295 262.9498

4.9827 6.1706 8.4936 261.1500 5.0656 5.9716 8.3055 250.5347

4.8963 6.0797 8.3826 249.5166 5.0051 5.8988 8.2150 241.6301

4.8603 6.0448 8.3392 245.0000 4.9792 5.8684 8.1772 237.9397

4.8245 6.0139 8.2955 240.6833 4.9536 5.8447 8.1421 234.6472

4.7652 5.9600 8.2265 233.6500 4.9054 5.8085 8.0780 228.8875

4.7426 5.9412 8.2001 231.0500 4.8843 5.7969 8.0517 226.5996

4.7116 5.9145 8.1666 227.5833 4.8527 5.7982 8.0161 223.9448

4.6898 5.8925 8.1404 224.9666 4.7698 5:9135 7.9484 221.2719

4.6666 5.8733 8.1162 222.4500 4.7171 5.9801 7.8992 218.8226

0)

c

o

"5

c

o

E

a

a

u

£

CO

a)

y ^'^^w*, /

/

ttkfrr.

/

^■■---A-fr.frr:.--a-frr-i -4-4

"CD--) •0 S-3_$.$_-;

• • -<BC->

o a(MP)

T a (MD.R)

-V- a (MD, F)

• ■■ b(EXP)

o b(MP)

A b(MD. R)

-A- b(MD, F)

• • -(EXP)

o c(MP)

• c (MD. R)

-o- c (MD. F)

0 2

Ä 320

•"" 3CKH o

I 280-

Ö

> 260-

U 240

= 220-

200

fe*f=fa^

••• D(P

— - EXP2

o MP

T MD.R

-*?- MO.F

••• EXP

O MP

T MD.R

-V MD.F

Figure 7. Pressure variation of the lattice dimensions (a), unit cell

volume (b), and volume compression V/Vb on the external pressure.

The calculated data from molecular packing (MP), NPT-MD simulations

using the rigid-body approximation (MD,R) and the flexible model

(MDJO are represented. The corresponding experimental data from refs

9 (EXP) and 12 (EXP2), respectively, are also indicated.

and the experimental crystallographic parameters are small. By

considering as reference the experimental data determined by

Cromer et al.9 we find that in the low-pressure region (0.3 GPa)

the percentage errors for lattice dimensions a, b, and c are 0.3%,

3.4%, and 2.9%, respectively. A similar good correspondence

is found in the high-pressure region (6.0 GPa) where the

percentage errors are 0.4%, 3.5%, and 2.4% for the a, b, and c

lattice dimensions, respectively. Moreover, for the entire pressure range investigated the NPT-MD predicted volume compression closely follows the corresponding variation observed

experimentally (see Figure 7b). In Figure 7b we have also

represented the experimental values reported by Yarger and

Olinger.12 As can be seen in this case the agreement with our

calculated unit cell volumes is even better than that found with

data provided by Cromer et al.,9 with decreasing deviations from

about 3.6% at 1 GPa to 1.4% at 7 GPa.

The MP results obtained using the rigid body approximation

generally follow the predictions obtained using the flexible

potentials, particularly in the low-pressure region. Also, the

agreement between the NPT-MD data obtained using the rigid

molecular assumption is surprisingly good with NPT-MD data

using the flexible model. This fact indicates that for the pressure

region investigated, the compression of the lattice is almost

entirely due to a reduction of intermolecular distances. These

effects can be also seen from the inspection of the normalized

dependence ofthe unit cell volume on pressure given in Figure

7c. In this plot the experimental curve9and those predicted by

the flexible and rigid models are practically superimposed for

the entire range of pressures investigated. The values obtained

in the MP calculations start to deviate more significantly from

the experimental results, particularly for pressures above 5 GPa.

These facts indicate thatthe lattice compressibility becomes less

well represented in these calculations as the pressure is

increased. However, the present set of results support our

previous findings6that in the region oflow to moderate pressures

(~5 GPa) the MP calculations can be used as an alternative

tool to describe the changes of the unit cell geometrical

parameters but at a fraction of the computational cost involved

in MD simulations.

Further comparison ofthe predicted and experimental lattice

parameters can be obtained by extrapolating to zero pressure

and 293 K the results of different studies, which in principle

should coincide. As in previous experimental studies,7'

9,12 the

extrapolation is taken because at room temperature nitromethane

is in the liquid phase. The coefficients of cubic fits in pressure

of the predicted lattice parameters and unit cell volume to the

data predicted by our flexible potential are given in Table 7. In

Table 8 we compare the zero-pressure extrapolated data obtained

in this work to those of previous experimental investigations.

Particularly, we have considered both the low-temperature data

obtained by Trevino et aL7and extrapolated to room temperature

as well as the high-pressure values determined previously9'

12

and extrapolated to zero pressure. As can be seen in Table 8,

our crystal lattice dimensions are within 2.2%, 2.9%, and 3.4%

from the data given in refs 7, 9, and, 12, respectively. This level

of agreement is satisfactory taking into account the relative large

deviations between different sets of experimental values.
8416 /. Phys. Chem. B, Vol. 104, No. 35, 2000 Sorescu et al.

TABLE 7: Coefficients of the Cubic Fit in Pressure (GPa) of the Lattice Constants and Unit Cell Volume

CXo <Xi a2 a3

«(A)

b(A)

c(Ä)

V(A3)

5.246 A

6.537 Ä

8.911 Ä

305.17 A

-1.601 x KT1 A GPa-1

-2.475 x lO"1 A GPa"1

-2.711 x lO-'ÄGPa"1

-2.885 x 10+1 A

3 GPa"1

TABLE 8: Comparison of the Lattice Constants for Solid

Nitromethane Extrapolated to 293 K and 0 GPa and of the

Murnagham Equation Coefficients (eq 8) from Present

Calculations and Various Studies

parameter

ref7

(low

temperature)

ref9

(high

pressure)

reflO

(high

pressure)

present

work

oo(A) 5.270 5.227 5.197 5.238

*>o(A) 6.375 6.330 6.292 6.514

co(A) 8.832 8.758 8.747 8.889

v„(A3) 296.72 289.77 286.02 302.80

VWÄ3) 292.7 291.87 310.35

ßo(GPa) 7.0 9.14 6.78

&0 5.7 6.04 5.88

Another set of parameters we have used to assess the accuracy

of our potential is represented by the coefficients of the

Mumaghan equation:935

*»«-] (8)

which has been used to fit the dependence of the unit cell

volume on pressure. In eq 8, V is the volume at pressure P, V0

is the volume at P = 0 (and allowed to be varied in the fitting

procedure), BQ is the bulk modulus at zero pressure, and B^ =

dBo/dP. In Table 8 we present the best-fit parameters obtained

from our NPT-MD calculated values together with the corresponding experimental data. Coir results clearly support the

experimental findings obtained by Cromer et al.9 with relative

percentage deviations of 3.2% and 3.1% for the bulk modulus

and its pressure derivative values, respectively.

From the compression data for our lattice parameters we have

also determined the variation of cell edge moduli, -X 9P/9X,

as a function of pressure (see Figure 8). Our results indicate

that the compressibilities of the b and c axes are close to each

other, the c axis being the least compressible at high pressures.

These findings are similar to those obtained by Cromer et al.9

and represented for comparison in Figure 8. However, the

differences in compressibilities for the b and c axes with

pressure are more pronounced in the experimental case than

our data show. Additionally, we find that the a axis is the least

compressible at low pressures but becomes the most compressible for pressures above 2 GPa. This trend is also present in the

Cromer et al. data9and has been also noted previously by Yarger

and Olinger.12 These authors12 have suggested that the behavior

of the compressive response of the unit cell, particularly with

regard to the a axis, could be explained by the formation of

"chains" of molecules, linked by bonds between the nitro and

methyl groups of the neighboring molecules as they are

compressed. Unfortunately, there is no additional theoretical or

experimental support to identify such a mechanism and further

studies, particularly based on ab initio calculations, are necessary

to explore the formation of such chainlike structures. The only

available data are related to a density functional theoretical study

performed at the B3LYP/6-31G** level to describe the interaction between two nitromethane molecules separated by distances

ranging from 1.6 to 3.5 A.36 The molecules were initially

arranged with the methyl group on one molecule facing the nitro

1.704 x 10~2 A GPa"2

3.836 x 10"2ÄGPa-2

3.788 x 10~2 A GPa"2

4.123 A

3 GPa"2

140

120 -

CO

a. 100

2.

»

| 80

o

o

o 60

O)

•aUJ

= 40

ID

o

20

-0.858 x 10~3 A GPa"3

-2.379 x lO-3 A GPa"3

-2.210 x 10"3AGPa-3

-2.428 x 10"1 A

3 GPa"3

c

/A^

^^^ C (exp)

b(exp)

-

/s "

,-—"" a

-

^^^afexp)

- ß/ -/J

1 1 1 i 1 1

12 3 4 5 6 7 8

Pressure (GPa)

Figure 8. The variation of the cell edge moduli as a function of

pressure from the present calculation and from experimental data given

in ref 9.

group on the second (an arrangement that would correspond to

compression along the a axis). As the molecules were brought

closer together, the potential energy first showed a steep

increase, and the methyl group was "flattened" to the point that

the hydrogen and carbon atoms were located in the same plane.

Further decreases in the intermolecular distance led to reaction

to form products HONO, CH20, and CH3O. Thus, these

calculations do not support the hypothesis that "chains" of

molecules are formed through new bonding of intermolecular

methyl—nitro group interactions upon compression ofthe a axis.

However, more detailed analyses corresponding to the specific

interactions that take place in crystals are necessary. The change

in the ordering of the most compressible axes appears to occur

near the pressure at which the equilibrium orientation of the

methyl group has changed to 45° relative to the crystal at low

pressure.9This pressure-induced rotational reorientation provides

another test of our interaction potential, the results of which

we discuss next.

Changes in the rotational motion of the methyl group as a

function of pressure can be observed through RDFs for the

O—H intermolecular bonds. We have already seen from

previous sections that RDFs for O-H intermolecular bonds

indicate an increase of the degree of rotational disorder of the

methyl group as the temperature increases (at ambient pressure).

The RDFs for O—H intermolecular bonds follow a reverse trend

when the pressure is increased (at ambient temperature), as

shown in Figure 3d. In this case the spectrum becomes better

resolved with pressure increases, indicating a significant decrease of the amplitude of methyl rotations around the equilibrium positions.

An illustration of the dynamic behavior of the methyl group

librations with increasing pressure is presented in Figure 9. In

this figure, the time histories of the H,-C]—N2—O3 (i = 5—7)

torsional angles of one of the molecules in the crystal at two
Solid Nitromethane J. Phys. Chem. B,'Vol> 104, No. 35, 2000 Uli

a

200

100

CD

c

<

•a

a

-100

-200

-300

150

a)

^ 100

O)

c

<

o

o

SZ -100

-150

-200

Time (ps)

Figure 9. The time history of the variation of H,—Q—N2-O3 (i =

5,6,7) dihedral angles for one of the nitromethane molecules from

trajectories performed at P = 0.3 GPa (a), and P = 7.0 GPa (b), and

T = 293 K. For visualization purposes some of the angles have been

shifted by ±360°.

different pressures (7/ = 293 K) are represented. In the lowpressure case, the methyl group rotation is slightly hindered, as

methyl group "jumps" from one equilibrium configuration to

another are seen at times 0.5 and 5 ps in the time histories. In

the high-pressure regime the torsional angles oscillate around

well-defined equilibrium values and the jumps between the

equilibrium wells are absent. Also, the equilibrium configuration

for the methyl group changes with pressure, as illustrated by

the time histories in Figure 9. To determine these equilibrium

values we have analyzed the cumulative distributions of the

H,-C1-N2-03 (/ = 5-7) torsional angles for all the molecules

in the simulation box during trajectories executed at different

pressures, T = 293 K. The corresponding data obtained in these

analyses are represented in Figure 4b. When compared to the

distributions of these angles as a function of temperature (see

Figure 4a), the increase of pressure has a totally different effect

on the distribution of torsional angles. In this case there is a

continuous shift of the peak positions with pressure such that

between 0.3 and 5.4 GPa this shift amounts to about 41° while

between 0.3 GPa and 7.0 GP the corresponding variation is

about 50°. These data agree very well with the experimental

findings obtained by Cromer et al.9in which the methyl group

was found to be rotated by about 45° relative to the lowtemperature configuration.

We next explore the basis for the rotational reorientation with

pressure by comparing properties of two crystals, one in which

the methyl groups are rotated and one in which the methyl

groups are in the arrangement corresponding to the lowtemperature, low-pressure configuration. Both crystals have

P2i2i2i symmetry and use the structural parameters corresponding to averages obtained in our simulations. The only difference

between the two models is in the arrangement ofthe hydrogens

of the methyl group. In one model, the methyl group is rotated

by 50° relative to the low-temperature, low-pressure crystal. In

the second model, the methyl group is in the same arrangement

as the low-temperature, low-pressure crystal. We first compared

the total potential energy for the two systems using structural

parameters for T = 4.2 K at 0 atm pressure. The model in which

the methyl group is not rotated has a total potential energy that

is lower by 3.5 kJ/mol per molecule. In that model, the repulsive

N-H contributions to the total potential are greater (by 4.9 kJ/

mol), but the repulsive H-H interactions are lower by 1.4 kJ/

mol. The attractive C-H contribution is larger for the model

in which the methyl group is rotated by 0.9 kJ/mol. The largest

difference in contribution to the total potential is due to the

O-H contributions, which are more attractive by 6.1 kJ/mol

for the low-energy crystal. These combined effects account for

the difference in energy between the two models.

For model crystals using structural parameters corresponding

to a pressure of 7 GPa and temperature of 293 K, the model in

which the methyl groups are rotated by 50° has a total potential

energy that is lower by 35 kJ/mol per molecule than the model

in which the methyl groups are not rotated. The energy

difference (per molecule) is due mainly to the H-N, O—H, and

H-H interactions. For the model in which the methyl group is

rotated by 50°, there is a net decrease of 17.0 kJ/mol for the

N—H repulsive interactions, a net increase of 1.1 kJ/mol for

the C-H attractive interactions, and a net increase of 26.8 kJ/

mol for the O-H attractive interactions. H—H repulsions are

greater in the rotated crystal model by 10 kJ/mol. For both

pressures, the largest differences in the total potential energy

are due to the contributions of O—H and N-H interactions.

The results indicate that the methyl group rotation with pressure

allows for the enhancement of the O-H attractions while

reducing the N-H repulsions.

Structural molecular parameters predicted by the model and

experimental values (where available) as a function of pressure

are given Table 9. The most significant changes in molecular

structure with pressure are the compression of the C—N bond

and the distortion of the Ci-N2-04 angle. These trends were

also seen in experiment; however, the effects were more

pronounced. In the experimental results, the C—N bond is

compressed by 0.08 Ä at 3.5 GPa, and the heavy atom angles

have changed by ~4°. In our model, the corresponding C-N

bond compression is only 0.01 Ä, and the heavy atom angles

have not significantly changed.

A comparison of the average fractional coordinates and

orientational Euler parameters for the four molecules in the unit

cell is given in Table 9. All orientational parameters ofthe four

molecules in the unit cell were averaged over time and all unit

cells in the simulation space. We have also generated orientational parameters for the four molecules in the unit cell assuming

perfect P2i2i2i symmetry in the manner described earlier. The

only complete set of fractional coordinates for all atoms in the

Cromer et al.9study corresponded to 3.5 GPa, 293 K, with the

assumption that the C-H bonds and H-H distances are 1.0

and 1.63 Ä, respectively. The simulation results are in good

agreement with these experimental data, indicating that the

model reasonably represents the crystal under moderate pressure.

Comparison of the averages of the four molecules in the unit

cell with those with P2i2i2i symmetry indicate that the space
8418 J. Phys. Chem. B, Vol. 104, No. 35, 2000 Sorescu et al.

TABLE 9: Comparison of the Average Orientational Parameters (fractional coordinates sx, sy, sz and Euler angles 0, 4», *P)

for the Four Molecules in the Unit Cell and of the Internal Geometrical Parameters of Nitromethane Molecule as Obtained

from Trajectory Calculations at Different Pressures; Corresponding Experimental Values (where available) Are Also Indicated"

pressure(GPa)

0.3 0.3 3.5 3.5 3.5 7.0 7.0

parameter NPT-MD ideal expt*-' NPT-MD ideal NPT-MD ideal

sxl 0.3488 0.3453 0.3391 0.3316

syl 0.4284 0.4055 0.4311 0.4304

szl 0.3765 0.3844 0.3855 0.3938

sx2 0.1519 0.1512 0.1547 0.1609 0.1609 0.1684 0.1684

sy2 -0.4243 -0.4284 -0.4055 -0.4300 -0.4311 -0.4295 -0.4304

sz2 -0.1235 -0.1235 -0.1156 -0.1139 -0.1145 -0.1059 -0.1062

sx3 -0.3479 -0.3488 -0.3453 -0.3388 -0.3391 -0.3314 -0.3316

sy3 -0.0687 -0.0716 -0.0945 -0.0688 -0.0689 -0.0693 -0.0696

sz3 0.1233 0.1235 0.1156 0.1149 0.1145 0.1060 0.1062

sx4 -0.1525 -0.1512 -0.1547 -0.1614 -0.1609 -0.1685 -0.1684

sy4 0.0686 0.0716 0.0945 0.0700 0.0689 0.0706 0.0696

sz4 -0.3769 -0.3765 -0.3844 -0.3853 -0.3855 -0.3940 -0.3938

01 91.4 91.3 90.6 89.5

*1 54.7 55.8 56.2 57.0

¥1 245.3 242.9 248.0 249.3

02 -91.4 -91.4 -91.3 -89.9 -90.6 -89.1 -89.5

<J>2 54.7 54.7 55.8 56.3 56.2 57.0 57.0

»P2 245.3 245.3 242.9 247.7 248.0 249.0 249.3

03 88.4 88.6 88.7 89.6 89.4 90.7 90.5

4>3 -54.8 -54.7 -55.8 -56.3 -56.2 -57.1 -57.0

m 65.2 65.3 62.9 67.8 68.0 69.0 69.3

©4 -88.6 -88.6 -88.7 -89.9 -89.4 -90.7 -90.5

#4 -54.6 -54.7 -55.8 -56.2 -56.2 -57.0 -57.0 u/4 65.2 65.3 62.9 67.7 68.0 69.0 69.3

pressure(GPa)

0.3 0.3 3.5 3.5 7.0

parameter expt* NPT-MD expt* NPT-MD NPT-MD

C-N 1.55 1.500 1.47 1.489 1.481

C-H 1.092 1.088 1.084

N-03 1.21 1.231 1.21 1.228 1.224

N-O4 1.20 1.231 1.26 1.228 1.226

HCH 110.5 110.7 110.6 110.8

HCN 108.1 107.9 107.9 107.7

CNO3 117.7 117.7 118.4 117.0 115.9

CNO4 116.2 117.9 119.9 118.8 119.9

O3NO4 125.9 124.1 . 121.6 123.9 123.9

" In the case of molecules 2—4, the ideal orientational parameters obtained based on P2i2]2i space group symmetry operators and the orientation

parameters of molecule 1, positioned in the asymmetric unit, are also indicated. * Experimental values from ref 9. Since the hydrogen atom positions

could not be resolved at 0.3 GPa, the orientational parameters for nitromethane at this pressure are not available for comparison.cUnit cell fractional

values have been shifted to fall within the range from —1/2 to 1/2.

group symmetry is conserved with pressure, as observed in

experiment.9

Also, as opposed to the temperature effects, which tend to

increase the degree of disorder in system, the increase of

pressure has the opposite effect. As can be seen in Figure 3c

the peaks of RDFs for C*-C intermolecular distances shift

toward smaller distance values with increasing pressure, indicating the compression of the material. Moreover, more and more

peaks regain their individuality with pressure increases, indicating an increase in the degree ofrotational and translational order.

V. Conclusions

"We have developed a classical potential for simulation of

solid nitromethane containing both intra- and intermolecular

potential terms. The parametrization of the intramolecular

potential has been done on the basis of the results of ab initio

calculations performed on the isolated nitromethane molecule

at the B3LYP/6-31G* level. Both structural geometrical parameters as well as data about the vibrational frequencies and

their eigenvectors have been used in the fitting procedure. The

intermolecular potential used in these calculations is of the

Buckingham 6-exp form and was previously developed for RDX

crystals1and showed to be transferable to 30 nitramine crystals4

and to 51 other non-nitxamine crystals.5

Molecular packing calculations using the proposed set of

intermolecular parameters and the set of HF charges indicate

an accurate prediction of crystallographic parameters with

deviations less than 1.21% for the lattice edges. Additionally,

the predicted lattice energy of nitromethane is in acceptable

agreement with previous theoretical estimations.34

The tests of this potential in NPT-MD simulations indicate

that the prediction of the crystallographic parameters is also

well reproduced as a function of temperature with deviations

between 1.7 and 3.6% from the experimental data7-

9

for the

temperature range 4.2—228 K. Moreover, the present potential

is able to predict changes of the crystallographic parameters

with pressure similar to those observed experimentally.9 A

similar good agreement is found for the bulk modulus and its

pressure derivative. At zero pressure the predicted bulk modulus

is Bo = 6.78 GPa while the corresponding experimental value

is 7.0 GPa.9
Solid Nitromethane J. Phys. Chem. B, Vol. 104, No. 35, 2000 8419

The results ofMD and NPT-MD simulations performed using

the rigid-body approximation indicate that the predicted values

in these cases are similar to those obtained from the flexible

model particularly for pressures up to about 5 GPa.

Both RDFs as well as the analyses of the HCNO dihedral

angle distributions during trajectories calculations obtained using

the flexible model indicate two distinctive regimes for methyl

group rotational dynamics. There is an increase of the degree

of rotation of the methyl group with temperature increase (at

ambient pressure) without significant changes of the average

equilibrium positions. On the other hand, at ambient temperature,

the increase of pressure causes a continuous shift ofthe average

equilibrium positions of the H atoms relative to the C-NO2

plane. A comparison of the contributions of the individual

intermolecular interactions to the total potential energies for

crystals in which the methyl groups are in the low-temperature,

low-pressure configuration or arranged in the high-pressure

configuration show that at high pressures, the O-H attractive

interactions are enhanced upon methyl group rotation in the

crystal while the N—H repulsive interactions are decreased; the

remaining intermolecularinteractions do not change significantly

between the two models.

The success of the present potential model to describe the

prototypical explosive, nitromethane, conjugated with the

performances ofthe present intermolecular potential to describe

a large number ofimportant energetic materials, provides further

incentives to develop these models and to extend their applications to more subtle effects such as the energy transfer37'

38 and

reactions in condensed phases.

Acknowledgment We are pleased to acknowledge many

inspiring and helpful discussions with Dr. S. F. Trevino. This

work was supported by the Strategic Environmental Research

and Development Program (SERDP). D.L.T. gratefully acknowledges support by the U.S. Army Research Office under

Grant DAAG55-98-1-0089.

References and Notes

Rice, B. M.; Thompson, D. L. J. Phys. Chem. 1997,

Rice, B. M.; Thompson, D. L. J. Phys. Chem. 1998,

Rice, B. M; Thompson, D. L. J. Phys. Chem. 1998,

Rice, B. M.; Thompson, D. L. J. Phys. Chem. 1998,

Rice, B. M.; Thompson, D. L. J. Phys. Chem. 1999,

Rice, B. M; Thompson, D. L. J. Phys. Chem. 1999,

Prince, E.; Hubbard, C. R. J. Chem. Phys. 1980,

(1) Sorescu, D. C:

B101, 798.

(2) Sorescu, D. C.

B102, 948.

(3) Sorescu, D. C

B 102, 6692.

(4) Sorescu, D. C:

A 102, 8386.

(5) Sorescu, D. C:

A 103, 989.

(6) Sorescu, D. C.

B 103, 6783.

(7) Trevino, S. F.

73, 2996.

(8) Trevino, S. F.:

(9) Cromer, D. T.

2315.

(10) Grosev, V. M.:

181.

Rymes, W. H. J. Chem. Phys. 1980, 73, 3001.

Ryan, R. R.; Schiferl, D. J. Phys. Chem. 1985,89,

Stelzer, F.; Jocham, D. J. Mol Structure 1999,476,

(11) Jones, W. M; Giauque, W. F. J. Am. Chem. Soc. 1947, 69, 983.

(12) Yarger, F. L.; Olinger, B. J. Chem. Phys. 1986, 85, 1534.

(13) Piermarini, G. J.; Block, S.; Miller, P. /. Phys. Chem. 1989, 93,

457.

(14) Cavagnat, D.; Magert, A.; Vertier, C; Anderson, I. S.; Trevino, S.

F. Phys. Rev. Lett. 1985, 54, 193.

(15) Tannebaum, E.; Myers, R. J; Gwinn, W. D. J. Chem. Phys. 1956,

25, 42.

(16) Remizov, A. B.; Musayakaeva, R. H. Opt. Spektrosk. 1975, 38,

226.

(17) 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.: Pittsburgh, PA, 1995.

(18) Möller, C. M. S. Phys. Rev. 1934, 46, 618.

(19) Rice, B. M.; Thompson, D. L. /. Chem. Phys. 1990, 93, 7986.

(20) (a) Becke, A. D. J. Chem. Phys. 1993,98,5648. (b) Lee, C; Yang,

W.; Parr, R. G. Phys. Rev. 1988, B41, 785.

(21) Hehre, W.; Radom, L.; Schleyer, P. v. R.; Pople, J. A. Ab Initio

Molecular Orbital Theory, Wiley-Interscience: New York, 1986.

(22) Rice, B. M.; Chabalowski, C. F.; Adams, G. F.; Mowrey, R. C;

Page, M. Chem. Phys. Lett. 1991, 184, 335.

(23) Forseman, J. B.; Frisch, A. In Exploring Chemistry with Electronic

Structure Methods; Gaussian, Inc.: Pittsburgh, PA, 1996.

(24) (a) Wells, A. J.; WUson, E. B. /. Chem. Phys. 1941, 9, 314. (b)

Smith, D. C; Pan, C; Nielsen, J. R. J. Chem. Phys. 1950, 18, 706. (c)

Jones, W. J.; Sheppard, N. Proc. R. Soc. London Ser. A 1968, 304, 135.

(d) Trinquecoste, C.; Rey-Lafon, M.; ForeL M.-T. Spectrochim. Acta A

1974 30,813. (e) McKean, D. C; Watt, R. A. J. Mol. Spectrosc. 1976, 61,

184.

(25) Pertsin, A, J.; Kitaigorodsky, A. L The Atom-Atom Potential

Method, Applications to Organic Molecular Solids; Springer-Verlag: Berlin,

1987.

(26) Desiraju, G. R. Crystal Engineering: The Design ofOrganic Solids;

Elsevier Amsterdam, 1989.

(27) Gibson, K. 0.; Scheraga, H. A. J. Phys. Chem. 1995, 99, 3752.

(28) Gibson, K. D.; Scheraga, H. A. LMIN: A Program for Crystal

Packing, QCPE, No. 664.

(29) Melchionna, S.; Ciccotti, G.; Holian, B. L. Mol Physics 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 ofthe Research Councils, Daresbury Laboratory at Daresbury,

Nr. Warrington, 1996.

(31) Allen, M. P.; Tildesley, D. J. Computer Simulation of Liquids;

Oxford University Press: New York, 1989.

(32) Fincham, D. Molecular Simulation 1992, 8, 165.

(33) Fundamentals of Crystallography, Giacovazzo, C, Ed.; Oxford

University Press: New York, 1992.

(34) The enthalpy of sublimation for nitromethane has been communicated by P. Politzer and has been obtained based on the method

described in Politzer, P.; Murray, J. S.; Grice, M. E.; DeSalvo, M.; Miller,

E. M. Mol Phys. 1997, 92, 923.

(35) Murnaghan, F. D. In Finite Deformation ofan Elastic Solid; Dover

Publications: New York, 1951; p 73.

(36) Cook, M. D.; Fellows, J.; Haskins, P. J. In Decomposition,

Combustion and Detonation Chemistry of Energetic Materials; Brill, T.

B., Russell, T. P., Tao, W. C, Wardle, R. B., Eds.; Mater. Res. Soc. Symp.

Proc. 1995, 418, 267-275.

(37) Aubuchon, C. M.; Rector, K. D.; Holmes, W.; Fayer, M. D. Chem.

Phys. Lett. 1999, 200, 84.

(38) Deak, J. C; Iwaki, L. K.; Dlott, D. D. J. Phys. Chem. A 1999,

103, 971.
INTENTIONALLY LEFT BLANK.
NO. OF NO. OF

COPIES ORGANIZATION COPIES ORGANIZATION

2 DEFENSE TECHNICAL

INFORMATION CENTER 1 DIRECTOR

DTIC DDA

8725 JOHN J KINGMAN RD

STE0944

US ARMY RESEARCH LAB

AMSRLDD

2800 POWDER MILL RD

FT BELVOJR VA 22060-6218 ADELPHI MD 20783-1197

1 HQDA

DAMOFDT

400 ARMY PENTAGON

WASHINGTON DC 20310-0460

1 DIRECTOR

US ARMY RESEARCH LAB

AMSRL CI AI R (RECORDS MGMT)

2800 POWDER MILL RD

ADELPHI MD 20783-1145

1 OSD

OUSD(A&T)/ODDDR&E(R) 3 DIRECTOR

RJTREW US ARMY RESEARCH LAB

THE PENTAGON AMSRL CILL

WASHINGTON DC 20301-7100 2800 POWDER MILL RD

ADELPHI MD 20783-1145

1 DPTYCGFORRDA

US ARMY MATERIEL CMD 1 DIRECTOR

AMCRDA US ARMY RESEARCH LAB

5001 EISENHOWER AVE AMSRL a AP

ALEXANDRIA VA 22333-0001 2800 POWDER MILL RD

ADELPHI MD 20783-1197

1 INST FOR ADVNCD TCHNLGY

THE UNTV OF TEXAS AT AUSTIN

PO BOX 202797

AUSTIN TX 78720-2797

4

ABERDEEN PROVING GROUND

DIRUSARL

1 DARPA

B KASPAR

3701 N FAIRFAX DR

ARLINGTON VA 22203-1714

AMSRL a LP (BLDG 305)

1 US MILITARY ACADEMY

MATH SCI CTR OF EXCELLENCE

MADNMATH

MAJHUBER

THAYERHALL

WEST POINT NY 10996-1786

1 DIRECTOR

US ARMY RESEARCH LAB

AMSRLD

* DR SMITH

2800 POWDER MILL RD

ADELPHI MD 20783-1197
NO. OF

COPES ORGANIZATION

ABERDEEN PROVING GROUND

22 DIRUSARL

AMSRLWM

BRINGERS

AMSRLWMBD

BEFORCH

WR ANDERSON

SWBUNTE

CFCHABALOWSKI

ACOHEN

R DANIEL

D DEVYNCK

RAFDFER

BEHOMAN

AJKOTLAR

KLMCNESBY

MMCQUAID

MSMILLER

AWMIZIOLEK

JB MORRIS

R A PESCE-RODRIGUEZ

BMRICE

RCSAUSA

MASCHROEDER

JAVANDERHOFF

AMSRLWM MB

BFINK
REPORT DOCUMENTATION PAGE Fom Approved

OUB No. 0704-0188

kiUcrapof»ngburaWiW»U«c<>«»ct^ol'lnWi«»OT^

gwhftna mi natal*** mm *m mi Hit, nid compUfcig id i»»H»>m ft» eo—aten of hrtonmtton. 9m* eowimnt» r»gM<nB «y« bardwi Mini»» »yntr in el»»»

iirttre Svnlost, MncMrati lor hvumtMon Opsntfona tn4 Rkpom, umtflMon

DtvllHlotlwtv.SuHl1204.Artlnaleil.VA 2220*4302, id to 0» Otlk» ot M«n«o»m»nHnd-ajdojiJP»a»r^n»<lu«lonPro^070«-01OT, W«»nhKlton, DC2MC1

1. AGENCY USE ONLY (LMW AtafkV? 2. REPORT DATE

February 2001

3. REPORT TYPE AND DATES COVERED

Reprint, Jun 1999-ffeb 2000

4. TITLE AND SUBTITLE

Theoretical Studies ofSolid Nitromethane

6. AUTHOR(S)

Dan C. Sorescue, * Betsy M. Rice, and Donald L. Thompson*

5. FUNDING NUMBERS

622618.H80

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES)

U.S. Army Research Laboratory

ATTN: AMSRL-WM-BD

Aberdeen Proving Ground, MD 21005-5066

8. PERFORMING ORGANIZATION

REPORT NUMBER

ARL-RP-14

9. SPONSORING/MONITORING AGENCY NAMES(S) AND ADDRESS(ES) lO-SPONSORING/MONITORING

AGENCY REPORT NUMBER

11. SUPPLEMENTARY NOTES

* Department of Chemistry, Oklahoma State University, Stillwater, OK 74078

A reprint from the Journal ofPhysical Chemistry B, vol. 104, no. 35, pp. 8406-8419.

12a. DISTRIBUTION/AVAILABILITYSTATEMENT

Approved for public release; distribution is unlimited.

12b. DISTRIBUTION CODE

13. ABSTRACTßfaxftnum 200 words)

A classical potential to simulate the dynamics of a nitromethane crystal as a function oftemperature and pressure is

described. The intramolecular part of the potential was taken as superposition of bond stretching, bond bending, and

torsional angles terms. These terms were parametrized on the basis of the geometric and spectroscopic (vibrational

frequencies and eigenvectors) data obtained using ab initio molecular orbital calculations performed at the

B3LYP/6-31G level on an isolated molecule. The intermolecular potential used is ofthe Buckingham 6-exp form plus

charge-charge Coulombic interactions and has been previously developed by us (Sorescu, D. C; Rice, B. M.;

Thompson, D. L. J. Phys. Chem. 1997, B101, 798) to simulate crystals containing nitramine molecules and several

other classes of nitro compounds. The analyses performed using constant pressure and temperature molecular

dynamics simulations and molecular packing calculations indicate that the proposed potential model is able to

reproduce accurately the changes of the structural crystallographic parameters as functions oftemperature or pressure

for the entire range of values investigated. In addition, the calculated bulk modulus of nitromethane was found in

excellent agreement with the corresponding experimental results. Moreover, it was determined that the present

potential predicts correctly an experimentally observed 45° change in methyl group orientation in the high-pressure

regime relative to the low-temperature configuration.

14. SUBJECT TERMS

nitromethane, molecular dynamics, intermolecular interaction potential

15. NUMBER OF PAGES

20

16. PRICE CODE

17. SECURITY CLASSIFICATION

OF REPORT

UNCLASSIFIED

18. SECURITY CLASSIFICATION

OF THIS PAGE

UNCLASSIFIED

19. SECURITY CLASSIFICATION

OF ABSTRACT

UNCLASSIFIED

20. LIMITATION OF ABSTRACT

UL

NSN 7540-01-280-5500 Standard Form 298 (Rev. 2-89)

Prescribed by ANSI Std. 239-18 298-102
INTENTIONALLY LEFT BLANK.
USER EVALUATION SHEET/CHANGE OF ADDRESS

This Laboratory undertakes a continuing effort to improve the quality ofthe reports it publishes. Your comments/answers to

the items/questions below will aid us in our efforts.

j 1. ART. Report Number/Author ARL-RP-14 (POC: Rice) Date ofReport February 2001

2. Date Report Received : ___

3. Does this report satisfy a need? (Comment on purpose, related project, or other area ofinterest for which the report will be

used.) _—

4. Specifically, how is the report being used? (Information source, design data, procedure, source ofideas, etc.).

5. Has the information in this report led to any quantitative savings as far as man-hours or dollars saved, operating costs

avoided, or efficiencies achieved, etc? Ifso, please elaborate.

6. General Comments. What do you think should be changed to improve future reports? (Indicate changes to organization,

technical content, format, etc.)

Organization

CURRENT Name E-mail Name

ADDRESS

Street or P.O. Box No.

City, State, Zip Code

7. Ifindicating a Change ofAddress or Address Correction, please provide the Current or Correct address above and the Old or

Incorrect address below.

Organization

OLD Name

ADDRESS

Street or P.O. Box No.

City, State, Zip Code

(Remove this sheet, fold as indicated, tape closed, and mail.)

(DO NOT STAPLE)
DEPARTMENT OF THE ARMY

OfFICIAL BUSINESS

BUSINESS REPLY MAIL

FIRST CLASS PERMIT NO 0001,APG.MD

POSTAGE WILL BE PAID BY ADDRESSEE

DIRECTOR

US ARMY RESEARCH LABORATORY

ATTN AMSRL WM BD

ABERDEEN PROVING GROUND MD 21005-5066

NO POSTAGE

NECESSARY

IF MAILED

IN THE

UNITED STATES