ARL-TN-0784 ● SEP 2016

US Army Research Laboratory

Molecular Simulations of Shear-Induced 

Dynamics in Nitromethane

by Brad A Steele and James P Larentzos

Approved for public release; distribution is unlimited.
NOTICES

Disclaimers

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 of manufacturer’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.
ARL-TN-0784 ● SEP 2016

US Army Research Laboratory

Molecular Simulations of Shear-Induced 

Dynamics in Nitromethane

by Brad A Steele

University of South Florida, Department of Physics, Tampa, Florida

James P Larentzos

Weapons and Materials Research Directorate, ARL

Approved for public release; distribution is unlimited.

FOR OFFICIAL USE ONLY (delete if not FOUO)
ii

REPORT DOCUMENTATION PAGE Form Approved

OMB No. 0704-0188

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instructions, searching existing data sources, gathering and maintaining the 

data needed, and completing and reviewing the collection information. Send comments regarding this burden estimate or any other aspect of this collection of information, including suggestions for reducing the 

burden, to Department of Defense, Washington Headquarters Services, Directorate for Information Operations and Reports (0704-0188), 1215 Jefferson Davis Highway, Suite 1204, Arlington, VA 22202-4302. 

Respondents should be aware that notwithstanding any other provision of law, no person shall be subject to any penalty for failing to comply with a collection of information if it does not display a currently valid 

OMB control number.

PLEASE DO NOT RETURN YOUR FORM TO THE ABOVE ADDRESS.

1. REPORT DATE (DD-MM-YYYY)

September 2016

2. REPORT TYPE

Technical Note

3. DATES COVERED (From - To)

1 May–31 August 2016

4. TITLE AND SUBTITLE

Molecular Simulations of Shear-Induced Dynamics in Nitromethane

5a. CONTRACT NUMBER

5b. GRANT NUMBER

5c. PROGRAM ELEMENT NUMBER

6. AUTHOR(S)

Brad A Steele and James P Larentzos

5d. PROJECT NUMBER

HIP-FY16-018

5e. TASK NUMBER

5f. WORK UNIT NUMBER

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

US Army Research Laboratory

ATTN: RDRL-WML-B

Aberdeen Proving Ground, MD 21005-5069

8. PERFORMING ORGANIZATION REPORT NUMBER

ARL-TN-0784

9. SPONSORING/MONITORING AGENCY NAME(S) AND ADDRESS(ES) 10. SPONSOR/MONITOR'S ACRONYM(S)

11. SPONSOR/MONITOR'S REPORT NUMBER(S)

12. DISTRIBUTION/AVAILABILITY STATEMENT

Approved for public release; distribution is unlimited.

13. SUPPLEMENTARY NOTES

14. ABSTRACT

The atomistic dynamics of the shear response of single-crystal and bicrystal nitromethane (NM) are simulated using molecular 

dynamics simulations. The atomic interactions are described using a recently optimized ReaxFF-lg potential trained 

specifically for NM. We found that NM transforms to a disordered state upon shear at a pressure of 20 GPa and 298 K. Shear 

simulations of several different orientations show that the maximum shear stress and shear angle where the transformation to 

the disordered state occurs are highly dependent on the crystallographic orientation of NM. The dynamics that occur during 

shear in several different single-crystal orientations are elucidated. Shear simulations in bicrystal NM show more complex 

behavior, and the disordered state typically originates at the grain boundary and then grows. The shear response is simulated at 

pressures and temperatures where NM has been experimentally observed to undergo chemical decomposition when shear 

stress is applied. Although no reactions occur in these simulations, these results shed light on the short time scale reaction 

pathways that NM can take during shear-induced initiation.

15. SUBJECT TERMS

nitromethane, shear, molecular dynamics, LAMMPS, grain boundary

16. SECURITY CLASSIFICATION OF:

17. LIMITATION

OF ABSTRACT

UU

18. NUMBER

OF PAGES

24

19a. NAME OF RESPONSIBLE PERSON

James P Larentzos

a. REPORT

Unclassified

b. ABSTRACT

Unclassified

c. THIS PAGE

Unclassified

19b. TELEPHONE NUMBER (Include area code)

410-306-0809 Standard Form 298 (Rev. 8/98)

