Classical Reactive Molecular Dynamics Implementations: State of the Art

Study on Thermophysical Properties and Phase Change Regulation Mechanism of Optically-Controlled Phase Change Materials: Synthesis, Crystal Structure and Molecular Dynamics

Optically-controlled phase change materials, which are prepared by introducing molecular photoswitches into traditional phase change materials (PCMs), can convert and store solar energy into photochemical enthalpy and phase change enthalpy. However, the thermophysical properties of optically controlled PCMs, which are crucial in the practical, are rarely paid attention to. 4-(phenyldiazenyl)phenyl decanoate (Azo-A-10) is experimentally prepared as an optically-controlled PCMs, whose energy storage density is 210.0 kJ·kg-1, and the trans single crystal structure is obtained. The density, phase transition temperature, thermal conductivity, and other parameters in trans state are measured experimentally. Furthermore, a microscopic model of Azo-A-10 is established, and the thermophysical properties are analyzed based on molecular dynamics. The results show that the microstructure parameter (order parameters) and thermophysical properties (density, radial distribution function, self-diffusion coefficient, phase change temperature, and thermal conductivity) of partially or completely isomerized Azo-A-10, which are challenging to observe in experiments, can be predicted by molecular dynamics simulation. The optically-controlled phase change mechanism can be clarified according to the differences in microstructure. The optically-controlled switchability of thermophysical properties of an optically-controlled PCM is analyzed. This study provides ideas for the improvement, development, and application of optically-controlled PCMs in the future.

Surface-Initiated Polymerization with an Initiator Gradient: A Monte Carlo Simulation

Due to the difficulty of accurately characterizing properties such as the molecular weight (Mn) and grafting density (σ) of gradient brushes (GBs), these properties are traditionally assumed to be uniform in space to simplify analysis. Applying a stochastic reaction model (SRM) developed for heterogeneous polymerizations, we explored surface-initiated polymerizations (SIPs) with initiator gradients in lattice Monte Carlo simulations to examine this assumption. An initial exploration of SIPs with 'homogeneously' distributed initiators revealed that increasing σ slows down the polymerization process, resulting in polymers with lower molecular weight and larger dispersity (Đ) for a given reaction time. In SIPs with an initiator gradient, we observed that the properties of the polymers are position-dependent, with lower Mn and larger Đ in regions of higher σ, indicating the non-uniform properties of polymers in GBs. The results reveal a significant deviation in the scaling behavior of brush height with σ compared to experimental data and theoretical predictions, and this deviation is attributed to the non-uniform Mn and Đ.

Updates to Protex for Simulating Proton Transfers in an Ionic Liquid

The Python-based program Protex was initially developed for simulating proton transfers in a pure protic ionic liquid via polarizable molecular dynamics simulations. This method employs a single topology approach wherein deprotonated species retain a dummy atom, which is transformed into a real hydrogen atom during the protonation process. In this work, we extended Protex to include more intricate systems and to facilitate the simulation of the Grotthuss mechanism to enhance alignment with the empirical findings. The handling of proton transfer events within Protex was further refined for increased flexibility. In the original model, each deprotonated molecule contained a single dummy atom connected to the hydrogen acceptor atom. This model posed limitations for molecules with multiple atoms that could undergo protonation. To mitigate this issue, Protex was extended to execute a proton transfer when one of these potential atoms was within a suitable proximity for the transfer event. For the purpose of maintaining simplicity, Protex continues to utilize only a single dummy atom per deprotonated molecule. Another new feature pertains to the determination of the eligibility for a proton transfer event. A range of acceptable distances can now be defined within which the transfer probability is gradually turned off. These modifications allow for a more nuanced approach to simulating proton transfer events, offering greater accuracy and control of the modeling process.

Molecular Dynamic Simulations of Aqueous Micellar Organometallic Catalysis: Methane Functionalization as a Case Study

