Coarse-Grain Model Simulations of Nonequilibrium Dynamics in

Jun 5, 2014 - Weapons and Materials Research Directorate, U.S. Army Research ... 6189, Theoretical Chemistry Section, U.S. Naval Research Laboratory, ...
0 downloads 0 Views 656KB Size
Subscriber access provided by Oregon State University

Letter

Coarse-Grain Model Simulations of NonEquilibrium Dynamics in Heterogeneous Materials John Kennedy Brennan, Martin Lisal, Joshua D. Moore, Sergei Izvekov, Igor V Schweigert, and James Peter Larentzos J. Phys. Chem. Lett., Just Accepted Manuscript • Publication Date (Web): 05 Jun 2014 Downloaded from http://pubs.acs.org on June 6, 2014

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry Letters is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Coarse-Grain Model Simulations of Non-Equilibrium Dynamics in Heterogeneous Materials *John K. Brennan1, Martin Lísal2,3, Joshua D. Moore1, Sergei Izvekov1, Igor V. Schweigert4, and James P. Larentzos5 1

Weapons and Materials Research Directorate, U.S. Army Research Laboratory, Aberdeen Proving Ground, MD, USA 2 Laboratory of Chemistry and Physics of Aerosols, Institute of Chemical Process Fundamentals of the ASCR, v. v. i., Prague, Czech Republic 3 Department of Physics, Faculty of Science, J. E. Purkinje University, České Mládeže 8, 400 96 Ústí n. Lab., Czech Republic 4 Code 6189, Theoretical Chemistry Section, U.S. Naval Research Laboratory, Washington, DC, USA 5 Engility Corporation, High Performance Technologies Group at the U.S. Army Research Laboratory, Aberdeen Proving Ground, MD, USA Abstract A suite of computational tools is described for particle-based mesoscale simulations of the non-equilibrium dynamics of energetic solids, including mechanical deformation, phase transitions, and chemical reactivity triggered by shock or thermal loading. The method builds upon our recent advances both in generating coarse-grain models under high strains and in developing a variant of dissipative particle dynamics (DPD) that includes chemical reactions. To describe chemical reactivity, a coarse-grain particle equation-of-state was introduced into the constant-energy DPD variant that rigorously treats complex chemical reactions and the associated chemical energy release. As illustration of these developments, we present simulations of shock compression of an RDX crystal and its thermal decomposition under high temperatures. We also discuss our current efforts towards a highly-scalable domaindecomposition implementation that extends applicability to micron-size simulations. With appropriate parameterization, the method is applicable to other materials whose dynamic response is driven by microstructural heterogeneities.

Keywords: coarse-grain, simulation, modeling, dissipative particle dynamics, non-equilibrium, dynamics John K. Brennan: [email protected] (corresponding author) Martin Lísal: [email protected] Joshua D. Moore: [email protected] Sergei Izvekov: [email protected] Igor V. Schweigert: [email protected] James P. Larentzos: [email protected] 1

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 12