Prescribed by ANSI Std. Z39.18
Approved for public release; distribution is unlimited.

iii

Contents

List of Figures iv

List of Tables iv

Acknowledgments v

1. Introduction 1

2. Computational Details 2

3. Results 2

4. Conclusions 8

5. References 10

Appendix. Impact of Summer Research Experience 13

List of Symbols, Abbreviations, and Acronyms 15

Distribution List 16
Approved for public release; distribution is unlimited.

iv

List of Figures

Fig. 1 Schematic of the shear simulations on NM using the single crystal 

〈011�〉 orientation sheared along the xy angle as an example. Initially 

a) each angle of the simulation box is 90° and then b) the xy angle is 

reduced to simulate shearing the material. In b) NM has transformed 

into a disordered state. ...........................................................................3

Fig. 2 The a) shear stress Si (xy, xz, or yz) and b) potential energy as a 

function of shear angle (xy, xz, or yz) during the shear simulation for 

the 〈011�〉 orientation of NM. This orientation is chosen to be displayed 

as a representative example of a large shear stress and large shear 

angle as well as a low shear stress and low shear angle that occurs in 

NM during shear. ...................................................................................4

Fig. 3 Zoomed-in snapshots of the NM molecules for the 〈011�〉 orientation 

sheared along the xy angle. The NM molecules are a) initially 

staggered along the y axis, and then after the crystal is sheared b), the 

NM molecules are rotated to orient themselves in the same direction. .4

Fig. 4 Zoomed-in snapshots of the NM molecules for the 〈011�〉 orientation 

sheared along the yz angle. The NM molecules are initially 

a) staggered along the y-axis, and the adjacent planes of NM 

molecules have approximately the same y coordinate. After shear b), 

the NM molecules are oriented along the same direction, and the 

adjacent plane of NM has shifted along the y axis. ...............................5

Fig. 5 The a) shear stress and b) potential energy given as a function of shear 

angle for a bicrystalline simulation of NM............................................7

Fig. 6 Snapshots during the shear stress along the yz angle with a bicrystal 

GB. A line is shown in a) and b) as a guide to the eye to shown the 

rotation of the grain................................................................................8

List of Tables

Table 1 Summary of the maximum shear stress and the shear angle at 

maximum shear stress for each single crystal orientation of NM 

simulated................................................................................................6

Table 2 Summary of the calculated grain boundary (GB) orientations, 

maximum shear stress, shear angle where the shear stress occurs, and 

GB energy ..............................................................................................7
Approved for public release; distribution is unlimited.

v

Acknowledgments

This work is supported by the DOD High Performance Computing Modernization 

Program under grant number HIP-FY16-018. We would like to thank Shawn 

Coleman for providing his grain boundary builder code, which made most of these 

simulations possible. We also thank Betsy Rice, Tim Jenkins, and William Matson 

for their insightful discussions.
Approved for public release; distribution is unlimited.

vi

INTENTIONALLY LEFT BLANK.
Approved for public release; distribution is unlimited.

1

1. Introduction

Nitromethane (NM) is one of the most widely studied insensitive energetic

molecules.

1–8 The main reason for its popularity is that it is one of the simplest 

nitro-organic molecules containing only one nitro group (NO2) bonded to one 

methyl group (CH3) through a C-N bond. This makes it an ideal candidate to study 

chemical reaction mechanisms associated with conventional explosive initiation 

and subsequent detonation.3,6,7 However, understanding the ultra-fast reaction 

mechanisms for initiation is challenging.

The study of initiation and chemical reaction mechanisms of mechanically stressed 

materials such as shock compression, hydrostatic compression, uniaxial 

compression, or shear is called mechanochemistry. The broad field of 

mechanochemistry applied to energetics covers a large range of conditions and 

interactions that can be potentially exploited for advanced energy release 

applications. In addition, mechanochemistry is also relevant for nonenergetics for 

phase transitions and chemical reactions that occur when mechanically stressed.

NM is a liquid at ambient pressure, but when compressed to about 0.3 GPa, it

crystallizes in the P212121 space group.1 When hydrostatically compressed to 

28 GPa, NM slowly decomposes to what is reported to be a 3-D disordered 

