Maximum Entropy Theory of Multiscale Coarse-Graining via Matching Thermodynamic Forces: Application to a Molecular Crystal (TATB).

The MSCG/FM (multiscale coarse-graining via force-matching) approach is an efficient supervised machine learning method to develop microscopically informed coarse-grained (CG) models. We present a theory based on the principle of maximum entropy (PME) enveloping the existing MSCG/FM approaches. This theory views the MSCG/FM method as a special case of matching the thermodynamic forces from the extended ensemble described by the set of thermodynamic (relevant) system coordinates. This set may include CG coordinates, the stress tensor, applied external fields, and so forth, and may be characterized by nonequilibrium conditions. Following the presentation of the theory, we discuss the consistent matching of both bonded and nonbonded interactions. The proposed PME formulation is used as a starting point to extend the MSCG/FM method to the constant strain ensemble, which together with the explicit matching of the bonded forces is better suited for coarse-graining anisotropic media at a submolecular resolution. The theory is demonstrated by performing the fine coarse-graining of crystalline 1,3,5-triamino-2,4,6-trinitrobenzene (TATB), a well-known insensitive molecular energetic material, which exhibits highly anisotropic mechanical properties.

Hierarchical Machine Learning of Low-Resolution Coarse-Grained Free Energy Potentials.

A force-matching-based method for supervised machine learning (ML) of coarse-grained (CG) free energy (FE) potentials─known as multiscale coarse-graining via force-matching (MSCG/FM)─is an efficient method to develop microscopically informed CG models that are thermodynamically and statistically equivalent to the reference microscopic models. For low-resolution models, when the coarse-graining is at supramolecular scales, objective-oriented clustering of nonbonded particles is required and the reduced description becomes a function of the clustering algorithm. In the present work, we explore the dependence of the ML of the CG Helmholtz FE potential on the clustering algorithm. We consider coarse-graining based on partitional (k-means, leading to Voronoi diagram) and hierarchical agglomerative (bottom-up) clustering algorithms common in unsupervised ML and develop theory connecting the MSCG/FM learned CG Helmholtz potential and the clustering statistics. By combining the agglomerative clustering and the MSCG/FM learning in a recursive manner, we propose an efficient ML methodology to develop the fine-to-low resolution hierarchies of the CG models. The methodology does not suffer from degrading accuracy or increased computational cost to construct larger hierarchies and as such does not impose an upper size limitation of the CG particles resulting from the extended hierarchies. The utility of the methodology is demonstrated by obtaining the bottom-up agglomerative hierarchy for liquid nitromethane from all-atom molecular dynamics (MD) simulations. For agglomerative hierarchies, we prove the existence of renormalization group transformations that indicate self-similarity and allow for learning the low-resolution MSCG/FM potentials at low computational cost by rescaling and renormalizing the certain finer-resolution members of the hierarchy. The hierarchies of the CG models can be used to carry out simulations under constant-pressure conditions.

Building Chemical Property Models for Energetic Materials from Small Datasets Using a Transfer Learning Approach.

For many experimentally measured chemical properties that cannot be directly computed from first-principles, the existing physics-based models do not extrapolate well to out-of-sample molecules, and experimental datasets themselves are too small for traditional machine learning (ML) approaches. To overcome these limitations, we apply a transfer learning approach, whereby we simultaneously train a multi-target regression model on a small number of molecules with experimentally measured values and a large number of molecules with related computed properties. We demonstrate this methodology on predicting the experimentally measured impact sensitivity of energetic crystals, finding that both characteristics of the computed dataset and model architecture are important to prediction accuracy of the small experimental dataset. Our directed-message passing neural network (D-MPNN) ML model using transfer learning outperforms direct-ML and physics-based models on a diverse test set, and the new methods described here are widely applicable to modeling many other structure-property relationships.

Bottom-up coarse-grain modeling of plasticity and nanoscale shear bands in α-RDX.

