Properties of Erythritol Tetranitrate from Molecular Dynamics Simulation.

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.

Halogenated PETN derivatives: interplay between physical and chemical factors in explosive sensitivity.

Determining the factors that influence and can help predict energetic material sensitivity has long been a challenge in the explosives community. Decades of literature reports identify a multitude of factors both chemical and physical that influence explosive sensitivity; however no unifying theory has been observed. Recent work by our team has demonstrated that the kinetics of "trigger linkages" (i.e., the weakest bonds in the energetic material) showed strong correlations with experimental drop hammer impact sensitivity. These correlations suggest that the simple kinetics of the first bonds to break are good indicators for the reactivity observed in simple handling sensitivity tests. Herein we report the synthesis of derivatives of the explosive pentaerythritol tetranitrate (PETN) in which one, two or three of the nitrate ester functional groups are substituted with an inert group. Experimental and computational studies show that explosive sensitivity correlates well with Q (heat of explosion), due to the change in the number of trigger linkages removed from the starting material. In addition, this correlation appears more significant than other observed chemical or physical effects imparted on the material by different inert functional groups, such as heat of formation, heat of explosion, heat capacity, oxygen balance, and the crystal structure of the material.

SCC-DFTB Parameters for Fe-C Interactions.

We present an optimized density-functional tight-binding (DFTB) parameterization for iron-based complexes based on the popular trans3d set of parameters. The transferability of the original and optimized parameterizations is assessed using a set of 50 iron complexes, which include carbonyl, cyanide, polypyridine, and cyclometalated ligands. DFTB-optimized structures predicted using the trans3d parameters show a good agreement with both experimental crystal geometries and density functional theory (DFT)-optimized structures for Fe-N bond lengths. Conversely, Fe-C bond lengths are systematically overestimated. We improve the accuracy of Fe-C interactions by truncating the Fe-O repulsive potential and reparameterizing the Fe-C repulsive potential using a training set of six isolated iron complexes. The new trans3d*-LANLFeC parameter set can produce accurate Fe-C bond lengths in both geometry optimizations and molecular dynamics (MD) simulations, without significantly affecting the accuracy of Fe-N bond lengths. Moreover, the potential energy curves of Fe-C interactions are considerably improved. This improved parameterization may open the door to accurate MD simulations at the DFTB level of theory for large systems containing iron complexes, such as sensitizer-semiconductor assemblies in dye-sensitized solar cells, that are not easily accessible with DFT approaches because of the large number of atoms.

Insight into the Chemistry of PETN Under Shock Compression Through Ultrafast Broadband Mid-Infrared Absorption Spectroscopy.

Thin films of pentaerythritol tetranitrate (PETN) were shock compressed using the laser driven shock apparatus at Los Alamos National Laboratory (LANL). Two spectroscopic probes were available to this apparatus: visible white light transient absorption spectroscopy (VIS) from 400 to 700 nm and mid-infrared transient absorption spectroscopy (MIR) from 1150 to 3800 cm-1. Important PETN vibrational modes are the symmetric and antisymmetric NO2 stretches at 1280 and 1650 cm-1, respectively, as well as CH stretches at ∼2900 cm-1. Shock strength was varied from approximately 3 to 55 GPa to span from the chemically unreactive regime to the regime in which fast chemical reaction took place on the 250 ps time scale of the measurements. VIS and MIR results suggest irreversible chemistry was induced in PETN at pressures above 30 GPa. At lower shock pressures, the spectroscopy showed minimal changes attributable to pressure induced effects. Under the higher-pressure reactive conditions, the frequency region at the antisymmetric NO2 stretch mode had a significantly increased absorption while the region around the symmetric NO2 stretch did not. No observable increased absorption occurred in the higher frequency regions where CH-, NH-, and OH- bond absorptions would be observed. A broad absorption appeared on the shoulder at the red-edge of the CO2 vibrational band around 2200 cm-1. In addition to the experiments, reactive molecular dynamics were carried out under equivalent shock conditions to correlate the evolution of the infrared spectrum to molecular processes. The simulations show results consistent to experiments up to 30 GPa but suggest that NO and NO2 related features provided the strongest contributions to the shocked infrared changes. Proposed mechanisms for shocked PETN chemistry are analyzed as consistent or inconsistent with the data presented here. Our experimental data suggests C≡O or N2O bond formation, nitrite formation, and absence of significant hydroxyl or amine concentrations in the initial chemistry steps in PETN shocked above 30 GPa.