polymer.2 In contrast, when NM is compressed to 28 GPa and then sheared, Raman 

studies show that the decomposition is sudden and explosive.5 In addition, 

structural modifications are found to occur at a lower hydrostatic pressure when 

sheared in comparison to hydrostatic pressure alone.5 One model that has been 

proposed to explain the explosive decomposition of NM under shear stress 

examines the orientation-dependent steric hindrance to shear flow at the molecular 

level.4 It was found that certain orientations of NM are highly sterically hindered 

and are therefore more sensitive to chemical reactions and initiation than others. 

However, this is a static model that does not capture the dynamics that occur in the 

crystal during shear.4

In this work we examine the dynamics that occur during shear of NM by performing 

classical molecular dynamics simulations. The crystallographic orientation 

dependence on the shear response of NM is calculated. In addition, since there is a 

strong correlation between grain boundaries (GBs) and chemical initiation in 

energetic materials, bicrystalline NM with a GB is also simulated.
Approved for public release; distribution is unlimited.

2

2. Computational Details

Molecular dynamics simulations are performed using the LAMMPS (Large-scale

Atomic/Molecular Massively Parallel Simulator) simulation package.9 The 

interatomic interactions are described by the ReaxFF-lg potential recently created 

using the MOES (Multiple Objective Evolutionary Strategies) algorithm10,11 and 

refitted to NM.10–13 This potential shows good agreement to the experimental 

equation of state and melting lines for NM. Fully 3-D periodic crystals of NM in 

various orientations are generated using the GB builder code developed by Shawn 

Coleman (US Army Research Laboratory, Weapons and Materials Research 

Directorate) and initially relaxed in the NPT (constant number of atoms, pressure, 

and temperature) ensemble at 20 GPa and 298 K. Special care is taken to ensure the 

crystal structure generated by the code is equivalent to the original crystal structure. 

The crystals are then sheared in the microcanonical ensemble with a constant strain 

rate equal to 0.0008 ps-1 over a 100-ps time period. The strain rate is defined as 

offset divided by the length, where the length is the initial box length in the direction 

perpendicular to shear (e.g., the y box length for xy deformation), and the offset is 

the change in the direction in the shear direction (e.g., the change in the x box length 

for xy deformation). Each crystal is sheared along all 3 angles in the simulation box 

(xy, xz, and yz). The effect the strain rate has on the results is analyzed by 

performing simulations with a strain rate an order of magnitude higher and lower 

than 0.0008, and no significant deviations are found. Each orientation has a 

different number of atoms in the unit cell to ensure 3-D periodicity; however, all 

simulations use at least 20,000 atoms in a simulation box bigger than 35 × 35 ×

35 Å3

. The number of atoms in the cell is also varied to ensure no significant size 

dependence on the results. Several hundred bicrystalline structures of NM are

constructed using similar orientations used in the single crystal simulations and are 

initially relaxed at 20 GPa and 298 K. The bicrystalline structures are 2-D periodic 

in the xz plane while the size of the grain along the y axis is 100 Å (both grains are 

100 Å). Simulations with 50 Å and 300 Å are found to give similar results.

3. Results

A schematic for the shear simulations is shown in Fig. 1a and 1b for orientation 

〈01�1�〉 deformed along the xy angle. The 〈011�〉 orientation is directed along the z

axis. Since NM has an orthorhombic unit cell, each cell angle is initially 90°. As 

the crystal is sheared, the xy angle decreases while the atoms move simultaneously 

according to Newton’s laws to minimize the energy and forces.
Approved for public release; distribution is unlimited.

3

Fig. 1 Schematic of the shear simulations on NM using the single crystal 〈0000�〉 orientation 

sheared along the xy angle as an example. Initially a) each angle of the simulation box is 90° 

and then b) the xy angle is reduced to simulate shearing the material. In b) NM has 

transformed into a disordered state.

The response of NM under shear is highly dependent on the crystallographic 

orientation. In Fig. 2a the shear stress is given as a function of shear angle for the 

〈011�〉 orientation for the 3 shear directions. This orientation is representative of the 

results found for many other orientations that are simulated. The shear stress profile 

is different for all 3 shear directions, but in all cases the shear stress builds up until 