Molecular Dynamics (MD) simulations constitute a powerful tool that provides a 3D perspective of the dynamical behavior of chemical systems. Herein the first MD study of the dynamics of a catalytic organometallic system, in micellar media, is presented. The challenging methane catalytic functionalization into ethyl propionate through a silver-catalyzed process has been targeted as the case study. The intimate nature of the micelles formed with the surfactants sodium dodecylsulfate (SDS) and potassium perfluorooctane sulfonate (PFOS) has been ascertained, as well as the relative distribution of the main actors in this transformation, namely methane, the diazo reagent and the silver catalyst, the latter in two different forms: the initial compound and a silver-carbene intermediate. Catalyst deactivation occurs with halide containing surfactants dodecyltrimethylammonium chloride (DTAC) and Triton X-100. Computed simulations allow explaining the experimental results, indicating that micelles behave differently regarding the degree of accumulation and the local distribution of the reactants and their effect in the molecular collisions leading to net reaction.

Elastocaloric Effect in Graphene Kirigami

Kirigami, a traditional Japanese art of paper cutting, has recently been explored for its elastocaloric effect (ECE) in kirigami-based materials (KMs), where an applied strain induces temperature changes. Importantly, the feasibility of a nanoscale graphene kirigami monolayer was experimentally demonstrated. Here, we investigate the ECE in GK representing the thinnest possible KM to better understand this phenomenon. Through molecular dynamics simulations, we analyze the temperature change and coefficient of performance (COP) of GK. Our findings reveal that while GKs lack the intricate temperature changes observed in macroscopic KMs, they exhibit a substantial temperature change of approximately 9.32 K (23 times higher than that of macroscopic KMs, which is about 0.4 K) for heating and -3.50 K for cooling. Furthermore, they demonstrate reasonable COP values of approximately 1.57 and 0.62, respectively. It is noteworthy that the one-atom-thick graphene configuration prevents the occurrence of the complex temperature distribution observed in macroscopic KMs.

Molecular Dynamics Simulations of Ideal Living Polymerization: Terminal Model and Kinetic Aspects

Living polymerization is an important synthetic approach to achieving precise control of synthesized polymers, which is crucial for their applications. The molecular weight distribution (MWD) prescribes the macroscopic properties of polymers and hence is a key feature to characterize polymerization. In this work, we present a systematic molecular dynamics simulation study of ideal living polymerization in bulk and surface-initiated systems based on a terminal stochastic reaction model. The evolution of polymer dispersity and MWD along with the polymerization process is examined. We demonstrate that MWD is generally well captured by the Schulz-Zimm distribution for bulk and surface-initiated systems with low grafting densities. However, as the grafting density in the surface-initiated case increases, heterogeneity in chain growth emerges due to the kinetic trapping of reactive sites, which causes the starving of short chains and the thriving of minority long chains such that a shoulder region shows up in MWD. This effect can be enhanced by kinetic compressing induced by polymerization. In addition, the interplay of bonding reaction kinetics and other kinetic properties (e.g., mass transfer and polymer relaxation) is further explored, alongside the influences of bonding probability and reactant concentration. We expect that this investigation will aid in our understanding of typical kinetic aspects of living polymerization.

A signal "on-off-on"-type electrochemiluminescence aptamer sensor for detection of sulfadimethoxine based on Ru@Zn-oxalate MOF composites

