Quantum mechanics based non-bonded force field functions for use in molecular dynamics simulations of materials and systems: The nitrogen and oxygen columns.

Accurate Force Fields (FFs) are essential for Molecular Dynamics (MD) simulations of the dynamics of realistic materials in terms of atomic-level interactions. The FF parameters of short-range valence interactions can be derived through Quantum Mechanical (QM) calculations on model systems practical for QM (<300 atoms). Similarly, the dynamic electrostatic interactions can be described with methods such as QEq or PQEq that allow charges and polarization to adjust dynamically. However, accurately extracting long-range van der Waals (vdW) interactions from QM calculations poses challenges due to the absence of a definitive method to distinguish between the different energetic components of electrostatics, polarization, vdW, hydrogen bonding, and valence interactions. To do this we use the Perdew-Burke-Ernzerhof flavor of Density Functional Theory, including empirical D3 vdW corrections, to predict the Equation of State for each element (keeping any covalent bonds fixed), from which we obtain the two-body vdW nonbond potential. Here, we extend these calculations to include non-bonded parameters for the N and O columns of the periodic table so that we now describe columns 15 (N), 16 (O), 17 (F), and 18 (Ne) of the periodic table. For these 20 elements, we find that the two-body vdW potentials can all be mapped to a single universal two-body curve, with just three scaling parameters: Re, De, and L. We refer to this as the Universal NonBond (UNB) potential. We expect this to be useful for new MD simulations and a helpful starting point to obtain UNB parameters for the remainder of the periodic table.

Topology induced crossover between Langevin, subdiffusion, and Brownian diffusion regimes in supercooled water.

Despite extensive studies of supercooled water, it remains challenging to understand its peculiar dynamic anomalous properties. In this work, we integrated full atomistic simulations of supercooled water over the temperature range of room temperature to 200 K using a quantum-mechanics-based polarizable force field with the dressed dynamics method that couples fast collision events and slow reorganization dynamics of hydrogen-bond networks. Our analysis unveils the salient multiscale features in the transient relaxation dynamics of supercooled water. A classical Langevin behavior dominates at fast timescales, while long-time relaxations unveil two different activation barriers in two temperature regions: below and above 230 K. The modulation of the entropy spectrum by temperature is elucidated in terms of a three-state model underlined by the complexity of the water dynamics associated with a topological transition of a strong hydrogen-bond network. This state-dependent network topology is quantitatively characterized by power-law exponents of inverse network connectivity from 200 to 298 K. This work provides valuable guidance for further studies on the transient relaxation dynamics of supercooled water.

Li-diffusion at the interface between Li-metal and [Pyr14][TFSI]-ionic liquid: Ab initio molecular dynamics simulations.

We previously reported comprehensive density functional theory-molecular dynamics (DFT-MD) at 400 K to determine the composition and structure of the solid electrolyte interface (SEI) between a Li anode and [Pyr14][TFSI] ionic liquid. In this paper, we examined diffusion rates in both the Li-electrode region and SEI compact layer in smaller 83Li/2[TFSI] and larger 164Li/4[TFSI] systems. At 400 K, the Li-diffusion constant in the Li-region is 1.35 × 10-10 m2/s for 83Li/2[TFSI] and 5.64 × 10-10 m2/s for 164Li/4[TFSI], while for the SEI it is 0.33 × 10-10 m2/s and 0.22 × 10-10 m2/s, thus about one order slower in the SEI compared to the Li-region. This Li-diffusion is dominated by hopping from the neighbor shell of one F or O to the neighbor shell of another. Comparing the Li-diffusion at different temperatures, we find that the activation energy is 0.03 and 0.11 eV for the Li-region in the smaller and larger systems, respectively, while for the SEI it is 0.09 and 0.06 eV.

Accurate non-bonded potentials based on periodic quantum mechanics calculations for use in molecular simulations of materials and systems.

