We present a hybrid semiempirical density functional tight-binding (DFTB) model with a machine learning neural network potential as a correction to the repulsive term. This hybrid model, termed machine learning tight-binding (MLTB), employs the standard self-consistent charge (SCC) DFTB formalism as a baseline, enhanced by the HIP-NN potential as an effective many-body correction for short-range pairwise repulsive interactions. The MLTB model demonstrates significantly improved transferability and extensibility compared to the SCC-DFTB and HIP-NN models. This work provides a practical computational framework for developing reliable SCC-DFTB models with additional many-body corrections that more closely approach the DFT level of accuracy. We illustrate this method with the development of an accurate model for the thorium-oxygen system, applied to the study of its nanocluster structures (ThO2)n.

A complete understanding of enzyme mechanisms requires atomistic details of chemical reactions. Quantum-based molecular dynamics simulations (QMD) are a potential source of this information, but tradeoffs between accuracy and computational cost have limited their use. Here, we develop a reactive QMD approach to investigate mechanisms of isocyanide hydratase (ICH) catalysis. In QMD simulations, molecular analogs of ICH active site residues reacted with para-nitrophenol isocyanide, forming a thioimidate intermediate. Analysis of simulated atomic configurational and charge dynamics revealed a pathway where protonation of the isocyanide carbon occurs prior to thioimidate formation. X-ray crystallography and functional assays of ICH mutants suggest this order of events might occur during enzyme catalysis. Mobile protons play essential roles in many enzymes, yet they are difficult to observe experimentally, making the ordering of proton-dependent steps ambiguous in many enzyme mechanisms. The ability to directly simulate reactions relevant to enzyme catalysis involving mobile protons demonstrates the significance of our reactive QMD approach and motivates further biological applications.

A complete understanding of enzyme mechanisms requires atomistic details of chemical reactions. Quantum-based molecular dynamics simulations (QMD) are a potential source of this information, but trade-offs between accuracy and computational cost have limited their use. We previously developed extended Lagrangian Born-Oppenheimer molecular dynamics (XL-BOMD) methods that leverage a negligible compromise in accuracy to substantially decrease the cost of QMD simulations. Here, we develop a reactive QMD approach using the latest XL-BOMD formulation, which enables efficient simulations of highly reactive systems, and use it to investigate mechanisms of intermediate formation in isocyanide hydratase (ICH) catalysis. In QMD simulations, molecular analogs of ICH active site residues reacted with para-nitrophenyl isocyanide, forming a thioimidate. Analysis of simulated atomic configurational and charge dynamics revealed a pathway where protonation of the isocyanide carbon occurs prior to thioimidate formation and suggested a possible role of Asp17 as a proton donor in the early phase of ICH catalysis. To test whether the pathway seen using the reactive QMD approach might be relevant to ICH catalysis, we performed X-ray crystallography and pre-steady-state enzyme kinetics studies of wild-type and D17N mutant ICH. Both the structure and kinetics are sensitive to the D17N mutation in a manner that is consistent with the order of the reaction steps seen in the simulations. Mobile protons play essential roles in many enzymes, yet they are difficult to observe experimentally, making the ordering of proton-dependent steps ambiguous in many enzyme mechanisms. The ability to directly simulate model reactions for the design of experiments that provide information about enzyme mechanisms involving mobile protons demonstrates the significance of our reactive QMD approach and motivates further biological applications.

Nanoindentation can be leveraged to aid in the high fidelity modeling of dislocation mediated plasticity in pentaerythritol tetranitrate (PETN), an anisotropic energetic molecular crystal. Moreover, nanoindentation tip parameters such as tip geometry, size, and degree of acuity can be utilized to target anisotropic behavior. In this work, nanoindentation was conducted across a range of orientations on the (110) face of PETN to characterize resultant yield behavior, mechanical property measurements, and resultant slip behavior and fracture initiation. Three different indentation tips were utilized: a 3-sided pyramidal Berkovich tip, a 4-sided high aspect ratio Knoop tip, and a 90° conical tip. Ultimately, indenter tip radius was documented to impact yield behavior, whereas tip geometry affected larger scale processes such as slip, and tip acuity was the dominating factor that led to fracture. The axisymmetric conical tip, serving as a baseline, showed the least amount of variation in mechanical property measurements but also the largest distribution of maximum shear stress at which initial yielding occurred. Its high degree of acuity, however, was more prone to induce fracture at higher loads. The Knoop tip was shown to be suitable for average measurements, but also for elucidation of certain anisotropic features. A distinctly higher perceived hardness at 45° was measured with the Knoop tip, indicating less dislocation motion in that direction also observed in this work via scanning probe microscopy. Lastly, the commonly used Berkovich tip was a good compromise whereby it provided a representative volume element describing the average behavior of the material. These results can be utilized to target desired anisotropic behavior in a wider range of molecular crystals, as well as to inform theoretical considerations for dislocation mediated plasticity in PETN.