An "on-off-on"-type electrochemiluminescence (ECL) aptamer sensor based on Ru@Zn-oxalate metal-organic framework (MOF) composites is constructed for sensitive detection of sulfadimethoxine (SDM). The prepared Ru@Zn-oxalate MOF composites with the three-dimensional structure provide good ECL performance for the "signal-on." The MOF structure with a large surface area enables the material to fix more Ru(bpy)₃²⁺. Moreover, the Zn-oxalate MOF with three-dimensional chromophore connectivity provides a medium which can accelerate excited-state energy transfer migration among Ru(bpy)₃²⁺ units, and greatly reduces the influence of solvent on chromophore, achieving a high-energy Ru emission efficiency. The aptamer chain modified with ferrocene at the end can hybridize with the capture chain DNA1 fixed on the surface of the modified electrode through base complementary pairing, which can significantly quench the ECL signal of Ru@Zn-oxalate MOF. SDM specifically binds to its aptamer to separate ferrocene from the electrode surface, resulting in a "signal-on" ECL signal. The use of the aptamer chain further improves the selectivity of the sensor. Thus, high-sensitivity detection of SDM specificity is realized through the specific affinity between SDM and its aptamer. This proposed ECL aptamer sensor has good analytical performance for SDM with low detection limit (27.3 fM) and wide detection range (100 fM-500 nM). The sensor also shows excellent stability, selectivity, and reproducibility, which proved its analytical performance. The relative standard deviation (RSD) of SDM detected by the sensor is between 2.39 and 5.32%, and the recovery is in the range 97.23 to 107.5%. The sensor shows satisfactory results in the analysis of actual seawater samples, which is expected to play a role in the exploration of marine environmental pollution.

Species-selective nanoreactor molecular dynamics simulations based on linear-scaling tight-binding quantum chemical calculations

Here, extensions to quantum chemical nanoreactor molecular dynamics simulations for discovering complex reactive events are presented. The species-selective algorithm, where the nanoreactor effectively works for the selected desired reactants, was introduced to the original scheme. Moreover, for efficient simulations of large model systems with the modified approach, the divide-and-conquer linear-scaling density functional tight-binding method was exploited. Two illustrative applications of the polymerization of propylene and cyclopropane mixtures and the aggregation of sodium chloride from aqueous solutions indicate that species-selective quantum chemical nanoreactor molecular dynamics is a promising method to accelerate the sampling of multicomponent chemical processes proceeding under relatively mild conditions.

Metal-organic frameworks for active food packaging. A review

Food wastage is a major concern for sustainable health and agriculture. To reduce food waste, classical preservation techniques such as drying, pasteurization, freeze-drying, fermentation, and microwave are available. Nonetheless, these techniques display shortcomings such as alteration of food and taste. Such shortcomings may be solved by active food packaging, which involves the incorporation of active agents into the packaging material. Recently, metal-organic frameworks, a class of porous hybrid supramolecular materials, have been developed as an active agent to extend food shelf life and maintain safety. Here, we review metal-organic frameworks in active packaging as oxygen scavengers, antimicrobials, moisture absorbers, and ethylene scavengers. We present methods of incorporation of metal-organic frameworks into packaging materials and their applications.

Controlling solvent quality by time: Self-avoiding sprints in nonequilibrium polymerization

A fundamental paradigm in polymer physics is that macromolecular conformations in equilibrium can be described by universal scaling laws, being key for structure, dynamics, and function of soft (biological) matter and in the materials sciences. Here we reveal that during diffusion-influenced, nonequilibrium chain-growth polymerization, scaling laws change qualitatively, in particular, the growing polymers exhibit a surprising self-avoiding walk behavior in poor and θ solvents. Our analysis, based on monomer-resolved, off-lattice reaction-diffusion computer simulations, demonstrates that this phenomenon is a result of (i) nonequilibrium monomer density depletion correlations around the active polymerization site, leading to a locally directed and self-avoiding growth, in conjunction with (ii) chain (Rouse) relaxation times larger than the competing polymerization reaction time. These intrinsic nonequilibrium mechanisms are facilitated by fast and persistent reaction-driven diffusion ("sprints") of the active site, with analogies to pseudochemotactic active Brownian particles. Our findings have implications for time-controlled structure formation in polymer processing, as in, e.g., reactive self-assembly, photocrosslinking, and three-dimensional printing.

Adaptive partitioning molecular dynamics using an extended Hamiltonian approach