Reaction Rates in Nitromethane under High Pressure from Density Functional Tight Binding Molecular Dynamics Simulations.

We use density functional tight binding (DFTB) molecular dynamics (MD) simulations to determine the reaction rates of nitromethane CH3NO2 (NM) under high pressure (P = 14-28 GPa), and temperature (T = 1450-1850 K). DFTB-MD simulations performed with the same initial conditions (P0, T0) reveal a stochastic behavior, both in terms of reaction times and chemical paths. By running series of MD simulations, we are able to obtain average reaction times with quantified errors and devise a simple two-step model for NM explosion: ignition/explosion. While our model bypasses the chemical complexity due to the numerous reaction paths and intermediates observed during reactions, the chemistry is accounted for via the accurate parameterization of the DFTB model, and our results suggest a single main reaction pathway for the pressure range considered here, dominated in the earlier stages by the formation of the aci-ion, CH2NOO-. By fitting our data to a Frank-Kamenetskii model, we extract prefactors and pressure-independent activation energies and volumes for the ignition and explosion stages. A two-step model is then built and compared to experimental observations. Single and two-step Arrhenius models are also provided for comparison with literature data. This work presents an efficient way of investigating the reactivity of high explosives by performing electronic structure-based MD simulations and provides reaction rates for simplified models that can be implemented into hydrocodes.

Tight-Binding Modeling of Uranium in an Aqueous Environment.

The first density functional tight-binding (DFTB) parameters for uranium, oxygen, and hydrogen chemistry are reported, which enable quantum molecular dynamics simulations that will be instrumental in understanding actinide speciation, reaction mechanisms, and kinetics. These parameters were fitted to atomization energies and forces obtained from density functional theory with a training set of small molecules that includes various oxidation states. The energetic results with these DFTB parameters for various reactions of hydration, hydrolysis, dimerization, and isomerization demonstrate that the DFTB method can qualitatively capture the correct chemistry with a small systematic deviation from the density functional theory reference values. Structural results on the molecules not in the training set, including dimers, show generally good agreement with the reference and demonstrate the transferability of these first DFTB parameters for uranium chemistry.

Development of Density Functional Tight-Binding Parameters Using Relative Energy Fitting and Particle Swarm Optimization.

We provide a strategy to optimize density functional tight-binding (DFTB) parameterization for the calculation of the structures and properties of organic molecules consisting of hydrogen, carbon, nitrogen, and oxygen. We utilize an objective function based on similarity measurements and the Particle Swarm Optimization (PSO) method to find an optimal set of parameters. This objective function considers not only the common DFTB descriptors of binding energies and atomic forces but also incorporates relative energies of isomers into the fitting procedure for more chemistry-driven results. The quality in the description of the binding energies and atomic forces is measured based on the Ballester similarity index and relative energies through a similarity index induced by the Levenshtein edit distance to quantify the correct energetic order of isomers. Training and testing datasets were created to include all relevant chemical functional groups. The accuracy of this strategy is assessed, and its range of applicability is discussed by comparison against our previous parameterization [A. Krishnapriyan, et al., J. Chem. Theory Comput. 13, 6191 (2017)]. The improved performance of the new DFTB parameterization is validated with respect to the density functional theory large datasets QM-9 [R. Ramakrishnan, et al., Sci. Data 1, 140022 (2014)] and ANI-1 [J. S. Smith, et al., Sci. Data 4, 170193 (2017)], where excellent agreement is found between the structures and properties available in these datasets, and the ones obtained with DFTB.

Ranking the Drop-Weight Impact Sensitivity of Common Explosives Using Arrhenius Chemical Rates Computed from Quantum Molecular Dynamics Simulations.

