^{1,a)}, Peter W. Chung

^{1}and Betsy M. Rice

^{1}

### Abstract

We describe the development of isotropic particle-based coarse-grain models for crystalline hexahydro-1,3,5-trinitro-s-triazine (RDX). The coarse graining employs the recently proposed multiscale coarse-graining (MS-CG) method, which is a particle-based force-matching approach for deriving free-energy effective interaction potentials. Though one-site and four-site coarse-grain (CG) models were parameterized from atomistic simulations of non-ordered (molten and ambient temperature amorphous) systems, the focus of the paper is a detailed study of the one-site model with a brief recourse to the four-site model. To improve the ability of the one-site model to be applied to crystalline phases at various pressures, it was found necessary to include explicit dependence on a particle density, and a new theory of local density-dependent MS-CG potentials is subsequently presented. The density-dependency is implemented through interpolation of MS-CG force fields derived at a preselected set of reference densities. The computationally economical procedure for obtaining the reference force fields starting from the interaction at ambient density is also described. The one-site MS-CG model adequately describes the atomistic lattice structure of *α*-RDX at ambient and high pressures, elastic and vibrational properties, pressure-volume curve up to *P* = 10 GPa, and the melting temperature. In the molten state, the model reproduces the correct pair structure at different pressures as well as higher order correlations. The potential of the MS-CG model is further evaluated in simulations of shocked crystalline RDX

The authors wish to thank Dr. John Brennan and Ms. Sarah Hamdan for helpful comments. This research was supported by the DoD High Performance Computing Modernization Program Software Application Institute for Multiscale Reactive Modeling of Insensitive Munitions. Computing support was provided by the DoD Supercomputer Resource Center.

I. INTRODUCTION

II. METHOD

A. MS-CG method

B. Density-dependent MS-CG potentials

C. Global density dependency

D. Local density dependency

E. Parameterization of density dependency

III. RESULTS AND DISCUSSION

A. Reference atomistic simulations and details of coarse graining

B. MS-CG potentials without density dependency

C. Density dependent one-site MS-CG model

D. Crystal structure

E. Pairwise structure in molten phase

F. Many-body effects

G. Melting

H. Elastic properties

I. Thermal properties

J. Vibrational spectrum

K. Shocked RDX

IV. CONCLUSIONS

### Key Topics

- Crystal structure
- 26.0
- Free energy
- 16.0
- Elasticity
- 12.0
- Molecular dynamics
- 11.0
- Thermodynamic properties
- 9.0

## Figures

Coarse graining of the RDX molecule into one-site and four-site representations.

Coarse graining of the RDX molecule into one-site and four-site representations.

Panel (a): One-site MS-CG force profiles as obtained from force-matching to 550 K (thin solid) and 350/550 K (thick solid) samplings. For 350/550 K force (MS-CG/0 model), possible choices of core interactions corresponding to core distances outlined in panel (b) in vertical lines are shown in dash-dash-dotted (hard-wall core) and dashed (exponential core). Panel (b): One-site MS-CG potentials at different pressures and *T* = 300 K. Vertical dashed lines mark possible selections of the core distance. Shown are potentials at *P* (GPa): 0 (filled circles), 0.5 (filled squares), 1 (filled triangles), 2 (empty circles), 3 (empty squares), 5 (empty triangles), 10 (stars). The *P* = 0 GPa potential corresponds to 350 K/550 K force in panel (a). Dotted line shows the potential corresponding to 550 K force in panel (a). On both panels dot-dashed line outlines the atomistic c.m.-c.m. RDF in arbitrary units.

Panel (a): One-site MS-CG force profiles as obtained from force-matching to 550 K (thin solid) and 350/550 K (thick solid) samplings. For 350/550 K force (MS-CG/0 model), possible choices of core interactions corresponding to core distances outlined in panel (b) in vertical lines are shown in dash-dash-dotted (hard-wall core) and dashed (exponential core). Panel (b): One-site MS-CG potentials at different pressures and *T* = 300 K. Vertical dashed lines mark possible selections of the core distance. Shown are potentials at *P* (GPa): 0 (filled circles), 0.5 (filled squares), 1 (filled triangles), 2 (empty circles), 3 (empty squares), 5 (empty triangles), 10 (stars). The *P* = 0 GPa potential corresponds to 350 K/550 K force in panel (a). Dotted line shows the potential corresponding to 550 K force in panel (a). On both panels dot-dashed line outlines the atomistic c.m.-c.m. RDF in arbitrary units.