A recently proposed extended Hamiltonian approach to switching interaction potentials is generalized to enable adaptive partitioning molecular dynamics simulations. Switching is performed along a fictitious classical degree of freedom whose value determines the mixing ratio of the two potentials on a time scale determined by its associated mass. We propose to choose this associated fictitious mass adaptively so as to ensure a constant time scale for all switching processes. For different model systems, including a harmonic oscillator and a Lennard-Jones fluid, we investigate the window of switching time scales that guarantees the conservation of the extended Hamiltonian for a large number of switching events. The methodology is first applied in the microcanonical ensemble and then generalized to the canonical ensemble using a Nosé-Hoover chain thermostat. It is shown that the method is stable for thousands of consecutive switching events during a single simulation, with constant temperature and a conserved extended Hamiltonian. A slight modification of the original Hamiltonian is introduced to avoid accumulation of small numerical errors incurred after each switching process.

Learning heterogenous reaction rates from stochastic simulations

Reaction rate equations are ordinary differential equations that are frequently used to describe deterministic chemical kinetics at the macroscopic scale. At the microscopic scale, the chemical kinetics is stochastic and can be captured by complex dynamical systems reproducing spatial movements of molecules and their collisions. Such molecular dynamics systems may implicitly capture intricate phenomena that affect reaction rates but are not accounted for in the macroscopic models. In this work we present a data assimilation procedure for learning nonhomogeneous kinetic parameters from molecular simulations with many simultaneously reacting species. The learned parameters can then be plugged into the deterministic reaction rate equations to predict long time evolution of the macroscopic system. In this way, our procedure discovers an effective differential equation for reaction kinetics. To demonstrate the procedure, we upscale the kinetics of a molecular system that forms a complex covalently bonded network severely interfering with the reaction rates. Incidentally, we report that the kinetic parameters of this system feature peculiar time and temperature dependences, whereas the probability of a network strand to close a cycle follows a universal distribution.

Machine Learning for Chemical Reactions

Machine learning (ML) techniques applied to chemical reactions have a long history. The present contribution discusses applications ranging from small molecule reaction dynamics to computational platforms for reaction planning. ML-based techniques can be particularly relevant for problems involving both computation and experiments. For one, Bayesian inference is a powerful approach to develop models consistent with knowledge from experiments. Second, ML-based methods can also be used to handle problems that are formally intractable using conventional approaches, such as exhaustive characterization of state-to-state information in reactive collisions. Finally, the explicit simulation of reactive networks as they occur in combustion has become possible using machine-learned neural network potentials. This review provides an overview of the questions that can and have been addressed using machine learning techniques, and an outlook discusses challenges in this diverse and stimulating field. It is concluded that ML applied to chemistry problems as practiced and conceived today has the potential to transform the way with which the field approaches problems involving chemical reactions, in both research and academic teaching.

Connecting the quantum and classical mechanics simulation world: Applications of reactive step molecular dynamics simulations

This article presents the application of the reactive step molecular dynamics simulation method [M. Biedermann, D. Diddens, and A. Heuer, J. Chem. Theory Comput. 17, 1074 (2021)] toward two different atomistic, chemically reactive systems. During reactive steps, transitions from reactant to product molecules are modeled according to physically correct transition probabilities based on quantum chemical information about the reactions such as molecular reaction rates via instant exchange of the employed force field and a subsequent, short relaxation of the structure. In the first application, we study the follow-up reactions of singly reduced ethylene carbonate (EC) radicals in EC solution, first, via extensive ab initio molecular dynamics simulations and, second, with the reactive step algorithm. A direct comparison of both simulation methods shows excellent agreement. Then, we employ the reactive step algorithm to simulate the enolate formation of 2-methylcyclopropanone with the base lithium diisopropylamine. Thereby, we can demonstrate that the reactive step algorithm is also capable of capturing effects from kinetic vs thermodynamic control of chemical reactions during simulation.

Reactive Molecular Dynamics at Constant Pressure via Nonreactive Force Fields: Extending the Empirical Valence Bond Method to the Isothermal-Isobaric Ensemble