Drop-weight impact tests are used routinely to characterize the handling safety of explosives. Numerous studies have sought to connect various physical and chemical properties of the energetic molecules and materials to their measured impact sensitivities. Wenograd in the early 1960s demonstrated that there is a strong dependency of the drop-heights on the critical temperatures required for explosives to undergo prompt reactions. Reactive quantum molecular dynamics simulations with the lanl31 density functional tight binding model have been used to compute the delay time before the thermal explosion of the secondary explosives erythritol tetranitrate (ETN), pentaerythritol tetranitrate (PETN), cyclotrimethylene trinitramine (RDX), cyclotetramethylene tetranitramine (HMX), trinitrotolune (TNT), and 3,3'-diamino-4,4'-azoxyfurazan (DAAF) as a function of the initial temperature and pressure. The delay time to explosion data are consistent with Arrhenius chemical kinetics, which is expected for thermally activated processes in materials and in accord with experimental measurements. The critical temperatures required for the materials to undergo prompt explosions display the same dependence on drop height as was observed by Wenograd. Hence, quantum-based reactive molecular dynamics simulations are potentially a tool for ranking the drop-weight impact sensitivity and handling safety of explosives.

Parallel replica dynamics simulations of reactions in shock compressed liquid benzene.

The study of the long-term evolution of slow chemical reactions is challenging because quantum-based reactive molecular dynamics simulation times are typically limited to hundreds of picoseconds. Here, the extended Lagrangian Born-Oppenheimer molecular dynamics formalism is used in conjunction with parallel replica dynamics to obtain an accurate tool to describe the long-term chemical dynamics of shock-compressed benzene. Langevin dynamics has been employed at different temperatures to calculate the first reaction times in liquid benzene at pressures and temperatures consistent with its unreacted Hugoniot. Our coupled engine runs for times on the order of nanoseconds (one to two orders of magnitude longer than traditional techniques) and is capable of detecting reactions that are characterized by rates significantly slower than we could study before. At lower pressures and temperatures, we mainly observe Diels-Alder metastable reactions, whereas at higher pressures and temperatures we observe stable polymerization reactions.

Probing ultrafast shock-induced chemistry in liquids using broad-band mid-infrared absorption spectroscopy.

We probe shock-induced chemistry in two organic liquids by measuring broadband, midinfrared absorption in the 800-1400 cm-1 frequency range. To test this new method and understand the signatures of chemical reactions in time resolved vibrational spectra, we compared liquid benzene shocked to unreactive conditions (shocked to a pressure of 18 GPa for a duration of 300 ps) to nitromethane under reactive conditions (25 GPa). We see clear signatures of shock-induced chemistry that are distinguishable from the pressure- and temperature-induced changes in vibrational mode shapes. While shocked benzene shows primarily a broadening and shifting of the vibrational modes, the nitromethane vibrational modes vanish once the shock wave enters the liquid and simultaneously, a spectrally broad feature appears that we interpret as the infrared spectrum of the complex mixture of product and intermediate species. To further interpret these measurements, we compare them to reactive quantum molecular dynamics simulations, which gives qualitatively consistent results. This work demonstrates a promising method for time resolving shock-induced chemistry, illustrating that chemical reactions produce distinct changes in the vibrational spectra.

Transferable density functional tight binding for carbon, hydrogen, nitrogen, and oxygen: Application to shock compression.

A new parameterization for density functional tight binding (DFTB) theory, lanl31, has been developed for molecules containing carbon, hydrogen, nitrogen, and oxygen. Optimal values for the Hubbard Us, on-site energies, and the radial dependences of the bond integrals and repulsive potentials were determined by numerical optimization using simulated annealing to a modest database of ab initio-calculated atomization energies and interatomic forces. The transferability of the optimized DFTB parameterization has been assessed using the CHNO subset of the QM-9 database [R. Ramakrishnan et al., Sci. Data 1, 140022 (2014)]. These analyses showed that the errors in the atomization energies and interatomic forces predicted by our model are small and in the vicinity of the differences between density functional theory calculations with different basis sets and exchange-correlation functionals. Good correlations between the molecular dipole moments and HOMO-LUMO gaps predicted by lanl31 and the QM-9 data set are also found. Furthermore, the errors in the atomization energies and forces derived from lanl31 are significantly smaller than those obtained from the ReaxFF-lg reactive force field for organic materials [L. Liu et al., J. Phys. Chem. A 115, 11016 (2011)]. The lanl31 DFTB parameterization for C, H, N, and O has been applied to the molecular dynamics simulation of the principal Hugoniot of liquid nitromethane, liquid benzene, liquid nitrogen, pentaerythritol tetranitrate, trinitrotoluene, and cyclotetramethylene tetranitramine. The computed and measured Hugoniot loci are in excellent agreement with experiment, and we discuss the sensitivity of the loci to the underestimated shock heating that is a characteristic of classical molecular dynamics simulations.