Comparison of RDX unit cells calculated at *T* = 4.5 K using atomistic and MS-CG-D models. For the atomistic model unit cell, red balls denote the molecular c.m. locations superimposed upon the ball-and-stick representations of the molecules, whereas green balls denote the CG interaction sites.

Comparison of RDX unit cells calculated at *T* = 4.5 K using atomistic and MS-CG-D models. For the atomistic model unit cell, red balls denote the molecular c.m. locations superimposed upon the ball-and-stick representations of the molecules, whereas green balls denote the CG interaction sites.

Comparison of the *α*-RDX crystal lattice structure from simulations using the atomistic (red balls denote molecular c.m.) and MS-CG-D (green balls) models. Panel (a): [001]-view of lattice from the MS-CG-D simulation; Panel (b): Same view as in panel (a) for the atomistic model; Panel (c): [100]-view of superimposed lattices; Panel (d): [010]-view. In panels (a), (b), and (c) dotted and dashed lines outline location of planes in the atomistic lattice, which correspond to planes shown in dotted and solid lines for MS-CG-D lattice. In panel (b), atomistic planes (dashed) move to new locations in MS-CG-D lattice (solid) in panel (a). In panel (c), the two distinct planes in the atomistic representation (dashed) collapse into a single plane in the MS-CG representation (solid line).

Comparison of the *α*-RDX crystal lattice structure from simulations using the atomistic (red balls denote molecular c.m.) and MS-CG-D (green balls) models. Panel (a): [001]-view of lattice from the MS-CG-D simulation; Panel (b): Same view as in panel (a) for the atomistic model; Panel (c): [100]-view of superimposed lattices; Panel (d): [010]-view. In panels (a), (b), and (c) dotted and dashed lines outline location of planes in the atomistic lattice, which correspond to planes shown in dotted and solid lines for MS-CG-D lattice. In panel (b), atomistic planes (dashed) move to new locations in MS-CG-D lattice (solid) in panel (a). In panel (c), the two distinct planes in the atomistic representation (dashed) collapse into a single plane in the MS-CG representation (solid line).

Panel (a): c.m. RDFs from atomistic (solid), MS-CG 350/550 K (dotted), and 550 K (dashed) simulations of liquid RDX at *T* = 550 K using one-site (black) and four-site (color) MS-CG models. Panel (b): c.m. RDFs from atomistic (solid) and one-site MS-CG-D (dashed) simulations at different *P*. Shown are *P* (GPa): 0 (cyan), 0.5 (black), 1 (blue), 2 (green), and 5 (red). Panel (c): c.m. RDFs from atomistic (thick solid) and MS-CG-D (thin solid) simulations of crystalline RDX at *T* = 300 K. The corresponding running coordination number *N* _{ c } shown in dashed lines.

Panel (a): c.m. RDFs from atomistic (solid), MS-CG 350/550 K (dotted), and 550 K (dashed) simulations of liquid RDX at *T* = 550 K using one-site (black) and four-site (color) MS-CG models. Panel (b): c.m. RDFs from atomistic (solid) and one-site MS-CG-D (dashed) simulations at different *P*. Shown are *P* (GPa): 0 (cyan), 0.5 (black), 1 (blue), 2 (green), and 5 (red). Panel (c): c.m. RDFs from atomistic (thick solid) and MS-CG-D (thin solid) simulations of crystalline RDX at *T* = 300 K. The corresponding running coordination number *N* _{ c } shown in dashed lines.

Distributions of the order parameter Eqs. (38) and (39) with different coordination number *N* _{ c } for atomistic (thick) and MS-CG-D (thin) models. The MS-CG-D results shown for reference crystal structure {R_{ IK }} from atomistic (solid) and MS-CG-D (dashed) models.