a maximum shear stress is reached. After the maximum shear stress is reached, the 

material suddenly transforms into a disordered state as shown in Fig. 1b and the 

shear stress suddenly drops (see Fig. 2a for shear along the xy angle). In addition 

to the drop in the shear stress, the transformation to the disordered state is 

accompanied with a drop in potential energy as shown in Fig. 2b, which is most 

easily seen for the shear of the xy angle at about 65°. The drop in potential energy 

corresponds to a sudden increase in kinetic energy. Nevertheless, the rapid increase 

in kinetic energy does not cause reactions to occur. The magnitude of the maximum 

shear stress, the shear angle, and the magnitude of the potential energy drop are all 

highly dependent on the crystallographic orientation and shear direction.
Approved for public release; distribution is unlimited.

4

Fig. 2 The a) shear stress Si (xy, xz, or yz) and b) potential energy as a function of shear angle 

(xy, xz, or yz) during the shear simulation for the 〈0000�〉 orientation of NM. This orientation 

is chosen to be displayed as a representative example of a large shear stress and large shear 

angle as well as a low shear stress and low shear angle that occurs in NM during shear.

The orientation dependence can be attributed to the atomistic dynamics that occur 

in the crystal during shear. For shear along the xy angle, the NM molecules are 

initially oriented so that they are staggered with respect to one another along the 

y-axis as shown in Fig. 3a. As the crystal is sheared, the NM molecules rotate so 

that they are roughly oriented in the same direction, as can be seen in Fig. 3b at 

65.77°, which is just before the transformation to the disordered state. The rotation 

of the molecules is what allows the crystal to be sheared to such a large angle, as 

shown in Fig. 2a. Once the molecules can no longer rotate, a large shear stress of 

about 4 GPa builds up, and then the transformation to the disordered state occurs.

Fig. 3 Zoomed-in snapshots of the NM molecules for the 〈0000�〉 orientation sheared along the 

xy angle. The NM molecules are a) initially staggered along the y axis, and then after the 

crystal is sheared b), the NM molecules are rotated to orient themselves in the same direction.

In contrast, when the 〈011�〉 orientation is sheared along the yz angle, only a small 

shear stress is built up (about 2 GPa), as shown in Fig. 2a. This is because in 

addition to the rotation of the NM molecules, the planes of the NM molecules also 
Approved for public release; distribution is unlimited.

5

shift in the shear direction, reducing the shear stress. The initial configuration of 

the NM molecules is displayed in Fig. 4a, where the NM molecules are once again 

staggered along the y axis. Notice that there is some empty space between the 

“planes” of NM molecules in the xy plane in Fig. 4a. This allows the top plane of 

NM molecules to shift with respect to the bottom plane during shear. This can be 

seen by comparing the relative positions of the NM molecules in Fig 4a with those 

in Fig. 4b, where the y-coordinate of adjacent NM molecules along the z axis are 

no longer equivalent. 

Fig. 4 Zoomed-in snapshots of the NM molecules for the 〈0000�〉 orientation sheared along the 

yz angle. The NM molecules are initially a) staggered along the y-axis, and the adjacent planes 

of NM molecules have approximately the same y coordinate. After shear b), the NM molecules 

are oriented along the same direction, and the adjacent plane of NM has shifted along the y

axis.

For all the orientations simulated, the maximum shear stress ranges between 1.68

and 5.24 GPa, and the shear angle ranges from 64.3° to 86.7° as shown in Table 1. 

The 〈011�〉 orientation given as an example has a range of 2.35–4.48 GPa and 

65.5°–81.5° for the 3 shear directions, so it is a good representative example of the 

variety of shear dynamics in NM. The orientations that can be sheared to the lowest 

angles likely have a significant amount of rotation in the NM molecules, such as 

what is observed in the 〈011�〉 orientation. In contrast, the orientations that reach the 

maximum stress near 90° likely do not allow for any rotation of the molecules. At 

this time it is not clear whether initiation is preferred for orientations with the 

highest/lowest shear stress or highest shear angle, or if another parameter such as 

maximum temperature controls the sensitivity in each orientation.
Approved for public release; distribution is unlimited.

6

Table 1 Summary of the maximum shear stress and the shear angle at maximum shear stress 