To date, computational resources and algorithms limit simulations of fully-atomistic models to nanometer length scales, where simulations of dense materials such as molecular solids are particularly costly. Therefore, computationally-reasonable particle-based simulations of material behavior governed by micro- and mesoscale structural heterogeneities require implementation of coarse-grain (CG) models and methodologies.1,2 In this Letter, we describe our development of a general modeling framework that is capable of simulating the dynamic response of CG models to various stimuli, such as mechanical or thermal shock. A number of challenges exist for building the framework, including the development of CG models and methods that can capture the following known thermo-mechanical responses: (i) phase transitions; (ii) structural rearrangements and mechanical deformation; and (iii) chemical reactivity. The CG models are built by utilizing the force-matching based multi-scale coarse-graining method,3 whereas the CG methodology framework is built upon the well-established dissipative particle dynamics (DPD) approach and includes: (a) a generalized expression for the particle internal energy in the isoenergetic (DPD-E) and isoenthalpic (DPD-H) methods;4-7 (b) both parallel and perpendicular dissipative forces acting on the CG particles; (c) chemical reactivity involving both inter- and intra-particle interactions; and (d) a generalized and efficient numerical integration scheme. We employ systematic parameterization from higher resolution models for all components, although higher resolution model input is not a requirement. Implementation of the tool suite into a highly-scalable domain decomposition code is ongoing.8 In this Letter, demonstrations of our modeling framework are presented for non-reactive shock of a one-site coarse-grain model of RDXa (CG-RDX), as well as the reactive response of this model to heating. Our CG models are built by applying an efficient and systematic bottom-up strategy for constructing effective pairwise and central CG force-fields.9,10 The method, termed multi-scale coarse-graining (MS-CG),3 requires applying least-squares approximations to the many-body forces as derived from the system free-energy potential. In this work, we have refined a pairwise-particle CG potential that has been previously built using the MS-CG method, where a single RDX molecule was mapped onto a one-site CG model.11 In the original work, the transferability of the model to higher pressures was improved by incorporating an explicit dependence on the local particle density, where direct mapping of the pressure-volume behavior of the atomistic model12,13 to the CG model was performed. Although no longer exclusively a bottom-up derived model, the resulting density-dependent potential of the original model was able to reasonably reproduce several atomistic model properties, including the crystal structure, elastic and vibrational properties, pressure-volume isotherms, and the melting point at ambient pressure. In this work, both the equilibrium and non-equilibrium properties were further assessed, where we have made the following improvements and refinements. First, the CG-RDX potential reported earlier was modified to eliminate spurious crystalline phases11 at elevated temperatures and pressures by removing a short-ranged interaction artifact. In the earlier model, sufficient sampling of the atomistic model trajectories could not be achieved in the crystalline phase, thus sampling was restricted to the molten phase, which resulted in a softer repulsive core characteristic of the liquid and low-temperature amorphous phases. Increasing the repulsion of the short-ranged separations significantly increased the transferability of the CG-RDX model for the crystalline phase at high temperature and under mechanical loading. Critical improvement in a

hexahydro-1,3,5-trinitro-s-triazine, C3H6N6O6

2

ACS Paragon Plus Environment

Page 3 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

the energy conservation of the model also has been made by: (a) fitting the force-field to a Chebyshev polynomial, thereby smoothing discontinuities; (b) modifying the local density calculation scheme; and (c) re-casting the functional form of the forces into separate two-body and many-body contributions with the appropriate normalization. For this latter effort, we have separated the conservative force for the CG-RDX model into two components:  = ,  + ,  ,  , 

(1)

where, ,  is a pair-wise force, which does not depend upon the local density ( ); and