Numerical Optimization of Density Functional Tight Binding Models: Application to Molecules Containing Carbon, Hydrogen, Nitrogen, and Oxygen.

New parametrizations for semiempirical density functional tight binding (DFTB) theory have been developed by the numerical optimization of adjustable parameters to minimize errors in the atomization energy and interatomic forces with respect to ab initio calculated data. Initial guesses for the radial dependences of the Slater-Koster bond integrals and overlap integrals were obtained from minimum basis density functional theory calculations. The radial dependences of the pair potentials and the bond and overlap integrals were represented by simple analytic functions. The adjustable parameters in these functions were optimized by simulated annealing and steepest descent algorithms to minimize the value of an objective function that quantifies the error between the DFTB model and ab initio calculated data. The accuracy and transferability of the resulting DFTB models for the C, H, N, and O system were assessed by comparing the predicted atomization energies and equilibrium molecular geometries of small molecules that were not included in the training data from DFTB to ab initio data. The DFTB models provide accurate predictions of the properties of hydrocarbons and more complex molecules containing C, H, N, and O.

Free Energy Based Equation of State for Pentaerythritol Tetranitrate.

An equation of state for the energetic molecular crystal pentaerythritol tetranitrate (PETN) has been developed from a parametrized model for its Helmholtz free energy. The ion motion contribution to the free energy is represented by a sum of Debye models for the vibrational modes of mainly lattice phonon and intramolecular character. The dependence of the frequencies of the normal modes on density is captured using the quasi-harmonic approximation whereby the Debye temperatures for both populations of modes depend explicitly on specific volume. The dependence of the Debye temperatures on specific volume was parametrized to normal-mode frequencies computed from solid state dispersion-corrected density functional theory. The model provides a good description of the thermophysical properties of PETN. The equation of state has been applied to the calculation of thermodynamic states along the principal Hugoniot of single crystal PETN.

Efficient parallel linear scaling construction of the density matrix for Born-Oppenheimer molecular dynamics.

We present an algorithm for the calculation of the density matrix that for insulators scales linearly with system size and parallelizes efficiently on multicore, shared memory platforms with small and controllable numerical errors. The algorithm is based on an implementation of the second-order spectral projection (SP2) algorithm [ Niklasson, A. M. N. Phys. Rev. B 2002 , 66 , 155115 ] in sparse matrix algebra with the ELLPACK-R data format. We illustrate the performance of the algorithm within self-consistent tight binding theory by total energy calculations of gas phase poly(ethylene) molecules and periodic liquid water systems containing up to 15,000 atoms on up to 16 CPU cores. We consider algorithm-specific performance aspects, such as local vs nonlocal memory access and the degree of matrix sparsity. Comparisons to sparse matrix algebra implementations using off-the-shelf libraries on multicore CPUs, graphics processing units (GPUs), and the Intel many integrated core (MIC) architecture are also presented. The accuracy and stability of the algorithm are illustrated with long duration Born-Oppenheimer molecular dynamics simulations of 1000 water molecules and a 303 atom Trp cage protein solvated by 2682 water molecules.

Extended Lagrangian formulation of charge-constrained tight-binding molecular dynamics.

The extended Lagrangian Born-Oppenheimer molecular dynamics formalism [Niklasson, Phys. Rev. Lett., 2008, 100, 123004] has been applied to a tight-binding model under the constraint of local charge neutrality to yield microcanonical trajectories with both precise, long-term energy conservation and a reduced number of self-consistent field optimizations at each time step. The extended Lagrangian molecular dynamics formalism restores time reversal symmetry in the propagation of the electronic degrees of freedom, and it enables the efficient and accurate self-consistent optimization of the chemical potential and atomwise potential energy shifts in the on-site elements of the tight-binding Hamiltonian that are required when enforcing local charge neutrality. These capabilities are illustrated with microcanonical molecular dynamics simulations of a small metallic cluster using an sd-valent tight-binding model for titanium. The effects of weak dissipation on the propagation of the auxiliary degrees of freedom for the chemical potential and on-site Hamiltonian matrix elements that is used to counteract the accumulation of numerical noise during trajectories was also investigated.

Extended Lagrangian Born-Oppenheimer molecular dynamics simulations of the shock-induced chemistry of phenylacetylene.