The ability to predict the handling sensitivity of new organic energetic materials has been a longstanding goal. We report the synthesis and characterization of six new nitropicramide energetic materials with mixed functional groups that mimic known explosives such as nitroglycerin, erythritol tetranitrate (ETN), and pentaerythritol tetranitrate (PETN). The molecules have been studied theoretically using quantum molecular dynamics (QMD) simulations and density functional theory (DFT) calculations to identify the weakest bond in the reactants - the trigger-linkages - which control handling sensitivity, and to quantify their specific enthalpies of explosion. In good accord with the drop weight impact sensitivity data, our calculations predict that the sensitivities of the molecules are very similar owing to the small variations of the energy output and rates of trigger linkage rupture. In addition, both the QMD and DFT calculations point to the nitropicramide N-NO2 bonds as the trigger linkages rather than the more typical O-NO2 bonds. We propose that the switch of the trigger linkage from the nitrate esters to the nitramine groups arises from the strongly electron withdrawing character of the adjacent trinitrobenzene groups.

We present a new integrated experimental and modeling effort that assesses the intrinsic sensitivity of energetic materials based on their reaction rates. The High Explosive Initiation Time (HEIT) experiment has been developed to provide a rapid assessment of the high-temperature reaction kinetics for the chemical decomposition of explosive materials. This effort is supported theoretically by quantum molecular dynamics (QMD) simulations that depict how different explosives can have vastly different adiabatic induction times at the same temperature. In this work, the ranking of explosive initiation properties between the HEIT experiment and QMD simulations is identical for six different energetic materials, even though they contain a variety of functional groups. We have also determined that the Arrhenius kinetics obtained by QMD simulations for homogeneous explosions connect remarkably well with those obtained from much longer duration one-dimensional time-to-explosion (ODTX) measurements. Kinetic Monte Carlo simulations have been developed to model the coupled heat transport and chemistry of the HEIT experiment to explicitly connect the experimental results with the Arrhenius rates for homogeneous explosions. These results confirm that ignition in the HEIT experiment is heterogeneous, where reactions start at the needle wall and propagate inward at a rate controlled by the thermal diffusivity and energy release. Overall, this work provides the first cohesive experimental and first-principles modeling effort to assess reaction kinetics of explosive chemical decomposition in the subshock regime and will be useful in predictive models needed for safety assessments.

The rate of discovery of new explosives with superior energy density and performance has largely stalled. Rapid property prediction through machine learning has the potential to accelerate the discovery of new molecules by screening of large numbers of molecules before they are ever synthesized. To support this goal, we assembled a 21,000-molecule database of experimentally synthesized molecules containing energetic functional groups. Using a combination of experimental density measurements and high throughput electronic structure and atomistic calculations, we calculated detonation velocities and pressures for all 21,000 compounds. Using these values, we trained machine learning models for the prediction of density, detonation velocity and detonation pressure. Notably, our model for crystal density surpassed the accuracy of all current models and decreased the root-mean square error (RMSE) of the previous best model by 20%. This improvement in model performance relative to past works is attributed to our handling of chiral-specified Simplified Molecular-Input Line-Entry System (SMILES) strings and introduction of a new molecular descriptor, MolDensity. To elucidate descriptor importance, we evaluated interpretable descriptors in terms of importance and compared the accuracy of a statistics-driven machine learning model against a model comprised of descriptors typically assumed to control material density. The inexpensive, yet highly accurate predictions from our models should enable creation of future artificial intelligence (AI) models that are able to screen large numbers (>106) of compounds to find the highest performing compounds in terms of crystal density, detonation velocity and detonation pressure.