for each single crystal orientation of NM simulated

Orientation 

〈hhhhhh〉

Max |Sxy| 

(GPa)

Shear angle

(°)

Max |Sxz| 

(GPa)

Shear angle

(°)

Max |Syz| 

(GPa)

Shear 

angle

(°)

〈01�1�〉 4.21 83.1 5.24 80.8 4.03 81.5

〈101�〉 2.15 69.1 2.70 85.2 4.67 82.1

〈111�〉 3.73 83.3 3.34 82.9 4.54 82.1

〈110〉 3.03 66.9 4.53 80.8 2.98 85.2

〈010〉 3.92 64.7 3.46 85.1 1.95 72.9

〈11�1〉 4.06 82.6 3.03 86.7 2.96 77.1

〈111〉 4.08 82.7 3.27 76.1 3.01 86.7

〈11�0〉 3.43 68.4 3.46 86.1 2.85 85.1

〈011�〉 4.38 65.5 4.48 81.5 2.35 71.6

〈101〉 2.19 64.3 2.79 84.9 3.51 78.3

〈011〉 3.07 67.4 2.94 86.6 1.68 72.9

〈001〉 2.60 67.5 2.77 72.7 4.27 64.9

〈100〉 4.26 64.6 3.59 66.2 2.54 64.2

In order to include GBsinto the simulations, the structure at the GB interface should 

be similar to what is present in experiment; however the experimental GB structure 

of NM is currently unknown. To generate realistic GB structures, the following 

method is used. Several hundred GB structures are generated by creating 3 GBs 

with different misorientation angles (along z axis) for grains with orientations 

〈h1kk1ll1〉/〈h2kk2ll2〉 (along the y axis) where h, k and l range from –1 to 1. A total of 

273 GBs are created; however, some of them have the same 〈hkk 〉 orientation and 

therefore represent a bulk structure with a defect at the interface rather than a real 

GB. Nevertheless, all the structures are relaxed at 20 GPa and 300 K and ranked 

according to the GB energy. The following formula is used for the GB energy: 

(UUGGGG − UUpppppppppppppp)/2AA, where UGB is the GB energy and Uperfect is the energy of a 

perfect, single crystal, and A is the area of the GB. The shear simulations are 

performed on the 15 best structures with the lowest GB energy.

When a bicrystalline GB is included in the shear simulations, the shear dynamics 

are more complex; however, the overall process is still similar to the single crystal. 

The shear stress still builds up followed by a transformation to a disordered state. 

In Fig. 5a and 5b, the shear stress and potential energy profile during shear are given 

as a function of shear angle for the 〈011〉/〈011〉 orientation along the z axis and 

〈0 27 �

14���〉/〈0 27���� 14〉 orientation along the y axis. This corresponds to a 

misorientation angle of 90° about the z axis. This orientation is chosen as an 

example because it has similar properties as many of the GBs that are simulated. 

As shown in Table 2, the maximum shear stresses and the GB energy for this 
Approved for public release; distribution is unlimited.

7

orientation are similar to many of the orientations simulated, aside from a few 

outliers.

Fig. 5 The (a) shear stress and (b) potential energy given as a function of shear angle for a 

bicrystalline simulation of NM

Table 2 Summary of the calculated grain boundary (GB) orientations, maximum shear 

stress, shear angle where the shear stress occurs, and GB energy

Z-axis orientation: 

〈hhhhhh〉/〈hhhhhh〉

Y-axis orientation:

〈hhhhhh〉/〈hhhhhh〉

GB 

energy 

(eV/Å2)

Max. 

|Sxy| 

(GPa)

Angle 

(°)

Max.

|Sxz| 

(GPa)

Angle 

(°)

Max. 

|Sxz| 

(GPa)

Angle 

(°)

〈011〉/〈011〉 〈0 27 14����〉/〈0 27 14〉 0.903 1.56 85.0 1.91 76.6 2.10 70.6

〈011〉/〈011�〉 〈0 27 14����〉/〈0 27 14〉 1.025 1.46 70.6 1.76 85.7 2.11 71.6

〈011〉/〈011�〉 〈0 27 14〉/〈0 27���� 14����〉 0.932 1.33 79.5 1.70 82.5 2.87 60.9