,  ,  ,  is a local-density dependent force-field, which follows the generalized manybody DPD form developed by Pagonabarraga and Frenkel,14 and subsequently used by others.1521 Further details of the CG-RDX model can be found in the Supporting Information (SI). For the simulation of materials that preclude the use of atomistic models, careful consideration must be taken for the choice of a CG model, but as critically, for the choice of an appropriate CG methodology. To determine properties other than static (equilibrium) properties, the molecular dynamics (MD) approach is not adequate due to the well-known speed-up of the dynamics of CG models compared to their atomistic counterparts, which is a direct consequence of losing molecular degrees-of-freedom (d.o.f.) during coarse graining.23,24 For our purposes, the DPD-E method4-7 is particularly suitable for non-equilibrium simulation scenarios and thermallyvariant conditions. DPD-E incorporates explicit treatment of the coarse-grained d.o.f. through both the dissipative forces and the particle internal energy term (ui), where the latter provides a computational means of ensuring energy conservation during the simulation. Moreover, ui provides an avenue to compensate for the coarse-grained d.o.f., which is essential for accurately recovering particular properties of the atomistic model (e.g., the shock behavior, as shown below in Figure 1). A characterizing feature of the DPD methodology is the dissipative forces that act between particles, which provide a means of accurately depicting the dynamics of CG models. We have developed several refinements of the dissipative contributions and incorporated them into this work, including: (i) a bottom-up parameterization of the dissipative force parameters;25 and (ii) within the DPD-E framework, incorporation of multi-directional dissipative interactions26,27 that are both parallel and perpendicular to the interparticle separation axis, which attempt to capture the coarse-grained d.o.f. that contribute to the molecular shape or polarity. (See SI for additional details.) Recognizing a limitation in the available DPD numerical integration schemes, our group recently developed an efficient integration algorithm for all DPD variants4 using Shardlow-like splitting algorithms (SSA).28 The SSA involves splitting the deterministic and stochastic dynamics, which allows straightforward extensions to all the variants and does not require modification to the existing deterministic integration schemes. Compared to the traditional DPD integrators, the SSA allows for larger time steps making simulations of realistic CG models now possible for the DPD-E and DPD-H variants. For the large-scale simulations reported in this work to be viable, we implemented the SSA integration scheme into the open-source LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) simulation package,29 utilizing the highly-scalable domain decomposition structure. The recursive nature of determining the i − j particle pair interactions in the SSA required special consideration when applying the scheme to 3

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 12

the LAMMPS framework, which involves spatial-decomposition of the simulation domain.8 While additional computational overhead was required during the simulation, a significant reduction in the overall time-to-solution O(10-102) was achieved. As an initial assessment of our tool set, we simulated crystalline RDX under shock compression using a fully-atomistic model12, 13 and our one-site CG-RDX model. Two types of simulations were performed: equilibrium simulations of the Hugoniot properties30 (Figure 1a) and non-equilibrium simulations of piston-supported shock compression (Figure 1b). Shown in Figure 1a is a comparison of the equilibrated temperature and relative volume of unreacted, shock-compressed RDX crystal for the models. MD simulations of the CG model (CG MD) were unable to reproduce the fully-atomistic model (Atomistic MD) Hugoniot temperature due to the absence of the CG intra-molecular d.o.f., which is accounted for in the DPD-E simulation via the CG particle equation-of-state (CG-EOS). For example, at 9 GPa, the temperature in the CG MD simulation was greater than 2000 K compared to 410 K31 in atomistic MD. In contrast, the CG DPD-E simulations, which incorporated explicit treatment of the CG d.o.f., was in excellent agreement with the atomistic model at pressures up to 2 GPa. Deviation at higher pressures (approximately 15% at 9 GPa) was attributed to inaccuracy of the implemented CG-EOS, where refinements of the CG-EOS are expected to reduce the deviation. In comparison to the atomistic MD result, the relative volumes for the CG DPD-E simulations were in better agreement; however, differences were less pronounced than for the temperature plane. This was evidence that the pressure-volume behavior of the CG-RDX model was not strongly dependent on temperature within the observed temperature range. Analogous behavior was observed in non-equilibrium simulations of piston-supported shock compression of an RDX crystal. Shown in Figure 1b are the temperature and density profiles along the shock direction at 15 ps after the crystal and piston impact. Since a non-reactive atomistic force-field was used as the reference, the piston velocity was kept low to avoid higher temperatures at which reactions may occur within the simulated time scales. In comparison to the atomistic MD result, the CG MD simulation significantly overestimated the temperature rise at and behind the shock front. However, the CG DPD-E simulation more accurately reproduced the temperature rise at the shock front due to the inclusion of the CG-EOS, which provides a mechanism to transfer the mechanical energy generated by the flyer-plate impact into the CG d.o.f. via the heat and momentum exchange between the particle internal energies. Similar to the volume plane of the Hugoniot, the CG DPD-E simulation resulted in a density profile that was in better agreement with the atomistic MD result.

4

ACS Paragon Plus Environment

Page 5 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

(a)

(b)