In Eulerian finite element simulations, the mesh moves relative to the material. After every change of position between the mesh and the material, the state variables are interpolated to the new mesh position, which is referred to as advection. Large strain crystal plasticity models are based on the multiplicative decomposition of the total deformation gradient. The stress is evaluated as a function of the thermoelastic strain, temperature, and other state variables. Advection of tensor quantities, such as the strain, is coupled with possibly significant advection errors. In an effort to reduce the advection errors, we develop two rate forms of an established dislocation density-based continuum model. To that end, we replace the multiplicative decomposition of the deformation gradient with the additive decomposition of the velocity gradient, and define the stress rate instead of the total stress. The Eulerian implementation is compared with Lagrangian calculations, and two numerical examples with severe deformation levels are presented.

The high temperature time-to-explosion experiment developed originally by Wenograd in the late 1950s has been revisited using modern standards and diagnostics. Our new High Explosives Initiation Time (HEIT) experiment heats small, 5 mg, samples of energetic materials contained within thin stainless steel hollow needles to temperatures up to 1673 K, the melting temperature of the steel, in a few 10s of microseconds using a 250 J desktop pulsed power system. Unlike Wenograd's original experiments, we detect the onset of thermal explosion via a high speed camera and optical fibers rather than electrically via the rupture of the needle. The test allows us to reliably obtain time-to-explosion as a function of input temperature, from which effective single step Arrhenius kinetics can be derived at the temperatures relevant to both common accident scenarios and shock initiation. Here, we describe in detail the calibration of the needle temperature as a function of input current and provide preliminary Arrhenius kinetics for the thermal explosion of RDX in the temperature interval 1500 > T > 850 K.

We have used molecular dynamics simulations to determine the transport properties of liquid pentaerythritol tetranitrate (PETN), an important energetic material. The density, rho, self-diffusion coefficient, D, thermal conductivity, kappa, and shear viscosity, mu, have been computed over pressures and temperatures relevant to the subshock regime (up to 1000 K and a few GPa), where PETN is known to melt prior to initiation. We find that the thermal conductivity kappa(P, T) can be represented by a simple analytical function that fits the data points with very good accuracy, even beyond the subshock regime, up to 2000 K and 20 GPa. The self-diffusion coefficient, D, exhibits nonmonotonic behavior, with notably the temperature-independent prefactor decreasing by several orders of magnitude between 0 and 2 GPa before remaining nearly constant after, and the activation energy varying little in the subshock regime before increasing linearly beyond. Lastly, the viscosity, mu, is well described by Nahme's law, which is fitted to the MD results and allows us to predict mu(P, T) for temperatures and pressures corresponding to the subshock regime. These results can be used to model the response of PETN to low-velocity impacts, where the material melts prior to the first reactions, and thermal conduction and viscosity play a crucial role.

Density functional tight binding (DFTB) models for f-element species are challenging to parametrize owing to the large number of adjustable parameters. The explicit optimization of the terms entering the semiempirical DFTB Hamiltonian related to f orbitals is crucial to generating a reliable parametrization for f-block elements, because they play import roles in bonding interactions. However, since the number of parameters grows quadratically with the number of orbitals, the computational cost for parameter optimization is much more expensive for the f-elements than for the main group elements. In this work we present a set of efficient approaches for mitigating the hurdle imposed by the large size of the parameter space. A novel group-by-orbital correction functions for two-center bond integrals was developed. With this approach the number of parameters is reduced, and it grows linearly with the number of elements, maintaining the accuracy and the number of parameters, in the case of f elements, by more than 40%. The parameter optimization step was accelerated by means of the mini-batch BFGS method. This method allows parameter optimizations with much larger training sets than other single batch methods. A stochastic optimizer was employed that helped overcome shallow local minima in the objective function. The proposed algorithm was used to parametrize the DFTB Hamiltonian for the Th-O system, which was subsequently applied to the study of ThO2 nanoparticles. The training set consisted of 6322 unique structures, which is barely feasible with conventional optimization methods. The optimized parameter set, LANL-ThO, displays good agreement with DFT-calculated properties such as energies, forces, and structures for both clusters and bulk ThO2. Benefiting from the fewer number of parameters and lower computational costs for objective function evaluations, this new approach shows its potential applications in DFTB parametrization for elements with high angular momentum, which present a challenge to conventional methods.