Molecular dynamics simulations require accurate force fields (FFs) to describe the physical and chemical properties of complex materials and systems. FF parameters for valence interactions can be determined from high-quality Quantum Mechanical (QM) calculations. However, it has been challenging to extract long-range nonbonded interaction potentials from QM calculations since there is no unambiguous method to separate the total QM energy into electrostatics (polarization), van der Waals (vdW), and other components. Here, we propose to use density functional theory with dispersion corrections to obtain the equation of state for single element solid systems (of H, C, N, O, F, Cl, Br, I, P, He, Ne, Ar, Kr, Xe, and Rn) from which we obtain the pure 2-body vdW nonbonded potentials. Recently, we developed the polarizable charge equilibration (PQEq) model based on QM polarization energy of electric probe dipoles with no contributions from vdW. Together, the vdW and PQEq interactions form the nonbonded potential of our new transferrable reactive FF (RexPoN). They may also be useful to replace the nonbonded parts of standard FFs, such as OPLS, Amber, UFF, and CHARMM. We find that the individual 2-body vdW potential curves can be scaled to a universal vdW potential using just three specific atomic parameters. This simplifies extension to the rest of the periodic table for atoms that do not exhibit molecular packing. We validate the accuracy of these nonbonded interactions for liquid water, energetic, and biological systems. In all cases, we find that our new nonbonded potentials provide good agreement with QM and experimental data.

Anomalies in Supercooled Water at ∼230 K Arise from a 1D Polymer to 2D Network Topological Transformation.

Puzzling anomalous properties of water are drastically enhanced in the supercooled region. However, the nature of these anomalies is not known. We report here molecular dynamics simulations using the RexPoN force field from 298 to 200 K along the 1 atm density curve. At 298 K, there are 2.1 strong hydrogen bonds (SHBs), leading to a dynamic branched one-dimensional (1D) polymer. Water remains 1D down to 240 K, but at and below 230 K, the number of SHBs becomes 3.0, leading to a two-dimensional (2D) network that persists to 200 K. We propose that this 1D-to-2D topological transition accounts for the anomalous properties of supercooled water. Near 230 K, the power spectra show dramatic increases in the angular vibrational frequency modes, while the diffusivity decreases dramatically, both arising from the 1D-to-2D transformation. This transition is not first order because free energy changes uniformly but fluctuations in the entropy near 230 K suggest a possible second-order transition.

Interface Structure in Li-Metal/[Pyr14][TFSI]-Ionic Liquid System from ab Initio Molecular Dynamics Simulations.

Ionic liquids (ILs) are promising materials for application in a new generation of Li batteries. They can be used as electrolyte or interlayer or incorporated into other materials. ILs have the ability to form a stable solid electrochemical interface (SEI), which plays an important role in protecting the Li-based electrode from oxidation and the electrolyte from extensive decomposition. Experimentally, it is hardly possible to elicit fine details of the SEI structure. To remedy this situation, we have performed a comprehensive computational study (density functional theory-based molecular dynamics) to determine the composition and structure of the SEI compact layer formed between the Li anode and [Pyr14][TFSI] IL. We found that the [TFSI] anions quickly reacted with Li and decomposed, unlike the [Pyr14] cations which remained stable. The obtained SEI compact layer structure is nonhomogeneous and consists of the atomized S, N, O, F, and C anions oxidized by Li atoms.

Liquid water is a dynamic polydisperse branched polymer.

We developed the RexPoN force field for water based entirely on quantum mechanics. It predicts the properties of water extremely accurately, with Tmelt = 273.3 K (273.15 K) and properties at 298 K: ΔHvap = 10.36 kcal/mol (10.52), density = 0.9965 g/cm3 (0.9965), entropy = 68.4 J/mol/K (69.9), and dielectric constant = 76.1 (78.4), where experimental values are in parentheses. Upon heating from 0.0 K (ice) to 273.0 K (still ice), the average number of strong hydrogen bonds (SHBs, rOO ≤ 2.93 Å) decreases from 4.0 to 3.3, but upon melting at 273.5 K, the number of SHBs drops suddenly to 2.3, decreasing slowly to 2.1 at 298 K and 1.6 at 400 K. The lifetime of the SHBs is 90.3 fs at 298 K, increasing monotonically for lower temperature. These SHBs connect to form multibranched polymer chains (151 H2O per chain at 298 K), where branch points have 3 SHBs and termination points have 1 SHB. This dynamic fluctuating branched polymer view of water provides a dramatically modified paradigm for understanding the properties of water. It may explain the 20-nm angular correlation lengths at 298 K and the critical point at 227 K in supercooled water. Indeed, the 15% jump in the SHB lifetime at 227 K suggests that the supercooled critical point may correspond to a phase transition temperature of the dynamic polymer structure. This paradigm for water could have a significant impact on the properties for protein, DNA, and other materials in aqueous media.