Figure 1: (a) Comparison of the shock Hugoniot properties of pure crystalline RDX: (Atomistic MD) fullyatomistic model using MD; (CG MD) one-site CG-RDX model using MD; and (CG DPD-E). Hugoniot properties are locus of points corresponding to the final shocked state originating from an initial un-shocked state (here P0 = 0 GPa, T0 = 300 K). Relative volume is the ratio of the shocked state volume to the initial un-shocked state volume. (b) For the same simulation approaches shown in (a), comparison of the non-equilibrium response due to mechanical shock generated in the [100] directionb by a 120a x 20b x 20c unit cell sample impacting a stationary piston comprised of rigid CG-RDX molecules at  = 0.5 km/s. Kinetic temperature and density profiles are snapshots taken 15.0 ps after impact. Relative density is the ratio of the shocked state density to the initial un-shocked state density. (See SI for further simulation details.)

A key feature of the DPD methodology is its ability to incorporate chemical reactions and thus describe the coupling between the material reactivity and its mechanical and thermodynamic properties. This feature is particularly appealing for simulating energetic materials under loading, wherein even a weak insult can trigger chemical reactions leading to material decomposition and a rapid release of chemical energy. The individual chemical reactions responsible for the energy release have characteristic time scales ranging from ns to µs, but remain much slower than the thermal motions of atoms (fs). A fully-atomistic model description of these reactions would require simulation over these disparate time scales, which currently is not computationally feasible. Furthermore, an atomistic simulation of these reactions must rely on interatomic potentials (e.g., see Refs. 33, 34) that properly describe the initial material, the reaction products, and the relevant transition structures, entirely within a single parameterization. DPD provides a framework for “coarse-graining” the explicit reaction dynamics via a reaction b

Coordinate system used for the crystallographic directions is the same as found in Ref. 32.

5

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 12

rate model, thus circumventing both the time scale discrepancy and the need for reactive potentials. We have developed an extension of DPD that incorporates chemical reactivity35,36 (notated DPD-RX) by effectively treating a CG particle as a continuous-stirred tank reactor,37 termed coarse-grain reactor (CG-reactor). Here CG-reactor infers the computational means of depicting the reactivity occurring within the local volume of the CG particle, whereby the species diffusion is assumed to be instantaneous at each time step. The chemical reactivity that results in composition variability of the CG-reactor is monitored by an extent-of-reaction variable assigned to each CG particle for each reaction in the reaction mechanism. In DPD-RX, the extent-ofreaction variable tracks the progress of material transformation according to the prescribed reaction rate model, similar to the manner in which chemical reactions are tracked within each volume cell in reactive flow simulations.38 Reaction rate models can entail reaction mechanisms that are single step or multi-step, and can be developed in a top-down fashion, by postulating a reduced reaction mechanism and calibrating it against experimental measurements, or in a bottom-up fashion, by developing a detailed mechanism from kinetic measurements and ab initio calculations. For the work presented here, the DPD-RX method has been generalized and extended from the original work of Maillet and co-workers35,36 to allow for additional flexibility with respect to the model chemistries. Most notably, a CG-EOS is introduced that couples the CG particle internal energy and internal partition functions of the reacting species (expressed via their standard molar chemical potentials). We derived the CG-EOS by relating the thermodynamic total internal energy of the reacting mixture and the DPD-E total energy.39 The internal energy of each CG-reactor is computed according to Hess’s law, accounting for the chemical energy release using the current concentrations of the reacting species and their standard enthalpies of formation and temperature-dependent heat capacities. (See SI for more details of the partitionfunction based CG-EOS.) As illustration of the DPD-RX method for treating chemical reactions in an energetic solid, for example, triggered by a thermal load, we developed a four-step reaction rate model describing the decomposition of RDX and the combustion of its decomposition products. The model was derived by combining a single-step rate of RDX decomposition above the melting temperature,40 and the detailed mechanism of nitramine combustion in the gas phase.41 The combined mechanism contained more than 30 chemical species, where a direct implementation would require the potentials and mixing rules for each species. Instead, the combined mechanism was reduced to a four-step rate model consisting of a unimolecular, irreversible reaction: R1: RDX→ 3HCN + 3/2(NO2 + NO + H2O)