〈11�0〉/〈11�0〉 〈3 2 0〉/〈49���� 33���� 0〉 0.867 1.80 85.3 3.64 73.8 1.64 59.7

〈11�1〉/〈11�1〉 〈8 15 5〉/〈8� 15���� 5�〉 0.819 2.55 77.2 2.48 82.2 2.06 68.7

〈11�1〉/〈111〉 〈8� 15���� 5�〉/〈8� 15 5�〉 1.246 2.09 78.1 4.15 83.1 1.92 69.7

〈11�1〉/〈111�〉 〈8 15 5〉/〈8� 15 5〉 3.111 2.09 80.9 2.55 64.2 1.84 82.9

〈11�1〉/〈111�〉 〈8 15 5〉/〈8 15���� 5�〉 0.921 2.14 77.5 2.21 76.6 2.21 65.0

〈11�1�〉/〈11�1�〉 〈8� 15���� 5〉/〈8 15 5�〉 0.698 1.98 82.6 2.60 82.0 1.98 69.3

〈11�1�〉/〈111�〉 〈8 15 5�〉/〈8 15���� 5�〉 0.987 1.71 85.7 3.76 83.4 2.17 67.5

〈11�1�〉/〈111�〉 〈8 15 5�〉/〈8� 15 5〉 0.773 2.00 81.6 2.34 74.9 1.92 69.5

〈111〉/〈111〉 〈8 15���� 5〉/〈8� 15 5�〉 0.770 1.99 81.8 2.62 81.8 1.90 67.5

〈111〉/〈111�〉 〈8� 15 5�〉/〈8� 15 5〉 0.928 1.85 83.9 2.29 73.5 1.90 71.0

〈111�〉/〈111�〉 〈8� 15 5〉/〈8 15���� 5�〉 0.809 2.39 77.9 2.42 70.8 1.91 68.1

〈011〉/〈011〉 〈0 27���� 14〉/〈0 27���� 14〉 0.449 1.49 85.0 1.73 72.6 2.5 58.4

As shown in Fig. 5a, as this orientation is sheared along the yz angle, the shear 

stress builds up slowly and even drops slightly around 85°, then builds up again. In 

this simulation, the grain is slowly rotating with the shear, as displayed in Fig. 6a 

and 6b. The molecules do not rotate much in the beginning of the simulation, but 

the grain rotates considerably. This can be seen more clearly by examining the 

C-N bonds in NM along the arrow shown in Fig. 6a and 6b.
Approved for public release; distribution is unlimited.

8

Fig. 6 Snapshots during the shear stress along the yz angle with a bicrystal GB. A line is 

shown in a) and b) as a guide to the eye to shown the rotation of the grain.

Also shown in Fig. 6a is that the disordered material originates at the GB interface 

and expands from there as shown in Fig. 6b. However, in stark contrast, when 

sheared along the xz angle, the entire crystal suddenly transforms into the 

disordered state, which causes the shear stress to drop rapidly, as shown in Fig. 5a 

at 85°. Interestingly, the shear stress then rises again. This is because the 

transformation to the disordered state is incomplete, and the material still maintains 

a small amount of order. The transformation to the disordered state is not complete 

until about 75°.

As shown in Table 1, the largest absolute shear stress is found to occur in the 

〈01�1〉/〈111〉 orientation along z and 〈8� 15���� 5�〉/〈8� 15 5�〉 along y sheared along the 

xz axis of 4.15 GPa at an angle of 83.1°. However, in spite of the large build-up of 

shear stress, no reactions are observed in this bicrystalline simulation, which is 

similar to what is found for the single crystal simulations.

4. Conclusions

The short time-scale dynamics during shear of NM have been revealed using 

molecular dynamics simulations. Shear simulations have been performed using 

single crystals with a range of orientations and shear directions as well as bicrystals,
Approved for public release; distribution is unlimited.

9

also with a range of orientations and misorientation angles. Although no chemical 

reactions are observed during the simulations, the dynamics presented are possibly 

associated with the short time-scale reaction pathway that leads to initiation. In this 

sense, these results are critical for understanding the complete mechanism for 

initiation of NM under shear stress.