Distributions of the order parameter Eqs. (38) and (39) with different coordination number *N* _{ c } for atomistic (thick) and MS-CG-D (thin) models. The MS-CG-D results shown for reference crystal structure {R_{ IK }} from atomistic (solid) and MS-CG-D (dashed) models.

Panel (a): Snapshots of the liquid-solid interface from simulation of slab melting at *T* = 490 K, at which close-packed structure formed at the interface (upper), and more common configuration (bottom). In the upper snapshot, the closed-packed domains of different from hP2 structure are contoured by dotted line. Such domains are artifact of the model reflecting limited transferability of the model to heterogeneous structures such as interfaces. These domains are similar to those shown in inset to Fig. 4(a) and discussed in the text. Panel (b): Equilibrium global order parameter, Eq. (39) with *N* _{ c } = 12 as a function of *T* for crystal. Insert shows for crystal at *T* = 450 and 550 K.

Panel (a): Snapshots of the liquid-solid interface from simulation of slab melting at *T* = 490 K, at which close-packed structure formed at the interface (upper), and more common configuration (bottom). In the upper snapshot, the closed-packed domains of different from hP2 structure are contoured by dotted line. Such domains are artifact of the model reflecting limited transferability of the model to heterogeneous structures such as interfaces. These domains are similar to those shown in inset to Fig. 4(a) and discussed in the text. Panel (b): Equilibrium global order parameter, Eq. (39) with *N* _{ c } = 12 as a function of *T* for crystal. Insert shows for crystal at *T* = 450 and 550 K.

Comparison of isotherms for crystalline RDX at *T* = 300 K (crosses) and molten RDX at *T* = 550 K (circles) from atomistic (solid) and MS-CG-D (dashed) models. Panel (a): *V*/*V* _{0} − *P* isotherm, where *V* _{0} is ambient volume; Panel (b) Δ*U* ^{ conf } − *P* isotherm, where is the change in the total internal potential energy with respect to the ambient value .

Comparison of isotherms for crystalline RDX at *T* = 300 K (crosses) and molten RDX at *T* = 550 K (circles) from atomistic (solid) and MS-CG-D (dashed) models. Panel (a): *V*/*V* _{0} − *P* isotherm, where *V* _{0} is ambient volume; Panel (b) Δ*U* ^{ conf } − *P* isotherm, where is the change in the total internal potential energy with respect to the ambient value .

Panel (a): All-atom VDOS (dashed), atomistic VDOS of c.m. (thin solid), and MS-CG-D (thick solid) in arbitrary units. The vertical dashed lines mark the parts of spectrum: (I) Only lattice modes associated with movements and rotations of molecular c.m.; (II) Mix of c.m. lattice modes with wagging and torsional modes of nitro groups. Panel (b): All-atom VDOS with portion of spectrum shown in panel (a) marked with a vertical solid line. The vertical dashed line (150 cm^{−1}) marks the low frequency region in which lattice modes are present.

Panel (a): All-atom VDOS (dashed), atomistic VDOS of c.m. (thin solid), and MS-CG-D (thick solid) in arbitrary units. The vertical dashed lines mark the parts of spectrum: (I) Only lattice modes associated with movements and rotations of molecular c.m.; (II) Mix of c.m. lattice modes with wagging and torsional modes of nitro groups. Panel (b): All-atom VDOS with portion of spectrum shown in panel (a) marked with a vertical solid line. The vertical dashed line (150 cm^{−1}) marks the low frequency region in which lattice modes are present.

Panel (a): Hugoniot curves of RDX from the atomistic model (circles), the MG-CG-D model (triangles), and experiment (solid line). Panel (b): Pressure vs. temperature for the calculated shock Hugoniot of RDX from atomistic (circles) and MS-CG-D (squares) models. Dashed line depicts the MS-CG-D melting curve using the Kraut-Kennedy relation as discussed in the text. Right *T*-axis is for the atomistic model and left *T*-axis is for the MS-CG-D model.