Computationally inexpensive particle-based coarse-grained (CG) models are essential for use in molecular dynamics (MD) simulations of mesoscopically slow cooperative phenomena, such as plastic deformations in solids. Molecular crystals possessing complex symmetry present enormous practical challenges for particle-based coarse-graining at molecularly resolved scales, when each molecule is in a single-site representation, and beyond. Presently, there is no published pairwise non-bonded single-site CG potential that is able to predict the space group and structure of a molecular crystal. In this paper, we present a successful coarse-graining at a molecular level from first principles of an energetic crystal, hexahydro-1,3,5-trinitro-s-triazine (RDX) in the alpha phase, using the force-matching-based multiscale coarse-graining (MSCG/FM) approach. The new MSCG/FM model, which implements an optimal pair decomposition of the crystal Helmholtz free energy potential in molecular center-of-mass coordinates, was obtained by force-matching atomistic MD simulations of liquid, amorphous, and crystalline states and in a wide range of pressures (up to 20 GPa). The MSCG/FM potentials for different pressures underwent top-down optimization to fine-tune the mechanical and thermodynamic properties, followed by consolidation into a transferable density-dependent model referred to as RDX-TC-DD (RDX True-Crystal Density-Dependent). The RDX-TC-DD model predicts accurately the crystal structure of α-RDX at room conditions and reproduces the atomistic reference system under isothermal (300 K) hydrostatic compression up to 20 GPa, in particular, the Pbca symmetry of α-RDX in the elastic regime. The RDX-TC-DD model was then used to simulate the plastic response of uniaxially ([100]) compressed α-RDX resulting in nanoscale shear banding, a key mechanism for plastic deformation and defect-free detonation initiation proposed for many molecular crystalline explosives. Additionally, a comparative analysis of the effect of core-softening of the RDX-TC-DD potential and the degree of molecular rigidity in the all-atom treatment suggests a stress-induced short-range softening of the effective intermolecular interaction as a fundamental cause of plastic instability in α-RDX. The reported RDX-TC-DD model and overall workflow to develop it open up possibilities to perform high quality simulation studies of molecular energetic materials under thermal and mechanical stimuli, including extreme conditions.

Machine learning transition temperatures from 2D structure.

A priori knowledge of physicochemical properties such as melting and boiling could expedite materials discovery. However, theoretical modeling from first principles poses a challenge for efficient virtual screening of potential candidates. As an alternative, the tools of data science are becoming increasingly important for exploring chemical datasets and predicting material properties. Herein, we extend a molecular representation, or set of descriptors, first developed for quantitative structure-property relationship modeling by Yalkowsky and coworkers known as the Unified Physicochemical Property Estimation Relationships (UPPER). This molecular representation has group-constitutive and geometrical descriptors that map to enthalpy and entropy; two thermodynamic quantities that drive thermal phase transitions. We extend the UPPER representation to include additional information about sp2-bonded fragments. Additionally, instead of using the UPPER descriptors in a series of thermodynamically-inspired calculations, as per Yalkowsky, we use the descriptors to construct a vector representation for use with machine learning techniques. The concise and easy-to-compute representation, combined with a gradient-boosting decision tree model, provides an appealing framework for predicting experimental transition temperatures in a diverse chemical space. An application to energetic materials shows that the method is predictive, despite a relatively modest energetics reference dataset. We also report competitive results on diverse public datasets of melting points (i.e., OCHEM, Enamine, Bradley, and Bergström) comprised of over 47k structures. Open source software is available at https://github.com/USArmyResearchLab/ARL-UPPER.

Heuristics for chemical species identification in dense systems.

A new approach to identify chemical species from molecular dynamics (MD) simulations of reacting materials under extreme temperatures and pressures is presented. The approach is based on bond-distance and vibrational criteria, derived from the examination of atomic behavior during a density functional theory MD simulation of an overdriven shock of the explosive pentaerythritol tetranitrate. For comparison, the trajectory was analyzed using popular bonding criteria commonly used in analysis of reactive MD simulations, including distance, distance-time, and bond-order criteria. Cluster analyses using the new time-dependent bond definition approach presented here and a bond-order approach revealed that species and their corresponding lifetimes were strongly dependent on the chosen approach, indicating significant implications for the development of chemical mechanisms and chemical kinetics models using the results of reactive MD simulations.