The Empirical Valence Bond (EVB) method offers a suitable framework to obtain reactive potentials through the coupling of nonreactive force fields. In this formalism, most of the implemented coupling terms are built using functional forms that depend on spatial coordinates, while parameters are fitted against reference data to model the change of chemistry between the participating nonreactive states. In this work, we demonstrate that the use of such coupling terms precludes the computation of the stress tensor for condensed phase systems and prevents the possibility to carry out EVB molecular dynamics in the isothermal-isobaric (NPT) ensemble. Alternatively, we make use of coupling terms that depend on the energy gaps, defined as the energy differences between the participating nonreactive force fields, and derive a general expression for the EVB stress tensor suitable for computation. Implementation of this new methodology is tested for a model of a single reactive malonaldehyde solvated in nonreactive water. Mass densities and probability distributions for the values of the energy gaps computed in the NPT ensemble reveal a negligible role of the reactive potential in the limit of low concentrated solutions, thus corroborating for the first time the validity of approximations based on the canonical NVT ensemble, customarily adopted for EVB simulations. The presented formalism also aims to contribute to future implementations and extensions of the EVB method to research the limit of highly concentrated solutions.

Ten quick tips for homology modeling of high-resolution protein 3D structures

The purpose of this quick guide is to help new modelers who have little or no background in comparative modeling yet are keen to produce high-resolution protein 3D structures for their study by following systematic good modeling practices, using affordable personal computers or online computational resources. Through the available experimental 3D-structure repositories, the modeler should be able to access and use the atomic coordinates for building homology models. We also aim to provide the modeler with a rationale behind making a simple list of atomic coordinates suitable for computational analysis abiding to principles of physics (e.g., molecular mechanics). Keeping that objective in mind, these quick tips cover the process of homology modeling and some postmodeling computations such as molecular docking and molecular dynamics (MD). A brief section was left for modeling nonprotein molecules, and a short case study of homology modeling is discussed.

Reactive Molecular Dynamics Study of the Thermal Decomposition of Phenolic Resins

The thermal decomposition of polyphenolic resins was studied by reactive molecular dynamics (RMD) simulation at elevated temperatures. Atomistic models of the polyphenolic resins to be used in the RMD were constructed using an automatic method which calls routines from the software package Materials Studio. In order to validate the models, simulated densities and heat capacities were compared with experimental values. The most suitable combination of force field and thermostat for this system was the Forcite force field with the Nosé–Hoover thermostat, which gave values of heat capacity closest to those of the experimental values. Simulated densities approached a final density of 1.05–1.08 g/cm3 which compared favorably with the experimental values of 1.16–1.21 g/cm3 for phenol-formaldehyde resins. The RMD calculations were run using LAMMPS software at temperatures of 1250 K and 3000 K using the ReaxFF force field and employing an in-house routine for removal of products of condensation. The species produced during RMD correlated with those found experimentally for polyphenolic systems and rearrangements to form cyclopropane moieties were observed. At the end of the RMD simulations a glassy carbon char resulted.

Low temperature reaction dynamics for CH3OH + OH collisions on a new full dimensional potential energy surface

Is the rise of the rate constant measured in laval expansion experiments of OH with organic molecules at low temperatures due to the reaction between the reactants or due to the formation of complexes with the buffer gas? This question has importance for understanding the evolution of prebiotic molecules observed in different astrophysical objects. Among these molecules methanol is one of the most widely observed, and its reaction with OH has been studied by several groups showing a fast increase in the rate constant under 100 K. Transition state theory doesn't reproduce this behavior and here dynamical calculations are performed on a new full dimensional potential energy surface developed for this purpose. The calculated classical reactive cross sections show an increase at low collision energies due to a complex forming mechanism. However, the calculated rate constant at temperatures below 100 K remains lower than the observed one. Quantum effects are likely responsible for the measured behavior at low temperatures.

Docking of Covalent Ligands: Challenges and Approaches