Since no reactions are observed in the simulations but are observed in experiment,

the future work is mostly geared toward resolving this intriguing discrepancy as 

well as generally increasing the length and time scales, and the range of orientations 

studied. The initiation may be related to the presence of multiple GBs, so 

simulations with polycrystalline NM may be important. Investigating the sizedependence of the results more thoroughly with GBs may also be needed because 

the GB simulations require that the size of the grain be much larger than the size of 

the GB. In addition, initiation may be related to the slip and subsequent momentum 

transfer of grains much larger than what is simulated, so the size-dependence needs 

to be investigated thoroughly. The time scale of the simulations may also be an 

issue because it may be that the time scale for reactions to occur is simply larger 

than the time scales simulated. An analysis of activation energies and Arrhenius 

reaction rates for proposed reactions may yield an estimate of the time scale needed. 

A greater variety of GB orientations can also be simulated, which may yield 

interesting results. All the results presented are at 20 GPa; it would be interesting 

to examine the pressure dependence of the results, since there is a clear pressure 

dependence in the experiment related to the reaction threshold pressure.

This work was conducted as part of the DOD High Performance Computing 

Modernization Program’s (HPCMP’s) Workforce Development initiative and is 

aimed at providing the future science and engineering workforce the computational 

skills and experience necessary to close the gap between the technological 

capability and computational skills necessary to support the DOD’s future 

Warfighter needs. The HPCMP’s 2016 Internship Program was a valuable 

experience for Brad Steele, and the impact of the summer research experience on 

his professional career is discussed in the Appendix. 
Approved for public release; distribution is unlimited.

10

5. References

1. Citroni M, Datchi F, Bini R, Di Vaira M, Pruzan P, Canny B, Schettino V. 

Crystal structure of nitromethane up to the reaction threshold pressure. J Phys

Chem B. 2008;112(4):1095–1103.

2. Citroni M, Bini R, Pagliai M, Cardini G, Schettino V. Nitromethane 

decomposition under high static pressure. J Phys Chem B. 2010;114(29):

9420–9428.

3. Guo F, Cheng X-l, Zhang H. Reactive molecular dynamics simulation of solid 

nitromethane impact on (010) surfaces induced and nonimpact thermal 

decomposition. J Phys Chem A. 2012;116(14):3514–3520.

4. Dick JJ. Orientation-dependent explosion sensitivity of solid nitromethane. J 

Phys Chem. 1993;97(23):6193–6196.

5. Hébert PA, Isambert A, Petitet J-P, Zerr A. Raman spectroscopy study of 

nitromethane in a shear diamond anvil cell. High Pressure Research.

2010;30(1):24–27.

6. Manaa RM, Reed EJ, Fried LE, Galli G, Gygi F. Early chemistry in hot and 

dense nitromethane: molecular dynamics simulations. J Chem Phys.

2004;120(21):10146–10153.

7. Nelson T, Bjorgaard J, Greenfield M, Bolme C, Brown K, McGrane S, Scharff

RJ, Tretiak S. Ultrafast photodissociation dynamics of nitromethane. J Phys 

Chem A. 2016;120(4):519–526.

8. Courtecuisse S, Cansell F, Fabre D, Petitet JP. Phase transitions and chemical 

transformations of nitromethane up to 350 °C and 35 GPa. J Chem Phys.

1995;102(2):968–974.

9. Plimpton S. Fast parallel algorithms for short-range molecular dynamics. J 

Comp Phys. 1995;117(1):1–19.

10. Larentzos JP, Rice BM, Byrd EFC, Weingarten NS, Lill JV. Parameterizing 

complex reactive force fields using multiple objective evolutionary strategies 

(MOES). Part 1: ReaxFF models for cyclotrimethylene trinitramine (RDX) 

and 1,1-diamino-2,2-dinitroethene (FOX-7). J ChemTheory and Computation.

2015;11(2):381–391.
Approved for public release; distribution is unlimited.

11

11. Rice BM, Larentzos JP, Byrd EFC, Weingarten NS. Parameterizing complex 

reactive force fields using multiple objective evolutionary strategies (MOES): 

Part 2: transferability of ReaxFF models to C–H–N–O energetic materials. J