Transferable Reactive Force Fields: Extensions of ReaxFF-lg to Nitromethane.

Transferable ReaxFF-lg models of nitromethane that predict a variety of material properties over a wide range of thermodynamic states are obtained by screening a library of ∼6600 potentials that were previously optimized through the Multiple Objective Evolutionary Strategies (MOES) approach using a training set that included information for other energetic materials composed of carbon, hydrogen, nitrogen, and oxygen. Models that best match experimental nitromethane lattice constants at 4.2 K and 1 atm are evaluated for transferability to high-pressure states at room temperature and are shown to better predict various liquid- and solid-phase structural, thermodynamic, and transport properties as compared to the existing ReaxFF and ReaxFF-lg parametrizations. Although demonstrated for an energetic material, the library of ReaxFF-lg models is supplied to the scientific community to enable new research explorations of complex reactive phenomena in a variety of materials research applications.

Theoretical Study of Shocked Formic Acid: Born-Oppenheimer MD Calculations of the Shock Hugoniot and Early-Stage Chemistry.

Quantum and classical molecular dynamics simulations are used to explore whether chemical reactivity of shocked formic acid occurs at pressures greater than 15 GPa, a question arising from results of different shock compression experiments. The classical molecular dynamics simulations were performed using a quantum-based nonreactive pair additive interaction potential whereas the full resolution quantum mechanical molecular dynamics simulations allow chemical reactions. Although the shock Hugoniot curve calculated using nonreactive classical MD for formic acid is in reasonable agreement with one set of experimental results, shock Hugoniot points calculated using Born-Oppenheimer MD at 30 GPa are in agreement with the set of experimental data that suggests chemical reactivity at these elevated temperatures and pressures. Examination of atomic positions throughout the Born-Oppenheimer MD trajectories clearly indicates extensive and complex chemical reaction, chiefly involving hydrogen-atom transfer and intermolecular complexation.

Parameterizing complex reactive force fields using multiple objective evolutionary strategies (MOES). Part 1: ReaxFF models for cyclotrimethylene trinitramine (RDX) and 1,1-diamino-2,2-dinitroethene (FOX-7).

ReaxFF (van Duin, A.C.T.; Dasgupta, S.; Lorant, F.; Goddard, W.A. J. Phys. Chem. A, 2001, 105, 9396-9409) reactive potentials are parametrized for cyclotrimethylene trinitramine (RDX) and 1,1-diamino-2,2-dinitroethene (FOX-7) in a novel application combining data envelopment analysis and a modern self-adaptive evolutionary algorithm to optimize multiple objectives simultaneously and map the entire family of solutions. In order to correct the poor crystallographic parameters predicted by ReaxFF using its base parametrization (Strachan, A.; van Duin, A. C. T.; Chakraborty, D.; Dasgupta S.; Goddard, W. A. Phys. Rev. Lett., 2003, 91, 098301), we augmented the existing training set data used for parametrization with additional (SAPT)DFT calculations of RDX and FOX-7 dimer interactions. By adjusting a small subset of the ReaxFF parameters that govern long-range interactions, the evolutionary algorithm approach converges on a family of solutions that best describe crystallographic parameters through simultaneous optimization of the objective functions. Molecular dynamics calculations of RDX and FOX-7 are conducted to assess the quality of the force fields, resulting in parametrizations that improve the overall prediction of the crystal structures.

Parameterizing complex reactive force fields using multiple objective evolutionary strategies (MOES): Part 2: transferability of ReaxFF models to C-H-N-O energetic materials.

The Multiple Objective Evolutionary Strategies (MOES) algorithm was used to parametrize force fields having the form of the reactive models ReaxFF (van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396) and ReaxFF-lg (Liu, L.; Liu, Y.; Zybin, S. V.; Sun, H.; Goddard, W. A. J. Phys. Chem. A 2011, 115, 11016) in an attempt to produce equal or superior ambient state crystallographic structural results for cyclotrimethylene trinitramine (RDX). Promising candidates were then subjected to molecular dynamics simulations of five other well-known conventional energetic materials to assess the degree of transferability of the models. Two models generated through the MOES search were shown to have performance better than or as good as ReaxFF-lg in describing the six energetic systems modeled. This study shows that MOES is an effective and efficient method to develop complex force fields.