(2)

and three bimolecular, irreversible reactions: R2: HCN + NO2 → NO + ½(N2 + H2) + CO R3: HCN + NO → CO + N2 + ½H2

(3)

R4: NO + CO → ½N2 + CO2

6

ACS Paragon Plus Environment

Page 7 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

This mechanism was heuristically chosen based on comparison to perfectly-stirred reactor simulations of the original mechanism at temperatures from 1000 to 3000 K, where the corresponding global rates were fitted to reproduce the speciation profiles. The complete set of rate laws and reaction rate constants are given in the SI. For the application of the DPD-RX method presented here, the chemical character of any particle can vary from pure RDX to a mixture of product gases, and states in between. The interaction potential was built by first specifying the end state interaction potentials, where at any time, the total configurational energy of a particle is given by the appropriately weighted sum of each potential according to their number fractions. The cross-interactions between the RDX and gas-mixture components of a pair of particles are then approximated in a manner similar to sitesite interaction potentials, with each of the end states taken as a site and the cross-interactions between sites appropriately weighted. The interaction between particle i with an RDX fraction, niRDX , and a product species fraction, 1 − niRDX , and particle j with an RDX fraction, n RDX , and j

(

(

)

)

, is approximated as a product species fraction, 1 − n RDX j

(

)

uijCG = niRDX n RDX uijRDX -RDX + niRDX 1 − n RDX uijRDX -psf j j +

(1 − n )n RDX i

RDX j

uijpsf -RDX +

(1 − n )(1 − n )u RDX i

RDX j

(4) psf - psf ij

where uijRDX is the CG interaction potential of RDX, uijpsf is the interaction potential for the product gas mixture, and uijRDX- psf and uijpsf - RDX are the cross-interactions. The product mixture is treated as a pseudo-fluid mixture (psf) using the van der Waals one-fluid approximation,42 where each gas species within the mixture is modeled using the exponential-6 potential.43 In applying the one-fluid approximation, the exponential-6 potential parameters of the pseudo-fluid depend on the composition of the mixture, where standard mixing rules are applied.44 Corresponding expressions for the forces are analogous to the interaction potential energy expressions. Note that the CG interaction potential defined here is not unique and other treatments are possible. The results of a constant-volume and constant-energy DPD-RX simulation of RDX decomposition are shown in Figure 2. An ideal RDX crystal at its ambient density was represented by 3000 particles confined in a 200 Å x 57.9 Å x 52.9 Å periodic simulation box, equilibrated at 300 K. At the start of the simulation, a 50 Å wide interior slab of the crystal was instantaneously heated to 3000 K by setting the particle internal temperatures to 3000 K and by selecting the particle velocities from the corresponding Maxwell-Boltzmann distribution. The ensuing DPD-RX dynamics was followed for approximately 3 ns using a variable time step ranging from 0.005-0.5 fs. At 3000 K, the rate of RDX dissociation [Eq. (2)] is on the scale of 1014 s-1, thus, in the heated central slab, all RDX decomposed to HCN and other intermediates within the first 10 fs. This global dissociation reaction is slightly exothermic and the temperature in the central region rose above the initial 3000 K (Figure 2a). Between 1 and 10 ps, as the heat is dissipated into the adjacent slabs, the crystal started to melt, as indicated by some disruption of the repeating molecular layers (Figure 2b). The phase change was followed by RDX dissociation in the adjacent layers (Figure 2c). Between 10 and 100 ps, as the heat continued to dissipate into the adjacent slabs, this portion of the sample continued to melt and dissociate, while the central slab 7

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 12