First principles-based multiscale atomistic methods for input into first principles nonequilibrium transport across interfaces.

This issue of PNAS features "nonequilibrium transport and mixing across interfaces," with several papers describing the nonequilibrium coupling of transport at interfaces, including mesoscopic and macroscopic dynamics in fluids, plasma, and other materials over scales from microscale to celestial. Most such descriptions describe the materials in terms of the density and equations of state rather than specific atomic structures and chemical processes. It is at interfacial boundaries where such atomistic information is most relevant. However, there is not yet a practical way to couple these phenomena with the atomistic description of chemistry. The starting point for including such information is the quantum mechanics (QM). However, practical QM calculations are limited to a hundred atoms for dozens of picoseconds, far from the scales required to inform the continuum level with the proper atomistic description. To bridge this enormous gap, we need to develop practical methods to extend the scale of the atomistic simulation by several orders of magnitude while retaining the level of QM accuracy in describing the chemical process. These developments would enable continuum modeling of turbulent transport at interfaces to incorporate the relevant chemistry. In this perspective, we will focus on recent progress in accomplishing these extensions in first principles-based atomistic simulations and the strategies being pursued to increase the accuracy of very large scales while dramatically decreasing the computational effort.

Polarizable Charge Equilibration Model for Transition-Metal Elements.

The polarizable charge equilibration (PQEq) method was developed to provide a simple but accurate description of the electrostatic interactions and polarization effects in materials. Previously, we optimized four parameters per element for the main group elements. Here, we extend this optimization to the 24 d-block transition-metal (TM) elements, columns 4-11 of the periodic table including Ti-Cu, Zr-Ag, and Hf-Au. We validate the PQEq description for these elements by comparing to interaction energies computed by quantum mechanics (QM). Because many materials applications involving TM are for oxides and other compounds that formally oxidize the metal, we consider a variety of oxidation states in 24 different molecular clusters. In each case, we compare interaction energies and induced fields from QM and PQEq along various directions. We find that the original χ and J parameters (electronegativity and hardness) related to the ionization of the atom remain valid; however, we find that the atomic radius parameter needs to be close to the experimental ionic radii of the transition metals. This leads to a much higher spring constant to describe the atomic polarizability. We find that these optimized parameters for PQEq provide accurate interaction energies compared to QM with charge distributions that depend in a reasonable way on the coordination number and oxidation states of the transition metals. We expect that this description of the electrostatic interactions for TM will be useful in molecular dynamics simulations of inorganic and organometallic materials.

The quantum mechanics-based polarizable force field for water simulations.

We report here a new force field for water based solely on quantum mechanics (QM) calculations with no empirical data. The QM was at a high level, coupled cluster single double triple, for all orientations and distances for water dimer plus X3LYP density functional theory (DFT) on 19 larger water clusters. In addition, we included charge and polarization based on the polarizable charge equilibration method and nonbond interactions from DFT-D3 calculations on the H2 and O2 crystal. This model, denoted as RexPoN, provides quite excellent agreement with experimental (expr) data for the solid and liquid phase of water: T melt = 273.3 K (expr = 273.15 K) and properties at 298 K: ΔH vap = 10.36 kcal/mol (expr = 10.52), density = 0.9965 gr/cm3 (expr = 0.9965), entropy = 68.4 (J/mol)/K (expr = 69.9), dielectric constant = 76.1 (expr = 78.4), and ln D s (self-diffusion coef) = -10.08 (expr = -11.24). Such an accurate force field for water will, we believe, be useful for full solvent calculations of electrocatalysis, where we can restrict QM water to just the first one or two layers involving reactions, using RexPoN to provide the polarization for a more distant solvent. Also, RexPoN may provide a better description of the solvent for proteins, DNA, polymers, and inorganic systems for applications to biomolecular, pharma, electrocatalysis (fuel cells and water splitting), and batteries where interaction with explicit water molecules plays a significant role.