Covalent ligands have recently regained considerable attention in drug discovery. The rational design of such ligands, however, is still faced with particular challenges, mostly related to the fact that covalent bond formation is a quantum mechanical phenomenon which cannot adequately be handled by the force fields or empirical approaches typically used for noncovalent protein-ligand interactions. Although the necessity for quantum chemical approaches is clear, they cannot yet routinely be applied on large data sets of ligands or for a broader exploration of binding modes in docking calculations. On the other hand, technical solutions for performing docking calculations with covalent ligands are available, but their scope is normally quite limited. Scoring functions typically neglect the contribution from covalent bond formation completely. In this situation, the question arises how to approach covalent ligands and which methods to choose for their docking and design.

Modeling the free-radical polymerization of hexanediol diacrylate (HDDA): a molecular dynamics and graph theory approach

In the printing, coating and ink industries, photocurable systems are becoming increasingly popular and multi-functional acrylates are one of the most commonly used monomers due to their high reactivity (fast curing). In this paper, we use molecular dynamics and graph theory tools to investigate the thermo-mechanical properties and topology of hexanediol diacrylate (HDDA) polymer networks. The gel point was determined as the point where a giant component was formed. For the conditions of our simulations, we found the gel point to be around 0.18 bond conversion. A detailed analysis of the network topology showed, unexpectedly, that the flexibility of the HDDA molecules plays an important role in increasing the conversion of double bonds, while delaying the gel point. This is due to a back-biting type of reaction mechanism that promotes the formation of small cycles. The glass transition temperature for several degrees of curing was obtained from the change in the thermal expansion coefficient. For a bond conversion close to experimental values we obtained a glass transition temperature around 400 K. For the same bond conversion we estimate a Young's modulus of 3 GPa. Both of these values are in good agreement with experiments.

Full dimensional potential energy surface and low temperature dynamics of the H2CO + OH → HCO + H2O reaction

A new method is proposed to analytically represent the potential energy surface of reactions involving polyatomic molecules capable of accurately describing long-range interactions and saddle points, needed to describe low-temperature collisions. It is based on two terms, a reactive force field term and a many-body term. The reactive force field term accurately describes the fragments, long-range interactions among them and the saddle points for reactions. The many-body term increases the desired accuracy everywhere else. This method has been applied to the OH + H2CO → H2O + HCO reaction, giving a barrier of 27.4 meV. The simulated classical rate constants with this potential are in good agreement with recent experimental results [Ocaña et al., Astrophys. J., 2017, submitted], showing an important increase at temperatures below 100 K. The reaction mechanism is analyzed in detail here, and explains the observed behavior at low energy by the formation of long-lived collision complexes, with roaming trajectories, with a capture observed for very long impact parameters, >100 a.u., determined by the long-range dipole-dipole interaction.

Comparing the Mechanical Response of Di-, Tri-, and Tetra-functional Resin Epoxies with Reactive Molecular Dynamics

The influence of monomer functionality on the mechanical properties of epoxies is studied using Molecular Dynamics (MD) with the Reax Force Field (ReaxFF). From deformation simulations, the Young's modulus, yield point, and Poisson's ratio are calculated and analyzed. The results demonstrate an increase in stiffness and yield strength with increasing resin functionality. Comparison between the network structures of distinct epoxies is further advanced by the Monomeric Degree Index (MDI). Experimental validation demonstrates the MD results correctly predict the relationship in Young's moduli. Therefore, ReaxFF is confirmed to be a useful tool for studying the mechanical behavior of epoxies.

Dynamics of proteins aggregation. I. Universal scaling in unbounded media