The initial chemical events that occur during the shock compression of liquid phenylacetylene have been investigated using self-consistent tight binding molecular dynamics simulations. The extended Lagrangian Born-Oppenheimer molecular dynamics formalism enabled us to compute microcanonical trajectories with precise conservation of the total energy. Our simulations revealed that the first density-increasing step under shock compression arises from the polymerization of phenylacetylene molecules at the acetylene moiety. The application of electronic structure-based molecular dynamics with long-term conservation of the total energy enabled us to identify electronic signatures of reactivity via monitoring changes in the HOMO-LUMO gap, and to capture directly adiabatic shock heating, transient non-equilibrium states, and changes in temperature arising from exothermic chemistry in classical molecular dynamics trajectories.

Canonical density matrix perturbation theory.

Density matrix perturbation theory [Niklasson and Challacombe, Phys. Rev. Lett. 92, 193001 (2004)] is generalized to canonical (NVT) free-energy ensembles in tight-binding, Hartree-Fock, or Kohn-Sham density-functional theory. The canonical density matrix perturbation theory can be used to calculate temperature-dependent response properties from the coupled perturbed self-consistent field equations as in density-functional perturbation theory. The method is well suited to take advantage of sparse matrix algebra to achieve linear scaling complexity in the computational cost as a function of system size for sufficiently large nonmetallic materials and metals at high temperatures.

Computation of the Density Matrix in Electronic Structure Theory in Parallel on Multiple Graphics Processing Units.

The algorithm developed in Cawkwell, M. J. et al. J. Chem. Theory Comput. 2012 , 8 , 4094 for the computation of the density matrix in electronic structure theory on a graphics processing unit (GPU) using the second-order spectral projection (SP2) method [ Niklasson, A. M. N. Phys. Rev. B 2002 , 66 , 155115 ] has been efficiently parallelized over multiple GPUs on a single compute node. The parallel implementation provides significant speed-ups with respect to the single GPU version with no loss of accuracy. The performance and accuracy of the parallel GPU-based algorithm is compared with the performance of the SP2 algorithm and traditional matrix diagonalization methods on a multicore central processing unit (CPU).

Energy conserving, linear scaling Born-Oppenheimer molecular dynamics.

Born-Oppenheimer molecular dynamics simulations with long-term conservation of the total energy and a computational cost that scales linearly with system size have been obtained simultaneously. Linear scaling with a low pre-factor is achieved using density matrix purification with sparse matrix algebra and a numerical threshold on matrix elements. The extended Lagrangian Born-Oppenheimer molecular dynamics formalism [A. M. N. Niklasson, Phys. Rev. Lett. 100, 123004 (2008)] yields microcanonical trajectories with the approximate forces obtained from the linear scaling method that exhibit no systematic drift over hundreds of picoseconds and which are indistinguishable from trajectories computed using exact forces.

Computing the Density Matrix in Electronic Structure Theory on Graphics Processing Units.

The self-consistent solution of a Schrödinger-like equation for the density matrix is a critical and computationally demanding step in quantum-based models of interatomic bonding. This step was tackled historically via the diagonalization of the Hamiltonian. We have investigated the performance and accuracy of the second-order spectral projection (SP2) algorithm for the computation of the density matrix via a recursive expansion of the Fermi operator in a series of generalized matrix-matrix multiplications. We demonstrate that owing to its simplicity, the SP2 algorithm [Niklasson, A. M. N. Phys. Rev. B2002, 66, 155115] is exceptionally well suited to implementation on graphics processing units (GPUs). The performance in double and single precision arithmetic of a hybrid GPU/central processing unit (CPU) and full GPU implementation of the SP2 algorithm exceed those of a CPU-only implementation of the SP2 algorithm and traditional matrix diagonalization when the dimensions of the matrices exceed about 2000 × 2000. Padding schemes for arrays allocated in the GPU memory that optimize the performance of the CUBLAS implementations of the level 3 BLAS DGEMM and SGEMM subroutines for generalized matrix-matrix multiplications are described in detail. The analysis of the relative performance of the hybrid CPU/GPU and full GPU implementations indicate that the transfer of arrays between the GPU and CPU constitutes only a small fraction of the total computation time. The errors measured in the self-consistent density matrices computed using the SP2 algorithm are generally smaller than those measured in matrices computed via diagonalization. Furthermore, the errors in the density matrices computed using the SP2 algorithm do not exhibit any dependence of system size, whereas the errors increase linearly with the number of orbitals when diagonalization is employed.