On the importance of shear dissipative forces in coarse-grained dynamics of molecular liquids.

In this work we demonstrate from first principles that the shear frictions describing dissipative forces in the direction normal to the vector connecting the coarse-grained (CG) particles in dissipative particle dynamics (DPD) could be dominant for certain real molecular liquids at high-resolution coarse-graining. This is in contrast to previous works on bottom-up DPD modeling and indicates that such liquids cannot be simulated accurately using the conventional form of DPD which relies only on frictions in the radial direction. Specifically, we describe the development of fully bottom-up CG models for liquid hexahydro-1,3,5-trinitro-s-triazine (RDX) which are incorporated into the DPD method. Consistent with the microscopic foundation of DPD dynamics, the conservative part of the DPD models is obtained by the multi-scale coarse-graining (MS-CG) approach, which implements the pairwise decomposition of the atomistic potential of mean force (PMF) in CG coordinates. The radial and shear distant-dependent friction coefficients in a parameter-free form are derived systematically from microscopic velocity and force correlation data along system trajectories using a recently proposed approach [J. Chem. Phys., 2014, 140, 104104]. The shear dissipative forces for the reported system appear to be dominant. We discuss the implications of dominant shear dissipation on dynamical and transport properties of CG liquids such as diffusion and viscosity as revealed by simulations of liquid RDX using the new MS-CG/DPD models.

A new parameter-free soft-core potential for silica and its application to simulation of silica anomalies.

A core-softening of the effective interaction between oxygen atoms in water and silica systems and its role in developing anomalous thermodynamic, transport, and structural properties have been extensively debated. For silica, the progress with addressing these issues has been hampered by a lack of effective interaction models with explicit core-softening. In this work, we present an extension of a two-body soft-core interatomic force field for silica recently reported by us [S. Izvekov and B. M. Rice, J. Chem. Phys. 136(13), 134508 (2012)] to include three-body forces. Similar to two-body interaction terms, the three-body terms are derived using parameter-free force-matching of the interactions from ab initio MD simulations of liquid silica. The derived shape of the O-Si-O three-body potential term affirms the existence of repulsion softening between oxygen atoms at short separations. The new model shows a good performance in simulating liquid, amorphous, and crystalline silica. By comparing the soft-core model and a similar model with the soft-core suppressed, we demonstrate that the topology reorganization within the local tetrahedral network and the O-O core-softening are two competitive mechanisms responsible for anomalous thermodynamic and kinetic behaviors observed in liquid and amorphous silica. The studied anomalies include the temperature of density maximum locus and anomalous diffusivity in liquid silica, and irreversible densification of amorphous silica. We show that the O-O core-softened interaction enhances the observed anomalies primarily through two mechanisms: facilitating the defect driven structural rearrangements of the silica tetrahedral network and modifying the tetrahedral ordering induced interactions toward multiple characteristic scales, the feature which underlies the thermodynamic anomalies.

Multi-scale coarse-graining of non-conservative interactions in molecular liquids.

A new bottom-up procedure for constructing non-conservative (dissipative and stochastic) interactions for dissipative particle dynamics (DPD) models is described and applied to perform hierarchical coarse-graining of a polar molecular liquid (nitromethane). The distant-dependent radial and shear frictions in functional-free form are derived consistently with a chosen form for conservative interactions by matching two-body force-velocity and three-body velocity-velocity correlations along the microscopic trajectories of the centroids of Voronoi cells (clusters), which represent the dissipative particles within the DPD description. The Voronoi tessellation is achieved by application of the K-means clustering algorithm at regular time intervals. Consistently with a notion of many-body DPD, the conservative interactions are determined through the multi-scale coarse-graining (MS-CG) method, which naturally implements a pairwise decomposition of the microscopic free energy. A hierarchy of MS-CG/DPD models starting with one molecule per Voronoi cell and up to 64 molecules per cell is derived. The radial contribution to the friction appears to be dominant for all models. As the Voronoi cell sizes increase, the dissipative forces rapidly become confined to the first coordination shell. For Voronoi cells of two and more molecules the time dependence of the velocity autocorrelation function becomes monotonic and well reproduced by the respective MS-CG/DPD models. A comparative analysis of force and velocity correlations in the atomistic and CG ensembles indicates Markovian behavior with as low as two molecules per dissipative particle. The models with one and two molecules per Voronoi cell yield transport properties (diffusion and shear viscosity) that are in good agreement with the atomistic data. The coarser models produce slower dynamics that can be appreciably attributed to unaccounted dissipation introduced by regular Voronoi re-partitioning as well as by larger numerical errors in mapping out the dissipative forces. The framework presented herein can be used to develop computational models of real liquids which are capable of bridging the atomistic and mesoscopic scales.