Predicted detonation properties at the Chapman-Jouguet state for proposed energetic materials (MTO and MTO3N) from combined ReaxFF and quantum mechanics reactive dynamics.

The development of new energetic materials (EMs) with improved detonation performance but low sensitivity and environmental impact is of considerable importance for applications in civilian and military fields. Often new designs are difficult to synthesize so predictions of performance in advance is most valuable. Examples include MTO (2,4,6-triamino-1,3,5-triazine-1,3,5-trioxide) and MTO3N (2,4,6-trinitro-1,3,5-triazine-1,3,5-trioxide) suggested by Klapötke as candidate EMs but not yet successfully synthesized. We propose and apply to these materials a new approach, RxMD(cQM), in which ReaxFF Reactive Molecular Dynamics (RxMD) is first used to predict the reaction products and thermochemical properties at the Chapman Jouguet (CJ) state for which the system is fully reacted and at chemical equilibrium. Quantum mechanics dynamics (QMD) is then applied to refine the pressure of the ReaxFF predicted CJ state to predict a more accurate final CJ point, leading to a very practical calculation that includes accurate long range vdW interactions needed for accurate pressure. For MTO, this RxMD(cQM) method predicts a detonation pressure of PCJ = 40.5 GPa and a detonation velocity of DCJ = 8.8 km s-1, while for MTO3N it predicts PCJ = 39.9 GPa and DCJ = 8.4 km s-1, making them comparable to HMX (PCJ = 39.5 GPa, DCJ = 9.1 km s-1) and worth synthesizing. This first-principles-based RxMD(cQM) methodology provides an excellent compromise between computational cost and accuracy including the formation of clusters that burn too slowly, providing a practical mean of assessing detonation performances for novel candidate EMs. This RxMD(cQM) method that links first principles atomistic molecular dynamics simulations with macroscopic properties to promote in silico design of new EMs should also be of general applicability to materials synthesis and processing.

Extension of the Polarizable Charge Equilibration Model to Higher Oxidation States with Applications to Ge, As, Se, Br, Sn, Sb, Te, I, Pb, Bi, Po, and At Elements.

We recently developed the polarizable charge equilibration (PQEq) model to predict accurate electrostatic interactions for molecules and solids and optimized parameters for H, C, N, O, F, Si, P, S, and Cl elements to fit polarization energies computed by quantum mechanics (QM). Here, we validate and optimize the PQEq parameters for other p-block elements including Ge, As, Se, Br, Sn, Sb, Te, I, Pb, Bi, Po, and At using 28 molecular structures containing these elements. For these elements, we now include molecules with higher oxidation states: III and V for the As column, IV and VI for the Se column, and I, III, and V for the Br column. We find that PQEq predicts polarization energies in excellent agreement with QM.

Polarizable charge equilibration model for predicting accurate electrostatic interactions in molecules and solids.

Electrostatic interactions play a critical role in determining the properties, structures, and dynamics of chemical, biochemical, and material systems. These interactions are described well at the level of quantum mechanics (QM) but not so well for the various models used in force field simulations of these systems. We propose and validate a new general methodology, denoted PQEq, to predict rapidly and dynamically the atomic charges and polarization underlying the electrostatic interactions. Here the polarization is described using an atomic sized Gaussian shaped electron density that can polarize away from the core in response to internal and external electric fields, while at the same time adjusting the charge on each core (described as a Gaussian function) so as to achieve a constant chemical potential across all atoms of the system. The parameters for PQEq are derived from experimental atomic properties of all elements up to Nobelium (atomic no. = 102). We validate PQEq by comparing to QM interaction energy as probe dipoles are brought along various directions up to 30 molecules containing H, C, N, O, F, Si, P, S, and Cl atoms. We find that PQEq predicts interaction energies in excellent agreement with QM, much better than other common charge models such as obtained from QM using Mulliken or ESP charges and those from standard force fields (OPLS and AMBER). Since PQEq increases the accuracy of electrostatic interactions and the response to external electric fields, we expect that PQEq will be useful for a large range of applications including ligand docking to proteins, catalytic reactions, electrocatalysis, ferroelectrics, and growth of ceramics and films, where it could be incorporated into standard force fields as OPLS, AMBER, CHARMM, Dreiding, ReaxFF, and UFF.