Chemical Theory and Computation. 2015;11(2):392–405.

12. Strachan A, van Duin ACT, Chakraborty D, Dasgupta S, Goddard WA. Shock 

waves in high-energy materials: the initial chemical events in nitramine RDX. 

Physical Review Letters. 2003;91(9):098301.

13. Zhang L, Zybin SV, van Duin ACT, Dasgupta S, Goddard WA, Kober EM.

Carbon cluster formation during thermal decomposition of octahydro-1,3,5,7-

tetranitro-1,3,5,7-tetrazocine and 1,3,5-triamino-2,4,6-trinitrobenzene high 

explosives from ReaxFF reactive molecular dynamics simulations. J Phys 

Chem A. 2009;113(40):10619–10640.
Approved for public release; distribution is unlimited.

12

INTENTIONALLY LEFT BLANK.
Approved for public release; distribution is unlimited.

13

Appendix. Impact of Summer Research Experience
Approved for public release; distribution is unlimited.

14

The following remarks were made by Brad Steele in reference to his experience as 

a student intern at ARL under the 2016 HPCMP Internship Program.

Background

I am a PhD candidate at the University of South Florida, Physics Department, where 

I earned my MS in Applied Physics. I have a significant amount of experience 

performing density functional theory calculations and bash scripting and 

programming for running large amounts of simulations and data analysis. My past 

research experience dealt with predicting and characterizing new kinds of materials 

at high pressures. This involves zero-temperature energy-minimization calculations 

to determine thermodynamic and dynamical stability. I have much less experience 

performing simulations using LAMMPS and analyzing time-dependent dynamical 

results. My goals for participating in the internship were to gain more experience 

using LAMMPS and analyzing dynamical results. Furthermore, I wanted to gain 

experience conducting research at a government research lab to understand the 

techniques and how research challenges are tackled on a government level.

Impact

The HPCMP Internship Program has been incredibly useful for me. I have gained 

a tremendous amount of experience running molecular dynamics simulations and 

analyzing the results. I have also been able to further develop my programming and 

bash scripting skills. More significantly, I have a much better appreciation for what 

molecular dynamics is capable of and what its limitations are. Before this 

internship, I had never used ReaxFF; I now recognize its potential. Furthermore, I 

have a much better understanding of the technical terms used by the energetic 

materials community and what types of properties they are interested in. One thing 

that surprised me about my project was how difficult it is to get the metastable 

molecules to react. I had generally assumed that chemical reactions can be observed 

easily using molecular dynamics, but that is not the case. I have also gained an 

appreciation of what researchers mean when they say a material reacts or 

decomposes under some sort of stress, why the response of the material depends on 

the type of stress applied, and why the rate of the reaction is very important factor 

for energetic materials. On a personal level, this experience is of tremendous value 

to me because I want to work at a government lab when I graduate. It has been very 

good to meet the other researchers and learn about what they are working on and 

how they are trying to solve their research problems. This program also fulfills a 

necessary requirement for my PhD.
Approved for public release; distribution is unlimited.

15

List of Symbols, Abbreviations, and Acronyms

GB grain boundary

HPCMP High Performance Computing Modernization Program

LAMMPS Large-scale Atomic/Molecular Massively Parallel Simulator

MOES Multiple Objective Evolutionary Strategies

NM nitromethane (CH2NO3)
Approved for public release; distribution is unlimited.

16

1 DEFENSE TECHNICAL

(PDF) INFORMATION CTR

DTIC OCA

2 DIRECTOR

(PDF) US ARMY RESEARCH LAB

RDRL CIO L

IMAL HRA MAIL & RECORDS

MGMT

1 GOVT PRINTG OFC

(PDF) A MALHOTRA

2 UNIV OF SOUTH FLA

(PDF) DEPT OF PHYSICS

B STEELE

I OLEYNIK

4 HIGH PERFORM CMPTG

(PDF) MODERNIZATION PROG OFC

S BELL

C DAHL

L DAVIS

D SCHWARTZ

7 DIR USARL

(PDF) RDRL WML B

T JENKINS

J P LARENTZOS

W MATTSON

B RICE

N S WEINGARTEN

N TRIVEDI

RDRL WMM F

S COLEMAN