cooled down. Interestingly, as this occurred, the fluid density in the central region started to rise. This was attributed to the thermal expansion of molten RDX in the adjacent slabs compressing the mixture of the centrally-located dissociation intermediates (HCN, NO2, etc.). At 1 ns, the crystal completely melted and the temperature started to reach its final value across the entire simulation box. After 1 ns, the dissociation was completed and the temperature across the slab was ~1000 K. At this temperature, the rates of combustion reactions [Eq. (3)] were on the scale of 108 s-1, and were too slow to be captured within the simulated time of 3 ns.

Figure 2: Results of a constant-volume constant-energy DPD-RX simulation of RDX decomposition caused by a central slab heated to 3000 K at the start of the simulation. Shown are instantaneous time and position profiles for: (a) particle internal temperatures; (b) particle densities; and (c) chemical species compositions within the particles.

8

ACS Paragon Plus Environment

Page 9 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

Acknowledgements The authors acknowledge various scientists for insightful and stimulating discussions, namely, Jean-Bernard Maillet (Commissariat à l'énergie atomique et aux energies alternatives, DAM), Gabriel Stoltz (CERMICS-Ecole des Ponts Paris Tech), Timothy W. Sirk, William D. Mattson, Michael S. Sellers, and Betsy M. Rice (U.S. Army Research Laboratory). The authors also acknowledge financial support from various sources, including the Office of Naval Research, the Department of Defense High Performance Computing Modernization Program Software Application Institute for Multi-scale Reactive Modeling of Insensitive Munitions, ARL under Cooperative Agreement: W911NF-10-2-0039, and the Czech Science Foundation under project: 13-02938S. Computer time was partially provided by the METACentrum computing facility. Supporting Information Available: Descriptions concerning the implementation of the DPD-E and DPD-RX methods, the development of the MS-CG model for RDX and the CG-EOS, and various simulation details are available. This material is accessible free of charge via the Internet at http://pubs.acs.org.

9

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

References (1)

(2) (3) (4)

(5) (6) (7) (8)

(9)

(10)

(11)

(12) (13)

(14) (15)

(16) (17) (18)

Brennan, J. K.; Lísal, M. CECAM Workshop: Dissipative Particle Dynamics: Addressing Deficiencies and Establishing New Frontiers (16–18 July 2008, Lausanne, Switzerland). Molec. Simul. 2009, 35, 766-769. Venturoli, M.; Sperotto, M. M.; Kranenburg, M.; Smit., B. Mesoscopic Models of Biological Membranes. Phys. Rep. 2006, 437, 1-54. Izvekov, S.; Voth, G. A. A Multi-Scale Coarse-Graining Method for Biomolecular Systems. J. Phys. Chem. B 2005, 109, 2469-2473. Lísal, M.; Brennan, J. K.; Bonet Avalos, J. Dissipative Particle Dynamics at Isothermal, Isobaric, Isoenergetic, and Isoenthalpic Conditions using Shardlow-Like Splitting Algorithms. J. Chem. Phys. 2011, 135, 204105. Bonet Avalos, J.; Mackie, A. D. Dissipative Particle Dynamics with Energy Conservation. Europhys. Lett. 1997, 40, 141-146. Español, P. Dissipative Particle Dynamics with Energy Conservation. Europhys. Lett. 1997, 40, 631-636. Mackie, A. D.; Bonet Avalos, J.; Navas, V. Dissipative Particle Dynamics with Energy Conservation: Modelling of Heat Flow. Phys. Chem. Chem. Phys. 1999, 1, 2039-2049. Larentzos, J. P.; Brennan, J. K.; Moore, J. D.; Lísal, M.; Mattson, W. D. Parallel Implementation of Isothermal and Isoenergetic Dissipative Particle Dynamics Using Shardlow-Like Splitting Algorithms. Comp. Phys. Comm. 2014, 185, 1987. Izvekov, S.; Parrinello, M.; Burnham, C. J.; Voth, G. A. Effective Force Fields for Condensed Phase Systems from Ab Initio Molecular Dynamics Simulation: A New Method for Force-Matching. J. Chem. Phys. 2004, 120, 10896-10913. Noid, W. G.; Chu, J. W.; Ayton, G. S.; Krishna, V.; Izvekov, S.; Voth, G. A.; Das A.; Andersen, H. C. The Multiscale Coarse-Graining Method. I. A Rigorous Bridge between Atomistic and Coarse-Grained Models. J. Chem. Phys. 2008, 128, 244114. Izvekov, S.; Chung, P. W.; Rice, B. M. Particle-Based Multi-Scale Coarse Graining with Density-Dependent Potentials: Application to Molecular Crystals (Hexahydro-1,3,5trinitro-s-triazine). J. Chem. Phys. 2011, 135, 044112. Smith, G. D.; Bharadwaj, R. K. Quantum Chemistry Based Force-Field for Simulations of HMX. J. Phys. Chem. B 1999, 103, 3570-3575. Bedrov, D.; Ayyagari, C.; Smith, G. D.; Sewell, T. D.; Menikoff, R.; Zaug, J. M. Molecular Dynamics Simulations of HMX Crystal Polymorphs using a Flexible Molecule Force Field. J. Computer-Aided Materials Design, 2001, 8, 77-85. Pagonabarraga, I.; Frenkel, D. Dissipative Particle Dynamics for Interacting Systems. J. Chem. Phys. 2001, 115, 5015-5026. Tofimov, S. Y.; Nies, E. L. F.; Michels, M. A. J. Thermodynamic Consistency in Dissipative Particle Dynamics Simulations of Strongly Nonideal Liquids and Liquid Mixtures. J. Chem. Phys. 2002, 117, 9383-9394. Warren, P. B. Vapor-Liquid Coexistence in Many-Body Dissipative Particle Dynamics. Phys. Rev. E 2003, 68, 066702. Merabia, S.; Pagonabarraga, I. Density Dependent Potentials: Structure and Thermodynamics. J. Chem. Phys. 2007, 127, 054903. Merabia, S.; Bonet Avalos, J. Dewetting of a Stratified Two-Component Liquid Film on a Solid Substrate. Phys. Rev. Lett. 2008, 101, 208304. 10

ACS Paragon Plus Environment

Page 10 of 12

Page 11 of 12

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry Letters

(19) Ghoufi, A.; Malfreyt, P. Calculation of the Surface Tension from Multibody Dissipative Particle Dynamics and Monte Carlo Methods. Phys. Rev. E 2010, 82, 016706. (20) Ghoufi, A.; Emile, J.; Malfreyt, P. Recent Advances in Many Body Dissipative Particles Dynamics Simulations of Liquid-Vapor Interfaces. Eur. Phys. J. E 2013, 36, 10. (21) Warren, P. B. No-Go Theorem in Many-Body Dissipative Particle Dynamics. Phys. Rev. E 2013, 87, 045303. (22) Seaton, M. A.; Anderson, R. L.; Metz, S.; Smith, W. DL_MESO: Highly Scalable Mesoscale Simulations. Mol. Sim., 2013, 39, 796-821. (23) Nielsen, S. O.; Lopez, C. F.; Srinivas, G.; Klein, M. L. Coarse Grain Models and the Computer Simulation of Soft Materials. J. Phys.: Condens. Matter 2004, 16, R481–R512. (24) Depa, P.; Chen, C. X.; Maranas, J. K. Why are Coarse-Grained Force Fields too Fast? A Look at Dynamics of Four Coarse-Grained Polymers. J. Chem. Phys. 2001, 134, 014903. (25) Izvekov, S.; Rice, B. M. Multi-Scale Coarse-graining of Non-conservative Interactions in Molecular Liquids. J. Chem. Phys. 2014, 140, 104104. (26) Groot, R. D.; Stoyanov, S. D. Mesoscopic Model for Colloidal Particles, Powders, and Granular Solids. Phys. Rev. E 2008, 78, 051403. (27) Junghans, C.; Praprotnik, M.; Kremer, K. Transport Properties Controlled by a Thermostat: An Extended Dissipative Particle Dynamics Thermostat. Soft Matter 2008, 4, 156-161. (28) Shardlow, T. Splitting for Dissipative Particle Dynamics. SIAM J. Sci. Comput. 2003, 24, 1267-1282. (29) Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comp. Phys. 1995, 117, 1-19. (30) Hugoniot properties of a material represent the steady-state thermodynamic path a material undergoes from an initial un-shocked state to a final shocked state. The Hugoniot properties presented in this work have been calculated using the Erpenbeck method, Erpenbeck, J. J. Molecular Dynamics of Detonation I. Equation of State and Hugoniot Curve for a Simple Reactive Fluid. Phys. Rev. A 1992, 46, 6406. See SI for further details. (31) Munday, L. B.; Chung, P. W.; Rice, B. M.; Solares, S. D. Simulations of High-Pressure Phases in RDX. J. Phys. Chem. B., 2011, 115, 4378-4836. (32) Choi, C. S.; Prince, E. The Crystal Structure of Cyclotrimethylene-trinitramine. Acta Cryst. 1972, B28, 2857-2862. (33) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. A Second-Generation Reactive Empirical Bond Order (REBO) Potential Energy Expression for Hydrocarbons. J. Phys.: Condens. Matter 2002, 14, 783-802. (34) Strachan, A.; Kober, E.; van Duin, A. C. T.; Oxgaard, J.; Goddard III, W. A. Thermal Decomposition of RDX from Reactive Molecular Dynamics. J. Chem. Phys. 2005, 122, 054502. (35) Maillet, J. B.; Soulard, L.; Stoltz, G. A Reduced Model for Shock and Detonation Waves. II. The Reactive Case. Europhys. Lett. 2007, 78, 68001. (36) Maillet, J. B.; Bourasseau, E.; Desbiens, N.; Vallverdu, G.; Stoltz, G. Mesoscopic Simulations of Shock-to-Detonation Transition in Reactive Liquid High Explosive. Europhys. Lett. 2011, 96, 68007. (37) Fogler, H. S. Elements of Chemical Reaction Engineering; Prentice Hall: Englewood Cliffs, NJ; 1992. (38) Prod’homme, R. Flows of Reactive Fluids; Springer: New York, NY; 2010. 11

ACS Paragon Plus Environment

The Journal of Physical Chemistry Letters

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 12

(39) Smith, W. R.; Lísal, M. Direct Monte Carlo Simulation Methods for Nonreacting and Reacting Systems at Fixed Total Internal Energy or Enthalpy. Phys. Rev. E 2002, 66, 01114. (40) Kumbhakarna, N.; Thynell, S. T.; Chowdhury, A.; Lin, P.; Analysis of RDX-TAGzT Pseudo-Propellant Combustion with Detailed Chemical Kinetics. Combust. Theory and Mod. 2011, 15, 933-956. (41) Yetter, R. A; Dryer, F. L.; Allen, M. T.; Gatto, J. L. Development of Gas-Phase Reaction Mechanisms for Nitramine Combustion. J. Propul. Power 1995, 11, 683-697. (42) Leland, T. W.; Rowlinson, J. S.; Sather, G. A. Statistical Thermodynamics of Mixtures of Molecules of Different Sizes. Trans. Faraday Soc. 1968, 64, 1447-1460. (43) Massey, H. S. W.; Buckingham, R. A. The Low-Temperature Properties of Gaseous Helium. Proc. R. Soc. London, Ser. A 1938, 168, 378-389. (44) Ree; F. H. Simple Mixing Rule for Mixtures with exp-6 Interactions. J. Chem. Phys. 1983, 78, 409-415.

12

ACS Paragon Plus Environment