Panel (a): Hugoniot curves of RDX from the atomistic model (circles), the MG-CG-D model (triangles), and experiment (solid line). Panel (b): Pressure vs. temperature for the calculated shock Hugoniot of RDX from atomistic (circles) and MS-CG-D (squares) models. Dashed line depicts the MS-CG-D melting curve using the Kraut-Kennedy relation as discussed in the text. Right *T*-axis is for the atomistic model and left *T*-axis is for the MS-CG-D model.

The profiles of (top frame) molecular density ρ_{ mol }, (middle frame) effective translational temperature *T* _{ eff }, atomistic temperature *T* _{ atm } (filled circles), and (bottom frame) global order parameter [Eq. (39)] in shocked RDX from atomistic (thick) and MS-CG-D (thin) simulations. Panel (a): *V* _{ p } = 1 km/s; Panel (b): *V* _{ p } = 3 km/s. The MS-CG-D ρ_{ mol } curve is shifted up by difference Δρ_{0} = 0.215 in the values of atomistic and MS-CG-D molecular density of crystal at *T* = 4.5 K for a better comparison.

The profiles of (top frame) molecular density ρ_{ mol }, (middle frame) effective translational temperature *T* _{ eff }, atomistic temperature *T* _{ atm } (filled circles), and (bottom frame) global order parameter [Eq. (39)] in shocked RDX from atomistic (thick) and MS-CG-D (thin) simulations. Panel (a): *V* _{ p } = 1 km/s; Panel (b): *V* _{ p } = 3 km/s. The MS-CG-D ρ_{ mol } curve is shifted up by difference Δρ_{0} = 0.215 in the values of atomistic and MS-CG-D molecular density of crystal at *T* = 4.5 K for a better comparison.

## Tables

Coefficients *A* _{ n } of the least-squares fit of the one-site MS-CG force at *P* = 0 GPa using the polynomial equation (35) (second column) and Chebyshev series equations (36) and (37) (third column), with power of the term *n* shown in the first column. Atomic units for force and distance were used. For the Chebyshev series, the following radii were used in Eq. (37): *R* _{ core } = 8.5 a.u., *R* _{ cut } = 29.011 a.u. The expansions were switched to zero using the linear switching function as in Eq. (33) with δ_{ R } = 4.0 a.u. for a polynomial expansion and with δ_{ R } = 2.0 a.u. for a Chebyshev expansion. The fourth column shows the values of parameters for the density dependent contribution Eq. (33). At *R* < *R* _{ core } = 8.5 a.u., the force extrapolated as discussed in Sec. III B. A cutoff of 1.535 nm must be applied to the expansions.

Coefficients *A* _{ n } of the least-squares fit of the one-site MS-CG force at *P* = 0 GPa using the polynomial equation (35) (second column) and Chebyshev series equations (36) and (37) (third column), with power of the term *n* shown in the first column. Atomic units for force and distance were used. For the Chebyshev series, the following radii were used in Eq. (37): *R* _{ core } = 8.5 a.u., *R* _{ cut } = 29.011 a.u. The expansions were switched to zero using the linear switching function as in Eq. (33) with δ_{ R } = 4.0 a.u. for a polynomial expansion and with δ_{ R } = 2.0 a.u. for a Chebyshev expansion. The fourth column shows the values of parameters for the density dependent contribution Eq. (33). At *R* < *R* _{ core } = 8.5 a.u., the force extrapolated as discussed in Sec. III B. A cutoff of 1.535 nm must be applied to the expansions.

Lattice parameters *a*, *b*, *c* and density *ρ* of the RDX crystal at different *T* and *P* = 0 GPa from atomistic and MS-CG-D models.

Lattice parameters *a*, *b*, *c* and density *ρ* of the RDX crystal at different *T* and *P* = 0 GPa from atomistic and MS-CG-D models.

Young's (*E*) and shear (*G*) moduli at different *T* (K) and *P* (GPa). Units are GPa.

Young's (*E*) and shear (*G*) moduli at different *T* (K) and *P* (GPa). Units are GPa.

Article metrics loading...

Full text loading...

Commenting has been disabled for this content