Predicting and controlling the failure of brittle materials against impacts have important applications in defense, mining, and medicine. To that end, the key is understanding at the mesoscale the events of crack initiation, propagation, branching, multiple crack interactions, and coalescence. Therefore, we developed an x-ray phase contrast imaging-based technique to directly visualize and quantify the cracking process. We chose to test our technique on single-crystal quartz because it has well-defined material and mechanical properties which computational models can use to accurately simulate the cracking process. Also, quartz serves as an ideal model material to developing experimental techniques/analysis and high-fidelity damage models for energetic materials and heterogenous geomaterials. To achieve the micron and nanosecond resolution required to resolve and track cracks in real-time, we use the high brilliance, spatially coherent synchrotron source at the Dynamic Compression Sector (Advanced Photon Source, Argonne National Laboratory) and the 8-frame LANL/DCS detector system coupled to a 150-μm thick single crystal LYSO scintillator. Quartz samples are uniaxially compressed at 103–104 s−1 strain rates with a custom-built Kolsky bar and stress-strain histories are measured using PDV probes. To characterize the evolving crack morphology, a physics-based inverse model is developed that converts the phase contrast-enhanced image intensity of the cracks into crack volume orientation distributions inside the sample. Using this model, we study how sample surface finish affects the dynamic behavior of cracks.

The drop weight impact experiment is used routinely to provide an initial screening of the handling sensitivity of new explosives. Despite the ubiquity and simplicity of the drop weight impact test, the physical mechanisms responsible for generating temperatures sufficiently high for the initiation of explosive reactions are not well understood. Preliminary work has shown that temperatures sufficient for ignition are in theory achievable in steady-state Poiseuille flow at the flow velocities observed experimentally with temperature-dependent shear viscosities. However, it is far from certain that such steady-state flow conditions can be achieved in the drop weight impact test because of its short duration. Therefore, here we study heat generation both experimentally and numerically to quantify the temperature distributions in molten pentaerythritol tetranitrate (PETN). Drop weight experiments have been performed starting with PETN that has been melted on a transparent anvil. The experimental results are consistent with our finite element calculations that indicate that the temperature does not approach the level needed for ignition on the time scales of the experiment. Interestingly, ignition regularly occurs in solid samples under the same conditions, so we discuss the phenomenology of ignition of both solid and molten PETN in the drop weight test and provide constraints to numerical simulations.

The catastrophic failure response of brittle solids is governed by the mechanics of crack nucleation and growth, which have been observed to be rate- and orientation-dependent. One property that is characteristic to this process is the failure strength. At low to intermediate strain rates, the failure strength has been observed to be nearly constant and equal to the strength observed under quasi-static conditions; however, at high-enough strain rates, the failure strength has been observed to become rate-dependent. The main objective of the present work is to interrogate the effects of loading rate and orientation on the failure strength of uniaxially compressed α-quartz at very high strain rates to ascertain the transition into rate sensitivity. For doing this, a miniature Kolsky bar is used to perform dynamic compression experiments on α-quartz at strain rates in the order of 103–104/s, and X-ray phase contrast imaging (XPCI) is used to directly visualize and quantify the cracking process. Experiments are carried out on nominally 1 mm and 2 mm rectangular α-quartz specimens compressed on the {-1,-1,2,0} and {-2,2,0,3} family of planes, resulting in strain rates of approximately 2000 - 5000/s to 20,000/s at the time of failure. The results show no appreciable orientation effects, suggesting that the loading configuration rather than the crystal orientation relative to the loading direction controls the orientation of crack propagation. However, the stress history and XPCI reveal that the failure strength is appreciably rate-sensitive within the present loading rate regimes. The stress history for both configurations exhibits an increasing average failure strength from around 2 GPa to 3 GPa as loading rates increase from 2000 - 5000/s to 20,000/s. The stress history and XPCI data are expected to provide crucial insight into the rate-dependence of the damage mechanisms occurring in this material.