It is well understood that in some cases proteins do not fold correctly and, depending on their environment, even properly-folded proteins change their conformation spontaneously, taking on a misfolded state that leads to protein aggregation and formation of large aggregates. An important factor that contributes to the aggregation is the interactions between the misfolded proteins. Depending on the aggregation environment, the aggregates may take on various shapes forming larger structures, such as protein plaques that are often toxic. Their deposition in tissues is a major contributing factor to many neuro-degenerative diseases, such as Alzheimer's, Parkinson's, amyotrophic lateral sclerosis, and prion. This paper represents the first part in a series devoted to molecular simulation of protein aggregation. We use the PRIME, a meso-scale model of proteins, together with extensive discontinuous molecular dynamics simulation to study the aggregation process in an unbounded fluid system, as the first step toward MD simulation of the same phenomenon in crowded cellular environments. Various properties of the aggregates have been computed, including dynamic evolution of aggregate-size distribution, mean aggregate size, number of peptides that contribute to the formation of β sheets, number of various types of hydrogen bonds formed in the system, radius of gyration of the aggregates, and the aggregates' diffusivity. We show that many of such quantities follow dynamic scaling, similar to those for aggregation of colloidal clusters. In particular, at long times the mean aggregate size S(t) grows with time as, S(t) ∼ t^z, where z is the dynamic exponent. To our knowledge, this is the first time that the qualitative similarity between aggregation of proteins and colloidal aggregates has been pointed out.

Atomic and global mechanical properties of systems described by the Stillinger-Weber potential

Analytical expressions for the stress and elasticity tensors of materials, in which the interactions are described by the Stillinger-Weber potential, are derived in the context of the stress fluctuation formalism. The derived formulas can be used both in Monte Carlo and molecular dynamics simulations. As an example of possible applications, they are employed to calculate the influence of the temperature and system size on the mechanical properties of crystalline cubic boron nitride. The system has been studied by molecular dynamics simulations. The computed mechanical properties are in good agreement with available experimental data and first principle calculations. In the studied crystalline cubic boron nitride system, the employed formalism is of higher accuracy than the 'small-strain' non-equilibrium method. The dominant contributions to the elastic constants stem from the Born and stress fluctuation terms. An increase in the system size reduces the statistical uncertainties in the computation of the mechanical properties. A rise of the temperature leads to a slight increase in the observed uncertainties. The derived expressions for the stress and elasticity tensors are further decomposed into sums of atomic level stress and atomic level elasticity tensors. The developed factorization enables us (i) to quantify the contribution of the various chemical groups, in the case under consideration of the different atoms, to the observed mechanical properties and (ii) to determine the elastic constants with reduced computational uncertainties. The reason is that the exact values of some terms of the proposed factorization can be determined theoretically beforehand. Thus, they can be substituted in the derived formulas leading to an enhanced convergence.

First-Principles Monte Carlo Simulations of Reaction Equilibria in Compressed Vapors

Predictive modeling of reaction equilibria presents one of the grand challenges in the field of molecular simulation. Difficulties in the study of such systems arise from the need (i) to accurately model both strong, short-ranged interactions leading to the formation of chemical bonds and weak interactions arising from the environment, and (ii) to sample the range of time scales involving frequent molecular collisions, slow diffusion, and infrequent reactive events. Here we present a novel reactive first-principles Monte Carlo (RxFPMC) approach that allows for investigation of reaction equilibria without the need to prespecify a set of chemical reactions and their ideal-gas equilibrium constants. We apply RxFPMC to investigate a nitrogen/oxygen mixture at T = 3000 K and p = 30 GPa, i.e., conditions that are present in atmospheric lightning strikes and explosions. The RxFPMC simulations show that the solvation environment leads to a significantly enhanced NO concentration that reaches a maximum when oxygen is present in slight excess. In addition, the RxFPMC simulations indicate the formation of NO2 and N2O in mole fractions approaching 1%, whereas N3 and O3 are not observed. The equilibrium distributions obtained from the RxFPMC simulations agree well with those from a thermochemical computer code parametrized to experimental data.

Reactive symbol sequences for a model of hydrogen combustion