Molecular Simulation Study of Gas Solubility and Diffusion in a Polymer-Boron Nitride Nanotube Composite.

We study the possibility of using polymer composites made of a polymer and boron nitride nanotubes (BNNTs) as a new type of membranes for gas separation. The polymer used is amorphous poly(ether imide) (PEI), and zigzag BNNTs are used to generate the composites with the PEI. The solubilities and self-diffusivities of CO2 and CH4 in the PEI and its composites with the BNNTs are calculated by molecular dynamics (MD) simulations. The molecular models of the PEI and its composites with the BNNTs are generated using energy minimization and MD simulation, and the Universal Force Field is used to represent the interactions between all the atoms. The morhology of the composites are characterized and are compared with that of PEI. The accuracy of the computations is tested by calculating the gases' solubilities and self-diffsivities in the pure PEI and comparing them with the experimental data. Good agreement is obtained with the data. The computed diffusivities and solubilities in the polymer-BNNTs composites are much larger than those in the pure polymer, which are attributed to the changes that the BNNTs induce in the polymer composite's free-volume distribution. As the mechanical properties of the polymer-BNNTs composites are superior over those of the pure PEI, their use as a membrane for gas separation offers distinct advantages over the pure polymer. We also demonstrate that, calculating the diffusion coefficients with MD simulations in the NPT ensemble, as opposed to the common practice of utilizing the NVT ensemble, leads to much more accurate results.

First principles-based multiparadigm, multiscale strategy for simulating complex materials processes with applications to amorphous SiC films.

Progress has recently been made in developing reactive force fields to describe chemical reactions in systems too large for quantum mechanical (QM) methods. In particular, ReaxFF, a force field with parameters that are obtained solely from fitting QM reaction data, has been used to predict structures and properties of many materials. Important applications require, however, determination of the final structures produced by such complex processes as chemical vapor deposition, atomic layer deposition, and formation of ceramic films by pyrolysis of polymers. This requires the force field to properly describe the formation of other products of the process, in addition to yielding the final structure of the material. We describe a strategy for accomplishing this and present an example of its use for forming amorphous SiC films that have a wide variety of applications. Extensive reactive molecular dynamics (MD) simulations have been carried out to simulate the pyrolysis of hydridopolycarbosilane. The reaction products all agree with the experimental data. After removing the reaction products, the system is cooled down to room temperature at which it produces amorphous SiC film, for which the computed radial distribution function, x-ray diffraction pattern, and the equation of state describing the three main SiC polytypes agree with the data and with the QM calculations. Extensive MD simulations have also been carried out to compute other structural properties, as well the effective diffusivities of light gases in the amorphous SiC film.

General Multiobjective Force Field Optimization Framework, with Application to Reactive Force Fields for Silicon Carbide.

First-principles-based force fields prepared from large quantum mechanical data sets are now the norm in predictive molecular dynamics simulations for complex chemical processes, as opposed to force fields fitted solely from phenomenological data. In principle, the former allow improved accuracy and transferability over a wider range of molecular compositions, interactions, and environmental conditions unexplored by experiments. That is, assuming they have been optimally prepared from a diverse training set. The trade-off has been force field engines that are functionally complex, with a large number of nonbonded and bonded analytical forms that give rise to rather large parameter search spaces. To address this problem, we have developed GARFfield (genetic algorithm-based reactive force field optimizer method), a hybrid multiobjective Pareto-optimal parameter development scheme based on genetic algorithms, hill-climbing routines and conjugate-gradient minimization. To demonstrate the capabilities of GARFfield we use it to develop two very different force fields: (1) the ReaxFF reactive force field for modeling the adiabatic reactive dynamics of silicon carbide growth from an methyltrichlorosilane precursor and (2) the SiC electron force field with effective core pseudopotentials for modeling nonadiabatic dynamic phenomena with highly excited electronic states. The flexible and open architecture of GARFfield enables efficient and fast parallel optimization of parameters from quantum mechanical data sets for demanding applications like ReaxFF, electronic fast forward (or electron force field), and others including atomistic reactive charge-optimized many-body interatomic potentials, Morse, and coarse-grain force fields.