Assessing the Performances of Dispersion-Corrected Density Functional Methods for Predicting the Crystallographic Properties of High Nitrogen Energetic Salts.

Several density functional methods with corrections for long-range dispersion interactions are evaluated for their capabilities to describe the crystallographic lattice properties of a set of 26 high nitrogen-content salts relevant for energetic materials applications. Computations were done using methods that ranged from adding atom-atom dispersion corrections with environment-independent and environment-dependent coefficients, to methods that incorporate dispersion effects via dispersion-corrected atom-centered potentials (DCACP), to methods that include nonlocal corrections. Among the functionals tested, the most successful is the nonlocal optPBE-vdW functional of Klimes and Michaelides that predicts unit cell volumes for all crystals of the reference set within the target error range of ±3% and gives individual lattice parameters with a mean average percent error of less than 0.81%. The DCACP, Grimme's D3, and Becke and Johnson's exchange-hole (XDM) methods, when used with the BLYP, PBE, and B86b functionals, respectively, are also quite successful at predicting the lattice parameters of the test set.

Evaluation of electrostatic descriptors for predicting crystalline density.

This study evaluates the importance of electrostatic corrections to earlier quantum-mechanically based methods to predict crystal densities of neutral and ionic molecular energetic materials. Our previous methods (B. M. Rice et al., J. Phys. Chem. A 2007, 111, 10874) use the molecular volumes of the isolated molecule or formula unit to estimate the crystal density; this volume is defined to be that inside the quantum-mechanically determined 0.001 a.u. isosurface of electron density surrounding the isolated molecule. The electrostatic corrections to these volumetric estimates are based on features of the electrostatic potential mapped onto this isosurface of electron density, and have been parameterized using information from 180 neutral and 23 ionic CHNO molecular systems. The quality of the electrostatically corrected methods was assessed through application to 38 neutral and 48 ionic compounds not used in the parameterization. The root mean square (rms) percent deviation and average absolute error of predictions for the 38 neutral species relative to experiment are 2.7% and 0.035 g/cm(3), respectively, decreases of 0.9% and 0.015 g/cm(3) from the earlier predictions (3.6% and 0.050 g/cm(3), respectively). The rms percent deviation and average absolute error of predictions for the 48 ionic compounds relative to experiment are 3.7% and 0.045 g/cm(3), respectively, decreases of 2.6% and 0.043 g/cm(3) from the earlier predictions that used the formula unit volumes only. The results clearly show a significant improvement to the earlier method upon inclusion of electrostatic corrections.

New form of polymeric nitrogen from dynamic shock simulation.

For many years there has been significant interest in polymeric phases of nitrogen at low pressure for potential application as an energetic material. This was the result of years of theoretical work indicating potentially meta-stable polymeric nitrogen. Experimental evidence of both an amorphous phase and a cubic-gauche phase has added greatly to this interest [A. F. Goncharov, E. A. Gregoryanz, H. K. Mao, Z. Liu, and R. J. Hemley, Phys. Rev. Lett. 85, 1262 (2000); M. I. Eremets, R. J. Hemley, H. K. Mao, and E. Gregoryanz, Nature (London) 411, 170 (2001)]. While most of the theoretical work has been done on the many crystal phases of nitrogen, little work has been done on simulating amorphous polymeric nitrogen. The original goal of this work was to simulate amorphous polymeric nitrogen at low pressure; however, we unexpectedly found a new form of polymeric nitrogen. Starting from first principles dynamic shock simulation of cubic-gauche nitrogen [W. D. Mattson and R. Balu, Phys. Rev. B 83, 174105 (2011)] we demonstrate a new low pressure porous form that exhibits stability at low temperatures. We describe the detailed procedure of obtaining this structure as well as some of its physical characteristics. Finally, we explore composite structures of this new form of polymeric nitrogen and their possible relationship to an amorphous form.