Transient, macroscopic states of chemical disequilibrium are born out of the microscopic dynamics of molecules. As a reaction mixture evolves, the temporal patterns of chemical species encodes some of this dynamical information, while their statistics are a manifestation of the bulk kinetics. Here, we define a chemically-informed symbolic dynamics as a coarse-grained representation of classical molecular dynamics, and analyze the sequences of chemical species for a model of hydrogen combustion. We use reactive molecular dynamics simulations to generate the sequences and derive probability distributions for sequence observables: the reaction time scales and the chain length - the total number of reactions between initiation of a reactant and termination at products. The time scales and likelihood of the sequences depend strongly on the chain length, temperature, and density. Temperature suppresses the uncertainty in chain length for hydrogen sequences, but enhances the uncertainty in oxygen sequence chain lengths. This method of analyzing a surrogate chemical symbolic dynamics reduces the complexity of the chemistry from the atomistic to the molecular level and has the potential for extension to more complicated reaction systems.

Pointillist rendering of electron charge and spin density suffices to replicate trends in atomic properties

The monotonic and non-monotonic variations of atomic properties within and between the rows of the periodic table underlie our understanding of chemistry and accounting for these variations has been a signature strength of quantum mechanics (QM). However, the computational burden of QM motivates the development of more efficient means of describing electrons and reactivity. The recently developed LEWIS • model incorporates lessons learnt from QM into a force field that includes electrons as explicit pseudo-classical particles. Here, we extend LEWIS • across the 2 p and 3 p elements, and show that it is capable of reproducing both monotonic and non-monotonic variations of chemically important atomic properties in a cost-effective manner. An indicator of the strength of the construct is the ability of pairwise potentials trained on ionization energies and the order of spin configurations to predict atomic polarizabilities. In this manner, some insights of QM are uncoupled from its onerous computational burden.

Reactive force fields made simple

Combining QMDFF with EVB allows to generate reactive force fields of useful quality for widely varying reactions with minimal effort.

Multisurface Adiabatic Reactive Molecular Dynamics

Adiabatic reactive molecular dynamics (ARMD) simulation method is a surface-crossing algorithm for modeling chemical reactions in classical molecular dynamics simulations using empirical force fields. As the ARMD Hamiltonian is time dependent during crossing, it allows only approximate energy conservation. In the current work, the range of applicability of conventional ARMD is explored, and a new multisurface ARMD (MS-ARMD) method is presented, implemented in CHARMM and applied to the vibrationally induced photodissociation of sulfuric acid (H2SO4) in the gas phase. For this, an accurate global potential energy surface (PES) involving 12 H2SO4 and 4 H2O + SO3 force fields fitted to MP2/6-311G++(2d,2p) reference energies is employed. The MS-ARMD simulations conserve total energy and feature both intramolecular H-transfer reactions and water elimination. An analytical treatment of the dynamics in the crossing region finds that conventional ARMD can approximately conserve total energy for limiting cases. In one of them, the reduced mass of the system is large, which often occurs for simulations of solvated biomolecular systems. On the other hand, MS-ARMD is a general approach for modeling chemical reactions including gas-phase, homogeneous, heterogeneous, and enzymatic catalytic reactions while conserving total energy in atomistic simulations.

ReaxFF molecular dynamics simulations of CO collisions on an O-preadsorbed silica surface

A quasiclassical trajectory dynamics study was performed for carbon monoxide collisions over an oxygen preadsorbed β-cristobalite (001) surface. A reactive molecular force field (ReaxFF) was used to model the potential energy surface. The collisions were performed fixing several initial conditions: CO rovibrational states (v = 0-5 and j = 0, 20, 35), collision energies (0.05 ≤ E(col) ≤ 2.5 eV), incident angles (θ(v) = 0°, 45°) and surface temperatures (T(surf) = 300 K, 900 K). The principal elementary processes were the molecular reflection and the non-dissociative molecular adsorption. CO₂ molecules were also formed in minor extension via an Eley-Rideal reaction although some of them were finally retained on the surface. The scattered CO molecules tend to be translationally colder and internally hotter (rotationally and vibrationally) than the initial ones. The present study supports that CO + O(ad) reaction should be less important than O + O(ad) reaction over silica for similar initial conditions of reactants, in agreement with experimental data.