There are few techniques available for chemists to obtain time-to-explosion data with known temperature inputs at the early stages of the design and synthesis of new explosives. In the 1960s, a technique was developed to rapidly heat milligram-quantities of confined explosives to ∼1000 K on microsecond timescales. Wenograd [Trans. Faraday Soc. 57, 1612 (1961)] loaded explosives inside stainless steel hypodermic needles, connected them to a fireset and rapidly discharged a capacitor through the steel. He obtained the temperature by measuring the needle resistance in a Wheatstone bridge arrangement and the time to explosion from a needle rupture. However, owing to the narrow-gauge needles used in the original research, the experiment was only possible with melt-castable explosives; it was never replicated, and modern diagnostics are now available with advances beyond the 1960s. Here, we report the development of the High Explosives Initiation Time (HEIT) test, which utilizes a 250 J pulsed power system to heat the needles. This work extends the Wenograd approach by using optical diagnostics, computational modeling, and advanced techniques to measure needle resistance and needle rupture. Preliminary rate information for pentaerythritol tetranitrate (PETN) will be presented.

The nonpolarizable force field for alkyl nitrates developed by Borodin et al. [J. Phys. Chem. B, 2008, 112, 734-742] has been employed to calculate selected properties of crystalline and liquid erythritol tetranitrate (ETN). The set of partial charges proposed by Borodin for pentaerythritol tetranitrate (PETN) was used except for a small correction to the H atom charges to ensure charge neutrality owing to the absence of the neopentyl carbon in ETN. The force field was used to compute the isothermal compression curve, lattice parameters, heat capacity, thermal expansivity, single crystal elastic constants, and Gruneisen parameters of crystalline ETN. The density- and temperature-dependent viscosities of liquid ETN are also reported. We anticipate that these data will be of some utility to the development of equations of state and thermomechanical models for ETN.

A thermodynamically complete equation of state (EOS) and reactive flow model of nitromethane is developed. Empirical Davis reactant and products EOSs are based on experimental data as well as higher-fidelity models from molecular dynamics simulations and thermochemical codes. Although the Arrhenius-Wescott-Stewart-Davis (AWSD) reactive flow model was originally devised for heterogeneous plastic bonded explosives, it can capture the salient features of homogeneous shock initiation, including superdetonation prior to turnover, observed experimentally in shock-to-detonation transition experiments. This is due to the model rate form having Arrhenius sensitivity to local temperature with a suitable choice of parameters. A complete set of AWSD rate parameters are calibrated to embedded electromagnetic gauge experiments. The resulting shock-to-detonation characteristics, caveats about utilizing one-dimensional simulations are presented and effects of embedded gauges are discussed.

Tailoring the molecular properties that govern energetic material sensitivity is essential to improve safety and help develop new energetic materials. Despite this need, understanding the complex chemistry and physics of explosive initiation and propagation is still a challenge. Recent work by our group has reinforced the view that explosive sensitivity under sub-shock conditions is connected to the strength of the weakest covalent bond in the molecule, that is, its "trigger linkage." These correlations have been observed with different classes of energetic molecules and indicate that "trigger linkage" bond breaking, and heat of explosion are good indicators for the sensitivity trends. Herein we report the synthesis of aliphatic energetic materials with ethane, propane and neopentane backbones. Experimental and computational studies show that the trigger linkage model, based on results from quantum molecular dynamics simulations, correctly predicts trends observed in the impact sensitivity of the molecules. However, while the model predicts the impact sensitivities of the ethane series, the neopentane series has higher impact sensitivities than predicted, which is presumably influenced by crystal packing effects.

Tabular and Debye equations of state (EoS) have been developed for single crystal 3,3 '-diamino-4,4 '-azoxyfurazan (DAAF) from density functional theory calculations of the static lattice energy (i.e., the cold curve) and the dependence of the frequencies of the vibrational normal modes on specific volume. The cold curves are represented by the three parameter MACAW reference curve developed by Lozano and Aslam. The thermal ionic contributions to the Helmholtz free energy are either interpolated and tabulated in a SESAME table or represented by a sum of five Debye models with volume-dependent Debye temperatures. The agreement between the two EoS and with the available experiments is generally good. We provide predictions for the Hugoniot locus up to 40 GPa and the temperature dependence of the heat capacity up to 2000 K. Tensors of second and third order elastic constants calculated using density functional theory, which enable the calculation of the thermoelastic free energy of DAAF under large lattice strains that probe the anharmonicity of the interatomic forces, are also reported.