Shock Hugoniot calculations of polymers using quantum mechanics and molecular dynamics.

Using quantum mechanics (QM) and classical force-field based molecular dynamics (FF), we have calculated the principle shock Hugoniot curves for numerous amorphous polymers including poly[methyl methacrylate] (PMMA), poly[styrene], polycarbonate, as well as both the amorphous and crystalline forms of poly[ethylene]. In the FF calculations, we considered a non-reactive force field (i.e., polymer consistent FF). The QM calculations were performed with density functional theory (DFT) using dispersion corrected atom centered pseudopotentials. Overall, results obtained by DFT show much better agreement with available experimental data than classical force fields. In particular, DFT calculated Hugoniot curves for PMMA up to 74 GPa are in very good agreement with experimental data, where a preliminary study of chain fracture and association was also performed. Structure analysis calculations of the radius of gyration and carbon-carbon radial distribution function were also carried out to elucidate contraction of the polymer chains with increasing pressure.

Free-energy based pair-additive potentials for bulk Ni-Al systems: application to study Ni-Al reactive alloying.

We present new numerical pair-additive Al, Ni, and Al-Ni potentials by force-matching (FM) ionic force and virial data from single (bulk liquid) phase ab initio molecular dynamics (MD) simulations using the Born-Oppenheimer method. The potentials are represented by piece-wise functions (splines) and, therefore, are not constrained to a particular choice of analytical functional form. The FM method with virial constraint naturally yields a potential which maps out the ionic free-energy surface of the reference ensemble. To further improve the free energetics of the FM ensemble, the FM procedure is modified to bias the potentials to reproduce the experimental melting temperatures of the reference (FCC-Al, FCC-Ni, B2-NiAl) phases, the only macroscopic data included in the fitting set. The performance of the resultant potentials in simulating bulk metallic phases is then evaluated. The new model is applied to perform MD simulations of self-propagating exothermic reaction in Ni-Al bilayers at P = 0-5 GPa initiated at T = 1300 K. Consistent with experimental observations, the new model describes realistically a sequence of peritectic phase transformations throughout the reaction and at a realistic rate. The reaction proceeds through interlayer diffusion of Al and Ni atoms at the interface with formation of B2-NiAl in the Al melt. Such material responses have, in the past, been proven to be difficult to observe with then-existing potentials.

Mechanism of densification in silica glass under pressure as revealed by a bottom-up pairwise effective interaction model.

A new short-range pairwise numerical potential for silica is presented. The potential is derived from a single ab initio molecular dynamics (AIMD) simulation of molten silica using the force-matching method with the forces being represented numerically by piecewise functions (splines). The AIMD simulation is performed using the Born-Oppenheimer method with the generalized gradient approximation (BLYP) for the XC energy functional. The new effective potential includes a soft-repulsive shoulder to describe the interactions of oxygen ions at short separations. The new potential, despite being short-ranged and derived from single-phase data, exhibits a good transferability to silica crystalline polymorphs and amorphous silica. The importance of the O-O soft-repulsive shoulder interaction on glass densification under cold and shock compressions is assessed from MD simulations of silica glass under room and shock Hugoniot conditions, respectively. Results from these simulations indicate that the appearance of oxygen complexes (primarily pairs) interacting through soft-repulsive shoulder potential occurs at 8-10 GPa, and under cold compression conditions becomes notable at 40 GPa, essentially coinciding with the transition to a Si sixfold coordination state. An analysis of changes in system structure in compressed and shocked states reveals that the O ions interacting through the soft-repulsive shoulder potential in denser states of silica glass may create a mechanical multi-stability under elevated pressures and thus to contribute to the observed anomalous densification.