A Generalized Crystal-Cutting Method (GCCM) is developed that automates construction of three-dimensionally periodic simulation cells containing arbitrarily oriented single crystals and thin films, two-dimensionally (2D) infinite crystal–crystal homophase and heterophase interfaces, and nanostructures with intrinsic NN-fold interfaces. The GCCM is based on a simple mathematical formalism that facilitates easy definition of constraints on cut crystal geometries. The method preserves the translational symmetry of all Bravais lattices and thus can be applied to any crystal described by such a lattice including complicated, low-symmetry molecular crystals. Implementations are presented with carefully articulated combinations of loop searches and constraints that drastically reduce computational complexity compared to simple loop searches. Orthorhombic representations of monoclinic and triclinic crystals found using the GCCM overcome some limitations in standard distributions of popular molecular dynamics software packages. Stability of grain boundaries in ββ-HMX was investigated using molecular dynamics and molecular statics simulations with 2D infinite crystal–crystal homophase interfaces created using the GCCM. The order of stabilities for the four grain boundaries studied is predicted to correlate with the relative prominence of particular crystal faces in lab-grown ββ-HMX crystals. We demonstrate how nanostructures can be constructed through simple constraints applied in the GCCM framework. Example GCCM constructions are shown that are relevant to some current problems in materials science, including shock sensitivity of explosives, layered electronic devices, and pharmaceuticals.

Nanoindentation of the insensitive energetic molecular crystal 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) was studied using constant-temperature and constant-energy molecular dynamics simulations. Displacement-controlled indentations at constant velocity were performed using a rigid, spherical indenter on the three principal crystallographic planes, (100), (010), and (001). The force–displacement curve for the (001) (basal) plane exhibits a distinct elastic region in agreement with the analytical solution for indentation of an anisotropic half-space by a parabola of revolution. Stiffening precedes inelastic deformation, and the elastic–inelastic transition occurs with kinking and delamination of the layered basal planes and significant pile-up. The predicted nanoindentation hardness on the basal plane is 1.02 ± 0.09 GPa. Nanoindentation on the (100) and (010) (nonbasal) planes yields a non-Hertzian response; this is attributed to an effective “softening” due to elastic bending of the molecular layers. Significant heating occurs in the vicinity of the indenter during indentation with a higher temperature rise for the basal plane.

Relaxation of nanoscopic idealized hot spots in the layered molecular crystalline explosive 1,3,5-triamino-2,4,6-trinitrobenzene (TATB) was studied along two crystallographically relevant directions using all-atom molecular dynamics (MD) simulations and continuum-level models based on the diffusive heat equation. Characteristics of relaxation from initial one-dimensional, nonequilibrium temperature distributions in the crystal were determined using MD simulations. Results from these MD simulations were fit to and compared with solutions for the one-dimensional diffusive heat equation by treating the thermal diffusivity as a parameter to assess the validity of using continuum models to describe heat transport in TATB on length scales below approximately 10 nm. It is found that energy transport is predominantly diffusive both within and between the crystal layers, even on the length scale of the unit cell. Continuum-level predictions are in excellent agreement with MD results for relaxation of a 2.5 nm wide hot spot along the direction normal to the crystal layers. The corresponding best-fit diffusivity is found to agree with independent predictions for the thermal conductivity along that direction to within 5%. Modest discrepancies in the continuum predictions are present for the direction nominally within the layers, and the apparent conductivity along that direction is found to be approximately 25% lower than expected. While some of these discrepancies are reconciled by treating the thermal conductivity as a function of temperature, the lower apparent conductivity indicates that some phonon modes might transfer energy ballistically within the layers on length scales below 10 nm.

Molecular dynamics (MD) simulations were used to study shock wave passage with normal incidence through the equilibrium interface between (100)-oriented nitromethane and the melt. The simulations were performed using the fully flexible, nonreactive SRT force field (Sorescu, D. C.; Rice, B. M.; Thompson, D. L. The Journal of Physical Chemistry B 2000, 104, 8406–8419). The local kinetic energies (intermolecular, intramolecular, and total) and stress states differ significantly in the liquid and crystal regions, and depend on whether the shock is initiated in the crystal or liquid. The number and spatial distributions of shock-induced molecular disorientations in the crystal for shocks initiated in the crystal are similar to those obtained for analogous simulations for a completely crystalline sample; however, substantial differences in the extent and distribution of shock-induced molecular disorientations in the crystal region were observed when the shock was initiated in the liquid. All three measures of kinetic temperature in the crystal region are higher when the shock is initiated in the crystal than when it is initiated in the liquid. Kinetic temperature profiles exhibit features in the vicinity of the interface considerably different from those in either bulk phase. The shock-induced local mechanical states (von Mises stress) indicate that the crystal is less able to support shear stresses when the shock is initiated in the crystal than when it is initiated in the liquid. There is a strong reflection back into the liquid when the shock wave passes through the liquid and encounters the interface with the crystal. This causes a large increase in the potential energy of the liquid and limits the amount of energy transmitted into the crystal, which limits the molecular disorientations in the crystal. Thus, a shock from liquid to crystal yields less inelastic deformation in the crystal.

Molecular dynamics simulations were used to study the shock-induced collapse of cylindrical pores in oriented single crystals of the energetic material α-1,3,5-trinitroperhydro-1,3,5-triazine (α-RDX). The shock propagation direction was parallel to the [100] crystal direction and the cylinder axis of the initially 35.0 nm diameter pore was parallel to [010]. Features of the collapse were studied for Rankine–Hugoniot shock pressures Ps = 9.71, 24.00, and 42.48 GPa. Pore collapse for the weak shock is dominated by visco-plastic deformation in which the pore pinches shut without jet formation and with little penetration of the upstream material into the downstream pore wall. For the strong shock the collapse is hydrodynamic-like and results in the formation of a jet that penetrates significantly into the downstream pore wall. Material flow during collapse was characterized by examining the spread and mixing of sets of initially contiguous molecules and evolution of local velocity fields. Local disorder during collapse was assessed using time autocorrelation functions for molecular rotation. Energy deposition and localization was studied using spatial maps of temperature and pressure calculated as functions of time.

Terahertz infrared absorption spectra of crystalline pentaerythritol tetranitrate were obtained using classical molecular dynamics simulations for temperature T = 298 K and hydrostatic pressures P = 0, 1, 2, and 3 GPa. Two approaches were used to calculate the spectra; one based on combined normal-mode analysis and mode-relaxation calculations and the other on Fourier analysis of the dipole–dipole time autocorrelation function. The two methods yield similar, though not identical, positions and relative spectral amplitudes for all pressures studied. The predicted spectra exhibit an overall blue shift with increasing pressure, accompanied by a decrease in the integral absorption.Graphical abstractHighlights► THz IR absorption spectra of PETN crystal were calculated using two MD methods. ► Spectra were calculated at T = 298 K for pressures P = 0, 1, 2, and 3 GPa. ► The two approaches yield similar predictions for the spectra. ► The terahertz spectra exhibit an overall blue shift with increasing pressure. ► Pressure-dependent Raman spectra at higher frequencies provide indirect validation.