Source: https://www.nature.com/articles/npjcompumats201511?error=cookies_not_supported&code=ec75438b-63c4-4755-963c-ecf13ecd4b6a
Timestamp: 2019-04-21 09:41:06+00:00

Document:
The reactive force-field (ReaxFF) interatomic potential is a powerful computational tool for exploring, developing and optimizing material properties. Methods based on the principles of quantum mechanics (QM), while offering valuable theoretical guidance at the electronic level, are often too computationally intense for simulations that consider the full dynamic evolution of a system. Alternatively, empirical interatomic potentials that are based on classical principles require significantly fewer computational resources, which enables simulations to better describe dynamic processes over longer timeframes and on larger scales. Such methods, however, typically require a predefined connectivity between atoms, precluding simulations that involve reactive events. The ReaxFF method was developed to help bridge this gap. Approaching the gap from the classical side, ReaxFF casts the empirical interatomic potential within a bond-order formalism, thus implicitly describing chemical bonding without expensive QM calculations. This article provides an overview of the development, application, and future directions of the ReaxFF method.
Atomistic-scale computational techniques provide a powerful means for exploring, developing and optimizing promising properties of novel materials. Simulation methods based on quantum mechanics (QM) have grown in popularity over recent decades due to the development of user-friendly software packages making QM level calculations widely accessible. Such availability has proved particularly relevant to material design, where QM frequently serves as a theoretical guide and screening tool. Unfortunately, the computational cost inherent to QM level calculations severely limits simulation scales. This limitation often excludes QM methods from considering the dynamic evolution of a system, thus hampering our theoretical understanding of key factors affecting the overall behaviour of a material. To alleviate this issue, QM structure and energy data are used to train empirical force fields that require significantly fewer computational resources, thereby enabling simulations to better describe dynamic processes. Such empirical methods, including reactive force-field (ReaxFF),1 trade accuracy for lower computational expense, making it possible to reach simulation scales that are orders of magnitude beyond what is tractable for QM.
Atomistic force-field methods utilise empirically determined interatomic potentials to calculate system energy as a function of atomic positions. Classical approximations are well suited for nonreactive interactions, such as angle-strain represented by harmonic potentials, dispersion represented by van der Waals potentials and Coulombic interactions represented by various polarisation schemes. However, such descriptions are inadequate for modelling changes in atom connectivity (i.e., for modelling chemical reactions as bonds break and form). This motivates the inclusion of connection-dependent terms in the force-field description, yielding a reactive force-field. In ReaxFF, the interatomic potential describes reactive events through a bond-order formalism, where bond order is empirically calculated from interatomic distances. Electronic interactions driving chemical bonding are treated implicitly, allowing the method to simulate reaction chemistry without explicit QM consideration.
The classical treatment of reactive chemistry made available by the ReaxFF methodology has opened the door for numerous studies of phenomena occurring on scales that were previously inaccessible to computational methods. In particular, ReaxFF enables simulations involving reactive events at the interface between solid, liquid, and gas phases, which is made possible because the ReaxFF description of each element is transferable across phases. For example, an oxygen atom is treated with the same mathematical formalism whether that oxygen is in the gas phase as O2, in the liquid phase within an H2O molecule, or incorporated in a solid oxide. Such transferability, coupled with a lower computational expense allowing for longer simulation timescales, allows ReaxFF to consider phenomena dependent not only on the reactivity of the involved species, but also on dynamic factors, such as diffusivity and solubility, affecting how species migrate through the system. This allows ReaxFF to model complex processes involving multiple phases in contact with one another.
To demonstrate these capabilities, this article will review the development and application of the ReaxFF method. In ‘History of ReaxFF development’, we discuss initial development choices shaping the overall ReaxFF formalism, whereas ‘Current ReaxFF methodology’ and ‘Overview of available ReaxFF parameterisations and development branches’ outline the currently employed method and available parameter sets, respectively. Comparisons with other reactive potentials available in the literature are briefly discussed in ‘Comparison to similar methods’. Examples of various ReaxFF applications are provided in ‘Applications of ReaxFF’, with emphasis placed on demonstrating the breadth of systems that can be modelled with the method. Finally, plans for future extensions and improvements are discussed in ‘Future developments and outlook’.
The current functional form of the ReaxFF potential, best described in the Chenoweth et al.2 hydrocarbon combustion work (herein referred to as 2008-C/H/O),2 has demonstrated significant transferability across the periodic table. It is important to note, however, that the 2008 functional form is different from the original 2001 ReaxFF hydrocarbon description,1 as well as from the 2003 extension to silicon and silica.3 Although conceptually similar to the current 2008-C/H/O functional form, the 2001 hydrocarbon description employed the same dissociation energy for C–C single, double and triple bonds. This approach was reasonable for hydrocarbons, but could not be extended to treat Si–O single and double bonds. As such, the 2003 Si/O/H extension required separate parameters describing single-, double- and triple-bond dissociation. Furthermore a lone-pair energy term was introduced to handle formation and dissociation of oxygen lone-pairs. The 2003 Si/O functional form was further augmented by a three-body conjugation term introduced to handle −NO2 group chemistry in nitramines, where a triple-bond stabilisation term was added to improve the description of terminal triple bonds. This led to the 2003–2005 ReaxFF description for the RDX high-energy material employed by Strachan et al.4,5 to study RDX initiation.
Since 2005, the ReaxFF functional form has been stable, although optional additions, such as angular terms to destabilise Mg–Mg–H zero-degree angles6 or double-well angular terms necessary for describing aqueous transition metal ions,7 have occasionally been added to the potential. Goddard and co-workers implemented an additional attractive van der Waals term to improve performance for nitramine crystals (ReaxFF-lg).8 This concept, however, was not made transferable with previous or later ReaxFF parameter sets. The 2005 functional form developed for RDX is the current version of ReaxFF distributed by the van Duin group (commonly referred to as ‘standalone ReaxFF’), as well as integrated in the open-source LAMMPS code,9 supported through Nanohub, (http://www.nanohub.org) and available through the PuReMD (Purdue Reactive Molecular Dynamics) code.10,​11,​12 Apart from these open-source distributions, the ReaxFF method is also integrated in ADF13 (released by SCM (http://www.scm.com)) and in Materials Studio (released under license by Accelrys (http://www.accelrys.com)). The pre-2005 ReaxFF parameter sets, including the aforementioned 2001-C/H,1 2003-Si/O,3 and 2004-Al/O14 descriptions, are not supported by any codes curated by the van Duin group or its collaborators. The materials described by the three pre-2005 ReaxFF parameterisations are equally, if not better, described by later parameterisations. The 2008-C/H/O parameter set was trained against the entire 2001-C/H training set, while the 2010 and 2011 Si/O/H parameterisations (Fogarty et al.15 on the ‘aqueous branch’ and Neyts et al.16 on the ‘combustion branch’) were validated against the full 2003-Si/O/H training set. Finally, the 2008-Al/H ReaxFF description, and later applications to aluminium oxides and aluminosilicates,17,​18,​19,​20,​21,​22,​23 fully contain and extend the 2004-Al/O description. We are aware of other ReaxFF implementations, often developed by individual scientists based on the 2008-C/H/O formalism. Given the complexity of the ReaxFF functional form, it is advisable to validate ReaxFF implementations against the standalone ReaxFF code prior to applying them in production-scale simulations.
(a) Overview of the ReaxFF total energy components and (b) elements currently described in available parameter sets.
As depicted in Figure 1a, the potential is divided into bond-order-dependent and -independent contributions. Bond order is calculated directly from interatomic distance using the empirical formula: (2)BOij=BOijσ+BOijπ+BOijππ=exp[pbo1(rijroσ)pbo2]+exp[pbo3(rijroπ)pbo4]+exp[pbo5(rijroππ)pbo6] where BO is the bond order between atoms i and j, rij is interatomic distance, ro terms are equilibrium bond lengths, and pbo terms are empirical parameters. Equation (2) is continuous, containing no discontinuities through transitions between σ, π, and ππ bond character. This yields a differentiable potential energy surface, as required for the calculation of interatomic forces. This bond-order formula accommodates long-distance covalent interactions characteristic in transition state structures, allowing the force-field to accurately predict reaction barriers. This covalent range is typically taken to be 5 Angstrom—which is sufficient for most elements to capture even the weakest of covalent interactions—but can be extended beyond this range; this may occasionally be required for elements with very large covalent radii. This long-distance covalent bond feature, however, necessitates the addition of a bond-order correction to remove spurious bond character between non-bonded neighbours, such as neighbouring H atoms in a methane molecule. Terms in the potential that are dependent on bond order, such as bond energy and angle strains, are calculated directly from the corrected bond order. Finally, a charge equilibration scheme is applied at each iteration to calculate partial atomic charges (see ‘Charge description improvements’), which are then used to calculate Coulombic interactions.
Note that the non-bonded and bonded terms in ReaxFF are calculated independently—there is no information transfer between the bond-order-based terms and the van der Waals- and Coulomb-related terms. For all materials and molecules, both bond-order-based terms and nonbonded terms are calculated—without exclusions—enabling ReaxFF to be applied to both predominantly covalent and ionic materials without user input.
Each element of the periodic table for which a ReaxFF parameter set has been published is highlighted in Figure 1b. Although a ReaxFF description exists for these elements, one cannot simply use these parameter sets in any combination and expect to obtain satisfactory transferability. As shown in Figure 2, there are currently two major groupings (i.e., the ReaxFF branches) of parameter sets that are intra-transferable with one another: (1) the combustion branch and (2) the aqueous branch. In order to explain the existence of branches within a common functional form, we must consider that ReaxFF does not employ atom typing strategies, in contrast to popular non-reactive force-fields like AMBER25 and CHARMM.26 For example, there is only one oxygen type in ReaxFF regardless of the chemical environment in which the oxygen atom finds itself. This is quite helpful because it allows atoms to migrate seamlessly between phases during a simulation. At the same time, however, this results in a significantly more complex force-field development process. The lack of transferability between branches is evident in the performance of the 2008-C/H/O combustion force-field,2 which accurately describes water as a gas-phase molecule yet fails to describe water as a liquid. During development at that time, describing liquid water was not a particular aim for ReaxFF. Since all intended applications were at temperatures well above the water boiling point, this was not a major development concern. In 2009, efforts were initiated to redevelop ReaxFF for aqueous chemistry, and it became clear that the 2008-C/H/O combustion force field, which at that time was already extended to a significant range of metal oxide (Me=V/Bi/Mo/Nb/Si) materials and catalysts,27,​28,​29,​30,​31 could not be parameterised to treat liquid water without changing general-ReaxFF and atom-specific parameters. As such, the decision was made to initiate a new branch—the aqueous-branch—that employs the same functional form as the 2008-C/H/O description, but with different O/H atom and bond parameters. This resulted in the creation of a number of parameter sets that are not directly transferable with those on the pre-existing combustion branch, thus leading to the new aqueous branch (Figure 2).
ReaxFF development tree, where parameter sets on a common ‘branch’ are fully transferable with one another. Parameter sets are available in (a) van Duin et al.;1 (b) van Duin et al.;3 (c) Strachan et al.;4,5 (d) Chenoweth et al.;2 (e) Goddard et al.29 and Chenoweth et al.;28 (f) Castro-Marcano et al.163 and Kamat et al.;13 (g) Vasenkov et al.;164 (h) Weismiller et al.;32 (i) Mueller et al.;33 (j) Agrawalla et al.;165 (k) Rahaman et al.;116 (l) Monti et al.;53 (m) Fogarty et al.;15 (n) Raymand et al.;50 (o) Pitman et al.23 and Manzano et al.;166 and (p) van Duin et al.7 (Images adapted with permission from the provided references: (d, e, h–k, o, p) Copyright 2008, 2010, 2011 and 2012 American Chemical Society; (f, n) Copyright 2010 and 2012 Elsevier; (g, m) Copyright 2010 and 2012 AIP Publishing LLC; (l) Copyright 2013 Royal Society of Chemistry).
In general, empirical methods tend to be faster and scale better than QM-based approaches,79 although relatively few head-to-head comparisons have appeared in the literature so far. Within the empirical methods, ReaxFF’s origin in hydrocarbon chemistry has led to a focus on reproducing energy barriers, while methods developed by the materials community (e.g., EAM and COMB) tend to focus more on elastic properties. Long-range bond-order terms render ReaxFF suitable for transition states, but complicate bond-order descriptions in condensed systems when non-bonded neighbours may be in close proximity. This does not necessarily exclude the method from reproducing elastic properties (e.g., refs 80,81), although it makes these properties less straightforward to include in force-field parameterisation.
Here we highlight a few of the many studies that have employed the ReaxFF methodology, with the objective of demonstrating the diverse applicability of ReaxFF rather than providing a comprehensive review of the literature. In ‘Heterogeneous catalysis’ and ‘Atomic layer deposition’ we discuss the application of ReaxFF to heterogeneous catalysis and to atomic layer deposition (ALD), respectively, as these applications demonstrate the methodological strength of ReaxFF: modelling reactive chemistry at heterogeneous interfaces. A broader overview of other applications will be provided in ‘Other applications’. Finally, in Section 3.4 we discuss the development of computational methodologies that take advantage of the ReaxFF formalism to explore phenomena on time and length scales that are inaccessible to traditional MD simulation methodologies.
Ni-catalysed carbon nanotube (CNT) growth involves dissolution of C atoms into Ni nanoclusters, with CNT formation dependent on the dynamic reformation and transport of carbon through the cluster. Modelling such processes clearly requires some treatment of bond scission and formation, yet the time and length scales inherent to C diffusion in nm-sized clusters is not tractable with QM. For this reason, the ReaxFF force-field has been extensively applied to model catalytic processes involving CNT growth initiated from dissociative hydrocarbon adsorption on Ni surfaces. Mueller et al. developed82 and employed33 a Ni/C/H parameter set to investigate the onset of CNT formation via the dissociative adsorption of various hydrocarbons on Ni surfaces, where they found that surface defects likely play an essential role initiating CNT growth. ReaxFF enabled MD simulations to reach beyond the time and length scales available to ab initio methods, which proved essential in modelling defective Ni surfaces and clusters exposing active sites on irregular surface terminations and cluster edges. The Ni/C/H parameter set was further utilised by Neyts et al.83 to model the effect of Ar+ bombardment on CNT formation. Simulations showed that the energy of impinging Ar+ ions can be tuned to break weak C–C bonds between under-coordinated C atoms, thus healing CNT defects. In similar works,84,​85,​86 a force-biased Monte Carlo (MC) method, which will be further discussed in ‘Uniform-acceptance force-biased Monte Carlo’, was employed to demonstrate the growth of CNTs with well-defined chirality and the effect of applied electric fields on plasma-assisted CNT growth.
In addition to catalysis by metal surfaces, ReaxFF is also an effective methodology for modelling the catalytic properties of oxides. Among the first studies in this area, Chenoweth et al. developed and implemented a V/O/C/H description aimed at capturing the interaction of hydrocarbons with vanadium oxide.27,28 Mixed-metal-oxide (MMO) catalysts typically feature partial, mixed, and irregular metal occupations at various crystallographic sites. This characteristic is difficult to capture with QM models, as prohibitively large unit cells are required if the system is to be treated with periodic boundary conditions. In ref. 28, a combined MC/reactive dynamics (MC/RD) procedure was employed to explore possible MMO structures, in which the MC/RD scheme was used to systematically interchange metal atoms between crystallographic sites to determine the lowest energy structures. Thermodynamically favoured structures were then exposed to hydrocarbons through MD simulations, demonstrating that optimal catalysts contain several surface channels exposing active metal sites. This study demonstrates the utility of ReaxFF when modelling catalytic processes on oxides, which can be readily combined with metal parameter sets to model cluster–support interactions impacting catalysis on oxide-supported metal catalysts—this extension is an active area of ReaxFF development.
ReaxFF has recently been applied to study reactive events playing a critical role in the atomic layer deposition process (ALD), which is central to microchip manufacturing. ALD provides high conformity and excellent film thickness control in the ultrathin limit. Although ALD was introduced to the semiconductor industry quite some time ago as a cost-effective, scalable, low-temperature deposition process, a comprehensive and quantitative understanding of its kinetic mechanisms at the atomic scale has not yet been attained. This may be attributed to three main factors: (1) lack of suitable experimental probes to directly image or otherwise characterise the surface state with sufficient temporal and lateral resolution under realistic deposition conditions using spectroscopic or diffraction techniques; (2) challenges in computationally describing ALD relevant processes, such as adsorption, chemisorption, diffusion, steric hindrance, desorption and reaction dynamics; which are all strongly dependent on the ‘chemical state’ of the solid surface; and (3) lack of interest in acquiring such detailed knowledge, as it is not entirely necessary in order to utilise ALD for Si-based technology.
With the current drive to develop CMOS technology beyond Si into large-scale, manufacturable processes with high yield enabling high performance, low-power logic by replacing the channel using high mobility, low bandgap semiconductors (Ge or III–V),87,​88,​89 the situation has markedly changed. The integration of a thermodynamically stable high-ĸ dielectric film that is in contact with the semiconductor and forms an electrically well-behaved interface, allows for aggressive thickness scaling to achieve a large gate capacitance while limiting the gate leakage is a critical roadblock toward realizing this technology.90,​91,​92 Here an atomic scale understanding of the chemical state of the surface, suitable precursors and successful passivation strategies to avoid electrically active traps at the interface is critical. The ReaxFF methodology can help to rationalise and rapidly identify promising approaches, thus shortening the time required to develop this technology.
Recent studies have focused on ALD of Al2O3 using trimethylaluminum (TMA) and H2O cycles as a route for developing Ge-based metal oxide semiconductors (CMOS).89,93,94 To demonstrate the temperature dependence of TMA nucleation on a bare Ge(100) surface, which has been previously reported in the literature,95 MD simulations were performed at 23 °C, 227 °C, and 427 °C. The analysis of simulation trajectories revealed that TMA molecules non-dissociatively chemisorb on Ge dangling bonds at room temperature, whereas at higher temperature (427 °C) residual methyl groups on TMA molecules either recombine to form methane (and ethane) or bind to neighbouring Ge dangling bonds resulting in Ge-Al-Ge connections on the Ge(100) surface. These results are qualitatively consistent with scanning tunneling microscopy (Figure 3). In addition, the coverage of TMA sites with respect to the number Ge dangling bonds was evaluated as a function of temperature and compared to values obtained with X-ray photoelectron spectroscopy (XPS). Both theoretical and experimental approaches confirm that carbon coverage decreases with elevated temperature due to the dissociation of methyl groups from TMA, while Al coverage remains nearly constant, indicating that TMA sites become saturated under operating conditions.
Application of ReaxFF in MD simulations demonstrating TMA nucleation on a bare Ge(100) slab. (a) Initial configuration of 200 TMA molecules on a bare Ge(100) surface. (b) Comparison of TMA nucleation on a bare Ge(100) surface as modelled by ReaxFF and as observed by STM.95 (c) Estimated coverage of TMA sites obtained by ReaxFF simulations (left panel) and determined experimentally using XPS (right panel).95 (STM images adapted with permission from ref. 95 Copyright 2011 AIP Publishing LLC).
The recently developed ReaxFF description demonstrates the ability to correctly describe the TMA chemisorption on the Ge surface. This will enable future studies to assess ALD of Al2O3, as well as other feasible high-ĸ candidates and passivation strategies (such as hafnia, zirconia, tantalum, germanium oxynitride96,​97,​98,​99,​100 and bi/tri-layer dielectrics101,​102,​103,​104) on pristine and oxidised Ge surfaces. As such, ReaxFF can be instrumental in identifying optimal processing conditions to form high quality high-ĸ dielectrics/non-Si semiconductor interfaces.
The utility of ReaxFF lies in its ability to treat reactive processes at the interface between gas, liquid, and solid phases, which is pertinent when dealing with heterogeneous catalysis on metals and oxides (Figures 4a and b, respectively). This capability is not only relevant to catalysis, but rather is essential for modelling any nanoscale phenomena in complex systems. Onofrio et al.105 employed a modified ReaxFF potential to conduct atomistic-scale simulations of electrode-bridge-electrode systems (Figure 4c). By modifying ReaxFF to enable the direct simulation of an applied potential, the authors were able to model the formation and degradation of metallic filaments bridging electrodes in resistance-switching cells. Tavazza et al.106 conducted MD simulations of nanoindentation processes on Ni substrates, where ReaxFF was used to more accurately describe attractive interactions between diamond indenter tips and Ni/NiO surfaces. The utility of using ReaxFF for modelling reactive events at the interface between gas and solid phases is demonstrated in the work of Bagri et al.,107 who investigated defect formation during the reduction of graphene oxide (GO). MD simulations revealed reorganisation patterns in the graphene sheet during reduction (Figure 4d), allowing the authors to identify more efficient pathways for minimising defect formation during GO reduction. This is an important issue, as GO reduction is a promising synthesis route for the rapid production of graphene sheets. Similarly, Srinivasan et al.108,109 studied hyperthermal oxygen collisions with graphene, which replicates degradation processes affecting the performance of graphene-based heat shields during spacecraft re-entry. Oxidation of metal surfaces was investigated by Fantauzzi et al.110 and by Senftle et al.,111 studies in which MD simulations were employed to assess kinetic limitations affecting the initial oxidation of Pt and Pd surfaces, respectively (Figure 4e). ReaxFF has also been employed to model the kinetics of hydrogen adsorption and diffusion in various Pt,112 Pd,113 and Fe114 phases. The effect of solid-gas interactions is further demonstrated in the work of Raju et al.,52,115 in which the authors investigated oriented attachment mechanisms governing the organisation of TiO2 nanocrystals (Figure 4f). MD simulations revealed that, in vacuum, TiO2 crystals tend to form polycrystalline aggregates. In solution, approaching TiO2 crystals effectively dissociate water, thus creating a network of hydrogen bonds facilitating oriented attachment.
Application of the ReaxFF method to (a) Ni-catalysed CNT growth,83 (b) oxidative dehydrogenation over MMO catalysts,28 (c) electrometallisation cells,105 (d) reduction of graphene oxide,107 (e) Pd surface oxidation,111 (f) oriented attachment mechanisms in TiO2 nanocrystals,115 (g) conformational dynamics of biomolecules,53 (h) proton diffusion membranes51 and (i) capacitive mixing by double layer expansion.117 (Images adapted with permission from the provided references: a, Copyright 2013 American Physical Society; b, Copyright 2009 John Wiley and Sons; c, d and h, Copyright 2010 and 2015 Nature Publishing Group; e, Copyright 2013 AIP Publishing LLC; g, Copyright 2013 Royal Society of Chemistry; f and i, Copyright 2014 American Chemical Society).
In addition to gas phase chemistry, the development of the ‘aqueous branch’ parameter set has led to numerous ReaxFF studies involving reactive processes occurring in the liquid phase. Rahaman et al.54,116 developed a ReaxFF description capable of modelling glycine tautomerisation in water, which was later employed to investigate the interaction between glycine and TiO2 surfaces.54 Monti et al.53 further extended this parameter set to include amino acids and short peptide structures, allowing ReaxFF to simulate the conformational dynamics of biomolecules in solution (Figure 4g). Aqueous proton transfer across graphene was investigated by Achtyl et al.,51 where ReaxFF helped to establish that proton transfer is enabled by hydroxyl-terminated atomic defects in the graphene sheet (Figure 4h). Hatzell et al.117 used ReaxFF to determine the impact of strong acid functional groups on graphene electrodes in capacitive mixing devices, in which salinity gradients are used for power generation (Figure 4i). For numerous examples of water–solid interfaces that have been modelled with the ReaxFF potential, we refer the reader to references on the ‘aqueous branch’ shown in Figure 2.
MD approaches with reactive force-fields are tremendously useful for examining dynamic and kinetic properties of chemical processes on timescales beyond what is accessible with QM. However, numerous chemical processes, and in particular processes requiring thermal activation to overcome significant kinetic barriers, occur on timescales that are still beyond the simulation capabilities of empirical methods. Such limitations have provided the impetus for developing MC and accelerated molecular dynamics (AMD) approaches for the ReaxFF potential. MC methods are most useful when the properties of interest in a system are dominated by thermodynamic driving forces (i.e., when kinetic factors do not limit the deterministic evolution of the system). Conversely, AMD methods are suitable for modelling kinetic processes, as arguments from transition state theory are employed to accelerate system evolution between rare events while maintain system determinism. This section highlights four methods (two MC-based and two AMD-based) that have been recently introduced to ReaxFF: (1) hybrid Grand Canonical MC/Molecular dynamics (GC-MC/MD), (2) uniform-acceptance force-biased MC (UFMC), (3) parallel replica dynamics (PRD), and (4) adaptive accelerated ReaxFF reaction dynamics (aARRDyn).
MC methods in the grand canonical ensemble stochastically insert, remove, and displace MC-atoms in the host system. The equilibrated number of MC-atoms, as well as system structure, is dictated by the chemical potential of the MC-atom reservoir (i.e., the free energy per MC-atom in the original source, such as ½GO2 for modelling oxide formation upon exposure to O2). When applying GC-MC with the ReaxFF potential, we have introduced a force relaxation step (via either low-temperature MD or similar energy minimisation algorithms) that allows the structure of the host system to locally reorganise as MC-steps are executed (i.e., to find the nearest local energy minimum). The introduction of this MD step, motivated by the MC/MD methodology employed by Chenoweth et al.,28 is essential for reaching an equilibrated structure. Otherwise the majority of MC moves involving the addition of an atom into a bulk region will be rejected because such moves inevitably result in a very high energy. Slow acceptance rates are effectively mitigated by allowing the system to restructure as it accommodates the inserted MC-atoms, mimicking the kinetic restructuring of the system as the new phase is formed.
The GC-MC/MD approach was initially developed to model the extent of oxidation in Pd nanoclusters, as PdO formation significantly affects the performance of Pd-based catalysts.111 To validate the GC-MC/MD approach, we derived an oxidation phase diagram for Pd clusters, where the predicted T, P ranges of PdO to Pd transitions were in agreement with experimental observations (Figure 5a). Coupled with MD simulations capturing the extent of oxygen migration through various Pd surface facets, this study demonstrated the utility of ReaxFF for investigating both kinetic and thermodynamic driving forces affecting phase stability. A similar approach has been applied to investigate the initial stages of oxide formation on Pt surfaces.110 The method has additionally been applied to investigate hydrogen and carbon uptake by palladium, where it successfully captured phase transitions in bulk, cluster and oxide-supported systems.113,118,​119,​120 It has also been used to produce open-circuit voltage profiles for Li intercalation in graphite121 and sulfur122 electrodes during cell discharge.
New simulation techniques for the ReaxFF potential: (a) grand canonical Monte Carlo,111 (b) parallel replica dynamics,126 (c) uniform-acceptance force-biased Monte Carlo,85 and (d) adaptive accelerated molecular dynamics.127 (Images adapted with permission from the provided references: a, Copyright 2013 AIP Publishing LLC; b–d, Copyright 2010, 2013, and 2014 American Chemical Society).
Applied to catalysis, GC-MC/MD provides an approach for predicting the thermodynamic stability of phases likely to form upon exposure to the reactant mixture, which can be used to establish plausible surface models that better approximate the structure of the catalyst under reaction conditions. To investigate kinetics, the derived surface models can then be directly employed in MD simulations of the catalytic process, or can serve to motivate smaller scale models for more accurate DFT calculations. More importantly, it provides a means to explore the phase space of a system without a priori knowledge of the phase diagram in question—making it a powerful low-CPU-cost tool when searching for novel system properties.
Pure MC methods are stochastic in nature, and correspondingly cannot be used to evaluate the deterministic evolution of a system in time. However, stochastic approaches can be incorporated in deterministic MD to accelerate reaction steps past kinetic barriers, decreasing the number of iterations required to model state-to-state transitions in slower processes. Although information regarding the timescale of such processes may be lost, the deterministic nature of the simulation is preserved. This is valuable when one wishes to understand how the system evolves, rather than how fast the system evolves. In UFMC simulations,123 stochastic spatial displacements (accepted at every iteration) are applied along the force acting on each atom, where the extent and direction of the displacement is controlled by both the magnitude of the force and the overall system temperature. The displacement is restricted along the acting force when large forces are present or the system is at low temperature, whereas displacements are fully random when small forces are acting or the system is at high temperature. Uniform-acceptance prevents the system from becoming restricted to its current potential energy basin, while force-bias ensures a preferential evolution along the lowest-force reaction coordinate.
This strategy, called uniform-acceptance force-biased MC (UFMC), was adapted for ReaxFF by Neyts et al.84,​85,​86,124 to investigate the initial stages of CNT growth over Ni clusters. During simulations of CNT nucleation, UFMC simulation stages were incorporated between MD runs to allow deposited carbon atoms to relax into chiral CNT structures. Such relaxation steps with pure MD would necessitate unphysical carbon deposition rates, as well as intractable simulation timeframes. As described above, this methodology was also successfully employed to assess the effect of applied potentials and Ar+ ion bombardment on CNT nucleation and growth (Figure 5c).
aARRDyn, introduced to ReaxFF by Cheng et al.,127 employs a ‘bond boost’ algorithm to accelerate dynamic evolution across the potential energy surface. In this implementation, the bond-order formalism of ReaxFF is aptly suited for identifying elongated bonds that are approaching a transition state, where the boost potential is calculated using an envelope function128 that directs energy toward bonds that are closest to reacting. The authors demonstrated the aARRDyn methodology by modelling H2/O2 ignition in simulations reaching ~10 ns, which successfully reproduced previously determined reaction pathways (Figure 5d). Altogether, PRD and aARRDyn demonstrate the extension of ReaxFF to low-temperature applications, where kinetic processes occur on timescales inaccessible to traditional MD trajectory integration.
A polarisable charge calculation method is essential for a transferable reactive force-field method. ReaxFF typically employs the electronegativity equalisation method (EEM) developed by Mortier et al.129 Unfortunately, EEM has a few drawbacks,130,​131,​132,​133,​134 the worst being its inability to restrain long-range charge-transfer, even between molecular fragments that are well separated. This issue is especially apparent in low-density gas phase simulations, where EEM-ReaxFF predicts small, but certainly nonzero, charges on isolated molecular species, which significantly affects accommodation coefficients.135 In simulations of dense systems, unrealistic charge-transfer may also occur (e.g., between two dielectric phases with a different intrinsic electronegativity). To ameliorate this issue, we have recently incorporated atom-condensed Kohn-Sham DFT approximated to second order (ACKS2) in ReaxFF.136,137 ACKS2 is an extension of EEM that penalises long-range charge transfer with a bond-polarisation energy, in line with the split-charge equilibration (SQE) model.138 The bond polarisability is a function of interatomic distance, which slightly increases beyond the equilibrium bond length but then quickly decays to zero, effectively enforcing fragment neutrality. Transferring from EEM-ReaxFF to ACKS2-ReaxFF does require reparameterisation, although the EEM-ReaxFF parameters are typically a very good starting point for deriving ACKS2-ReaxFF parameters. As such, redevelopment is relatively straightforward. We believe that ACKS2-ReaxFF is most relevant for obtaining reliable accommodation coefficients, as well as for describing physisorption of neutral and ionic molecules on surfaces. It is also vital for incorporating explicit electronic degrees of freedom into ReaxFF, as described in the next section.
The treatment of explicit charge and polarisation is essential for force-field descriptions applied to systems such as rechargeable battery interfaces and ferro-/piezoelectric materials. Essentially, such descriptions require a classical treatment for an explicit electron. Recently, potentials including some form of explicit electron description have been introduced, such as the electron force-field139 and the LEWIS140 force-field. Nonetheless, these methods have not yet been demonstrated to accurately simulate complex materials and intricate chemistries.
To extend the ReaxFF description to include chemistry dependent on electron diffusion, we have introduced explicit electron-like and hole-like particles that carry negative (−1) and positive (+1) charges, respectively. The electron and hole particles interact with atomic centers through a single Gaussian function.139 We implemented charge-valence coupling, which allows the electron or hole particle to modify the number of valence electrons in a host atom, thus ensuring the appropriate change in valence when calculating the degree of over or under coordination in the host atom. To demonstrate the capability of the electron-explicit version of ReaxFF (eReaxFF), we trained our force-field to capture the electron affinity (EA) of various hydrocarbon species. Figure 6 summarises the performance of eReaxFF compared with both standard ReaxFF and available literature data.141,142 eReaxFF qualitatively reproduces the literature data, whereas standard ReaxFF entirely fails to capture the EA of most species considered in this training set. Still, eReaxFF significantly underestimates the EA of ethylene and propene compared with experimental data.142 To further investigate this discrepancy, we performed DFT calculations with the M06-2X functional (aug-cc-pVTZ basis set) implemented in Jaguar 7.5,143 which shows that DFT also underestimates the EA of these two species. Overall, eReaxFF provides a significant improvement in capturing EAs in comparison with standard ReaxFF. Future development of eReaxFF will focus on the description of interfacial chemistry in batteries and on structural/polarisation behaviour in piezo/ferro-electric materials.
Most ReaxFF parameter development is currently performed using a straightforward, single-parameter parabolic search method.144 This method has a high level of transparency—which, given the complexity of the ReaxFF functional form and the typically large size of its training set, is a significant advantage. However, this method is clearly not the most efficient for sampling the complex force-field error landscape, as it does not implicitly include parameter correlation and has a significant risk of becoming ensnared in local energy minima. As such, various groups have recently focused on developing more sophisticated parameter sampling methods. In particular, MC and genetic algorithm (GA) based approaches have risen in popularity.145,​146,​147,​148,​149 So far, these methods have only been used for isolated force-field development efforts, and as such they still have to demonstrate their transferability and user friendliness. They do hold promise for reducing force-field development time, as well as for making the process more accessible to nonexpert users.
A number of different ReaxFF implementations have been developed over the years. The first-generation ReaxFF implementation of van Duin et al.1 established the utility of the force field in the context of various applications. This serial, fortran-77, implementation was integrated into the publicly available, open-source LAMMPS code150 by Thompson et al.151 as the Reax package152—the first publicly available parallel implementation of ReaxFF. Nomura et al. have developed a parallel ReaxFF implementation, which has been used in a number of large-scale simulations, including high-energy materials, metal grain boundary decohesion, water bubbles and surface chemistry.153,​154,​155,​156,​157,​158,​159 This Nomura et al. code is not publicly available.
The PuReMD software package, and its integration into LAMMPS, which is available as the USER-REAXC package, represent the state-of-the-art open-source implementations of ReaxFF. Since PuReMD and USER-REAXC are, to the best of our knowledge, the most widely used open-source codes for ReaxFF simulations, we summarise algorithmic and software design choices, as well as future directions in PuReMD. PuReMD uses novel algorithms and data structures to achieve high performance in force computations, with a small memory footprint. An optimised binning-based neighbour generation method, elimination of the bond-order derivatives list in bonded interactions, lookup tables to accelerate non-bonded interaction computations, and preconditioned fast iterative solvers for the charge equilibration problem are the major algorithmic innovations in PuReMD.10,12,15 The dynamic nature of the bond, three-body and four-body interactions in a reactive molecular system presents challenges in terms of memory management and data structures for efficiently computing bonded interactions. PuReMD introduces novel data structures to store three-body and four-body interactions in a compact form. Its dynamic memory management system automatically adapts to the needs of the input system over the course of a simulation, significantly reducing the memory footprint and minimizing the effort to setup a simulation. Its parallel formulation, called PuReMD, uses conventional message passing (MPI) to enable the simulation of large molecular systems12 on scalable parallel platforms.
Recent efforts by us and others have focused on enabling fast ReaxFF simulations on GPUs. Zheng et al.160 reported the first GPU implementation of ReaxFF, called GMD-Reax. GMD-Reax is reported to be up to six times faster than the USER-REAXC package on small systems, and about 2.5 times faster on typical workloads using a quad core Intel Xeon CPU. However, this performance is significantly predicated on the use of single-precision arithmetic operations (PuReMD codebase is fully double precision) and low-precision transcendental functions, which can potentially lead to significant energy drifts in long NVE simulations. GMD-Reax can run on single GPUs only, and is not publicly available.
Recently, Kylasa et al. have developed a publically available GPU version of the PuReMD code, called PuReMD-GPU10, which uses CUDA along with an extensive set of optimisations to deliver significant speedups on single GPU systems. For a variety of benchmarks, this code shows up to a 16-fold speedup on an Nvidia C2075 GPU, compared with a single core of an Intel Xeon processor.
Due to device memory limitations, single GPU simulations are typically limited to 30 to 40 thousand atoms. Our more recent development efforts have resulted in the PuReMD-PGPU code, which builds on PuReMD-GPU to enable ReaxFF simulations on multiple GPUs. PuReMD-PGPU has been benchmarked on Texas Advanced Computing Center’s Maverick cluster, which has 132 nodes each with two Intel E5–2680 v2 Ivy Bridge processors running at 2.86 GHz and one NVIDIA Tesla K40 GPU. We ran NVT MD calculations of DNA strands suspended in water, having a total of 531,494 atoms. In Figure 7, we present the performance of this simulation on 6, 8, 12, 15 and 18 nodes, where we use either 1 GPU per node to perform the calculations using PuReMD-PGPU or the MPI-based CPU code (PuReMD) with a varying number of cores per node. In terms of speed, we find that using one GPU per node with PuReMD-PGPU is equivalent to using 8 to 10 cores when using PuReMD.
Performance of the PuReMD implementation of ReaxFF. Time required to run 1,000 time steps of a MD simulation on half a million atoms where either GPU or CPU is used for computing the interatomic forces.
More recently, a Hybrid-PuReMD, an MPI/pthreads/CUDA implementation by Kylasa et al., leverages all the available compute resources on GPU-equipped nodes. Hybrid-PuReMD uses an effective task-parallel implementation by overlapping the execution of bonded interactions on CPU and non-bonded and charge equilibration computations on the GPU. This code has been benchmarked on Michigan State University’s High Performance Computing Center GPU cluster, which has 20 nodes each with two 2.5 GHz 10-core Intel Xeon E5–2670v2 processors and two NVIDIA K20 GPUs. Figures 8 and 9 plot the speedup achieved by Hybrid-PuReMD compared with PuReMD (20 MPI processes per node) and total time per time-step for weak scaling scenarios. We notice that for bulk water system with 100 K atoms per node Hybrid-PuReMD effectively delivers a speedup of over 3.3-fold over PuReMD on a per-node basis (i.e., 3.3 times the performance of PuReMD running on all 20 cores of the node)—achieving near optimal resource utilisation.
Speedup achieved by Hybrid-PuReMD compared to PuReMD with 20 MPI processors per node.
Total time per time-step comparison between Hybrid-PuReMD and PuReMD for weak scaling scenarios.
Other ongoing developments with PuReMD includes performance optimisation for many-core architectures, such as IBM BGQ and Intel Xeon Phi processors. Using a hybrid MPI/OpenMP parallel implementation, Aktulga et al. have achieved speedups of 1.5–4.5× over the current USER-REAXC implementation for PETN crystal benchmarks of sizes ranging from 32 thousand to 16.6 million atoms161. This study has been conducted using up to 16,384 nodes (262 K cores) on Mira, an IBM BGQ supercomputer.
While promising, performance results obtained on GPUs and many-core CPUs suggest that there is still room for performance improvements. Our future work will focus on achieving improved performance on GPUs, other accelerators, and many-core CPUs by better managing on-node parallelism, exploiting data locality, and the development of scalable solvers for the charge equilibration problem, which represents a major bottleneck in massively parallel runs.
In this article, we reviewed the development, application, and future directions of the ReaxFF method. ReaxFF helps to bridge the gap in simulation scale separating QM and classical methods. By employing a bond-order formalism within a classical approach, ReaxFF implicitly describes chemical bonding without expensive QM calculations. We provided an overview of ReaxFF development history, lending insight into the development choices that shape the currently employed ReaxFF formalism. In addition, we discussed numerous applications of the method, which include processes ranging from the combustion of coal to the conformational dynamics of biomolecules. A particular emphasis was placed on studies that use ReaxFF to investigate reactive events at the interphase between solid, liquid and gas phases, thus demonstrating the diverse transferability of the ReaxFF formalism. In addition, we reviewed the recent development and application of MC and accelerated molecular dynamics methods tailored to the bond-length formalism of ReaxFF. Finally, we provided an overview of future research directions we will be taking, as we seek to improve the performance of ReaxFF, as well as to extend the range of chemical phenomena that can be studied with the method.
van Duin, A. C. T., Dasgupta, S., Lorant, F. & Goddard, III W. A. ReaxFF: a reactive force field for hydrocarbons . J. Phys. Chem. A 105, 9396–9409 (2001).
Chenoweth, K., van Duin, A. C. T. & Goddard, W. A. ReaxFF reactive force field for molecular dynamics simulations of hydrocarbon oxidation . J. Phys. Chem. A 112, 1040–1053 (2008).
van Duin, A. C. T. et al. ReaxFFSiO reactive force field for silicon and silicon oxide systems . J. Phys. Chem. A 107, 3803–3811 (2003).
Strachan, A., van Duin, A. C. T., Chakraborty, D., Dasgupta, S. & Goddard, W. A. Shock waves in high-energy materials: The initial chemical events in nitramine RDX . Phys. Rev. Lett. 91, 098301 (2003).
Strachan, A., Kober, E. M., van Duin, A. C. T., Oxgaard, J. & Goddard, W. A. Thermal decomposition of RDX from reactive molecular dynamics . J. Chem. Phys. 122, 054502 (2005).
Cheung, S., Deng, W. Q., van Duin, A. C. T. & Goddard, W. A. ReaxFF(MgH) reactive force field for magnesium hydride systems . J. Phys. Chem. A 109, 851–859 (2005).
van Duin, A. C. T. et al. Development and validation of a ReaxFF reactive force field for Cu cation/water interactions and copper metal/metal oxide/metal hydroxide condensed phases . J. Phys. Chem. A 114, 9507–9514 (2010).
Liu, L., Liu, Y., Zybin, S. V., Sun, H. & Goddard, W. A. ReaxFF-lg: correction of the ReaxFF reactive force field for London dispersion, with applications to the equations of state for energetic materials . J. Phys. Chem. A 115, 11016–11022 (2011).
Plimpton, S. J. & Thompson, A. P. Computational aspects of many-body potentials . MRS Bull. 37, 513–521 (2012).
Kylasa, S. B., Aktulga, H. M. & Grama, A. Y. PuReMD-GPU: A reactive molecular dynamics simulation package for GPUs . J. Comput. Phys. 272, 343–359 (2014).
Aktulga, H. M., Pandit, S. A., van Duin, A. C. T. & Grama, A. Y. Reactive molecular dynamics: numerical methods and algorithmic techniques . SIAM J. Sci. Comput. 34, C1–C23 (2012).
Kamat, A. M., van Duin, A. C. T. & Yakovlev, A. Molecular dynamics simulations of laser-induced incandescence of soot using an extended reaxff reactive force field . J. Phys. Chem. A 114, 12561–12572 (2010).
Zhang, Q. et al. Adhesion and nonwetting-wetting transition in the Al/α-Al2O3 interface . Phys. Rev. B 69, 045423 (2004).
Fogarty, J. C., Aktulga, H. M., Grama, A. Y., van Duin, A. C. T. & Pandit, S. A. A reactive molecular dynamics simulation of the silica-water interface . J. Chem. Phys. 132, 174704 (2010).
Neyts, E. C., Khalilov, U., Portois, G. & van Duin, A. C. T. Hyperthermal oxygen interacting with silicon surfaces: adsorption, implantation and damage creation . J. Phys. Chem. C 115, 4818–4823 (2011).
Castro-Marcano, F. & van Duin, A. Comparison of thermal and catalytic cracking of hydrocarbon fuel from ReaxFF reactive molecular dynamics simulations . Combust. Flame 160, 766–775 (2013).
Sen, F. G., Alpas, A. T., van Duin, A. C. T. & Qi, Y. Oxidation assisted ductility in aluminum nanowires . Nat. Commun. 5, 3959 (2014).
Russo, M., Li, R., Mench, M. & van Duin, A. C. T. Molecular dynamic simulation of aluminum-water reactions using the ReaxFF reactive force field . Int. J. Hydrogen Energ. 36, 5828–5835 (2011).
Abdolhosseini Qomi, M. J. et al. Combinatorial molecular optimization of cement hydrates . Nat. Commun. 5, 4960 (2014).
Joshi, K., Psofogiannakis, G., Raman, S. & van Duin, A. C. T. Reactive molecular simulations of protonation of water clusters and depletion of acidity in H-ZSM-5 zeolite . Phys. Chem. Chem. Phys. 16, 18433–18441 (2014).
Joshi, K. & van Duin, A. C. T. Molecular dynamics study on the influence of additives on the high temperature structural and acidic properties of ZSM-5 zeolite . Energ. Fuel. 27, 4481–4488 (2013).
Pitman, M. C. & van Duin, A. C. T. Dynamics of confined reactive water in smectite clay-zeolite composites . J. Am. Chem. Soc. 134, 3042–3053 (2012).
Russo, M. F. Jr & van Duin, A. C. T. Atomistic-scale simulations of chemical reactions: Bridging from quantum chemistry to engineering . Nucl. Instrum. Methods Phys. Res. B 269, 1549–1554 (2011).
Case, D. A. et al. The Amber biomolecular simulation programs . J. Comput. Chem. 26, 1668–1688 (2005).
Chenoweth, K., van Duin, A. C. T. & Goddard, W. A. The ReaxFF Monte Carlo reactive dynamics method for predicting atomistic structures of disordered ceramics: application to the Mo3VOx Catalyst . Angew. Chem. Int. Ed. 48, 7630–7634 (2009).
Goddard, W. A., Chenoweth, K., Pudar, S., van Duin, A. C. T. & Cheng, M. J. Structures, mechanisms, and kinetics of selective ammoxidation and oxidation of propane over multi-metal oxide catalysts . Top. Catal. 50, 2–18 (2008).
Buehler, M. J., van Duin, A. C. T. & Goddard, W. A. Multiparadigm modeling of dynamical crack propagation in silicon using a reactive force field . Phys. Rev. Lett., 96, 095505 (2006).
Buehler, M. J., Tang, H., van Duin, A. C. T. & Goddard, W. A. Threshold crack speed controls dynamical fracture of silicon single crystals . Phys. Rev. Lett. 99, 165502 (2007).
Weismiller, M. R., Duin, A. C. T. v., Lee, J. & Yetter, R. A. ReaxFF reactive force field development and applications for molecular dynamics simulations of ammonia borane dehydrogenation and combustion . J. Phys. Chem. A 114, 5485–5492 (2010).
Mueller, J. E., van Duin, A. C. T. & Goddard, W. A. Application of the ReaxFF reactive force field to reactive dynamics of hydrocarbon chemisorption and decomposition . J. Phys. Chem. C 114, 5675–5685 (2010).
Paupitz, R., Junkermeier, C. E., van Duin, A. C. T. & Branicio, P. Fullerenes generated from porous structures . Phys. Chem. Chem. Phys. 16, 25515–25522 (2014).
Verners, O. & van Duin, A. C. T. Comparative molecular dynamics study of fcc-Ni nanoplate stress corrosion in water . Surf. Sci. 633, 94–101 (2015).
Zou, C., Shin, Y. K., Liu, Z. K., Fang, H. Z. & van Duin, A. C. T. Molecular dynamics simulations of the vacancy effect on the nickel self diffusion, oxygen diffusion and oxidation initiation in nickel using the ReaxFF reactive force field . Acta Mater. 83, 102–112 (2015).
Fang, H. Z. et al. First-Principles Studies on Vacancy-modified Interstitial Diffusion Mechanism of Oxygen in Nickel, Associated with Large-Scale Atomic Simulation Techniques . J. Appl. Phys. 115, 043501 (2014).
Han, S.-p., Strachan, A., van Duin, A. C. T. & Goddard, W. A. Thermal decomposition of condense-phase nitromethane from molecular dynamics using reaxff reactive dynamics . J. Phys. Chem. B 115, 6534–6540 (2011).
Cherukara, M. J., Wood, M. A., Kober, E. M. & Strachan, A. Ultra-fast chemistry under non-equilibrium conditions and the shock to deflagration transition at the nanoscale . J. Phys. Chem. C 119, 22008–22015 (2015).
Guo, F., Zhang, H. & Cheng, X. Molecular dynamic simulations of solid nitromethane under high pressures . J. Theor. Comput. Chem. 9, 315–325 (2010).
Zhang, L. Z. et al. Carbon cluster formation during thermal decomposition of octahydro-1,3,5,7-tetranitro-1,3,5,7-tetrazocine and 1,3,5-triamino-2,4,6-trinitrobenzene high explosives from ReaxFF reactive molecular dynamics simulations . J. Phys. Chem. A 113, 10619–10640 (2009).
Zhang, L. Z., van Duin, A. C. T., Zybin, S. V. & Goddard, W. A. Thermal decomposition of hydrazines from reactive dynamics using the ReaxFF reactive force field . J. Phys. Chem. B 113, 10770–10778 (2009).
Rom, N. et al. Change of mechanism in the decomposition of hot dense liquid nitromethane as a function of density: MD simulations based on REAXFF . J. Phys. Chem. A 115, 10181–10202 (2011).
Chen, N., Lusk, M. T., van Duin, A. C. T. & Goddard, W. A. Mechanical properties of connected carbon nanorings via molecular dynamics simulation . Phys. Rev. B 72, 085416 (2005).
Wood, M. A., van Duin, A. C. T. & Strachan, A. Coupled Thermal and Electromagnetic Induced Decomposition in the Molecular Explosive αHMX; A Reactive Molecular Dynamics Study . J. Phys. Chem. A 118, 885–895 (2014).
Kim, S.-Y. et al. Development of a ReaxFF reactive force field for Titanium dioxide/water systems . Langmuir 29, 7838–7846 (2013).
Pitman, M. C. & van Duin, A. C. T. Dynamics of confined reactive water in smectic clay-zeolite composites . J. Am. Chem. Soc. 134, 3042–3053 (2012).
Gale, J. D., Ratieri, P. & van Duin, A. C. T. A reactive force field for aqueous-calcium carbonate systems . Phys. Chem. Chem. Phys. 13, 16666–16679 (2011).
Aryanpour, M., van Duin, A. C. T. & Kubicki, J. D. Development of a reactive force field for iron-oxyhydroxide systems . J. Phys. Chem. A 114, 6298–6307 (2010).
Raymand, D., van Duin, A. C. T., Spångberg, D., Goddard, W. A. III & Hermansson, K. Water adsorption on stepped ZnO surfaces from MD simulation . Surf. Sci. 604, 741–752 (2010).
Achtyl, J. L. et al. Aqueous proton transfer across single-layer graphene . Nat. Commun. 6, 6539 (2015).
Monti, S. et al. Exploring the conformational and reactive dynamics of biomolecules in solution using an extended version of the glycine reactive force field . Phys. Chem. Chem. Phys. 15, 15062–15077 (2013).
Monti, S., van Duin, A. C. T., Kim, S.-Y. & Barone, V. Exploration of the conformational and reactive dynamics of glycine and diglycine on TiO2: computational investigations in the gas phase and in solution . J. Phys. Chem. C 116, 5141–5150 (2012).
Yue, D.-C. et al. Tribochemistry of phosphoric acid sheared between quartz surfaces: a reactive molecular dynamics study . J. Phys. Chem. C 117, 25604–25614 (2013).
Rahaman, O. et al. Development of a ReaxFF reactive force field for aqueous chloride and copper chloride . J. Phys. Chem. A 114, 3556–3568 (2010).
Zhu, R. et al. Characterization of the active site of yeast RNA polymerase II by DFT and ReaxFF calculations . Theor. Chem. Acc. 120, 479–489 (2008).
Yusupov, M. et al. Atomic scale simulations of plasma species interacting with bacteria cell walls . New. J. Phys. 14, 093043 (2012).
Abolfath, R. M., van Duin, A. C. T., Biswas, P. & Brabec, T. Reactive Molecular Dynamics study on the first steps of DNA-damage by free hydroxyl radicals . J. Phys. Chem. A 115, 11045–11049 (2011).
Keten, S., Chou, C.-C., van Duin, A. C. T. & Buehler, M. J. Tunable nanomechanics of protein disulphide bond begets weakening in reducing and stabilization in oxidizing chemical microenvironments . J. Mech. Behav. Biochem. Mater. 5, 32–40 (2012).
Pettifor, D. G. New Many-Body Potential for the Bond Order . Phys. Rev. Lett. 63, 2480–2483 (1989).
Stewart, J. P. MOPAC: A semiempirical molecular orbital program . J. Comput. Aid. Mol. Des. 4, 1–103 (1990).
Porezag, D., Frauenheim, T., Köhler, T., Seifert, G. & Kaschner, R. Construction of tight-binding-like potentials on the basis of density-functional theory: Application to carbon . Phys. Rev. B 51, 12947–12957 (1995).
Lewis, J. P. et al. Advances and applications in the FIREBALLab initio tight-binding molecular-dynamics formalism . Phys. Stat. Solid. B 248, 1989–2007 (2011).
Lee, B.-J. & Baskes, M. I. Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals . Phys. Rev. B 62, 8564–8567 (2000).
Daw, M. S. & Baskes, M. I. Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals . Phys. Rev. B 29, 6443–6453 (1984).
Baskes, M. I. Modified embedded-atom potentials for cubic materials and impurities . Phys. Rev. B 46, 2727–2742 (1992).
Abell, G. C. Empirical chemical pseudopotential theory of molecular and metallic bonding . Phys. Rev. B 31, 6184–6196 (1985).
Tersoff, J. New empirical approach for the structure and energy of covalent systems . Phys. Rev. B 37, 6991–7000 (1988).
Tersoff, J. Empirical interatomic potential for carbon, with applications to amorphous carbon . Phys. Rev. Lett. 61, 2879–2882 (1988).
Brenner, D. W. Empirical potential for hydrocarbons for use in simulating the chemical vapor deposition of diamond films . Phys. Rev. B 42, 9458–9471 (1990).
Stuart, S. J., Tutein, A. B. & Harrison, J. A. A reactive potential for hydrocarbons with intermolecular interactions . J. Chem. Phys. 112, 6472 (2000).
Liang, T. et al. Reactive Potentials for Advanced Atomistic Simulations . Annu. Rev. Mater. Sci. 43, 109–129 (2013).
Shin, Y. K. et al. Variable Charge Many-Body Interatomic Potentials . MRS Rev. 37, 504–512 (2012).
Fonseca, A. F. et al. Reparameterization of the REBO-CHO potential for graphene oxide molecular dynamics simulations . Phys. Rev. B 84, 075460 (2011).
Shan, T.-R. et al. Second generation charge optimized many-body (COMB) Potential for Si/SiO2 and amorphous silica . Phys. Rev. B 82, 235302 (2010).
Phillpot, S. R. & Sinnott, S. B. Simulating multifunctional structures . Science 325, 1634–1635 (2009).
Akimov, A. V. & Prezhdo, O. V. Large-scale computations in chemistry: a bird’s eye view of a vibrant field . Chem. Rev. 115, 5797–5890 (2015).
Qian, H.-J., van Duin, A. C. T., Morokuma, K. & Irle, S. Reactive molecular dynamics simulation of fullerene combustion synthesis: ReaxFF vs DFTB potentials . J. Chem. Theor. Comput. 7, 2040–2048 (2011).
Raymand, D., van Duin, A. C. T., Baudin, M. & Hermansson, K. A reactive force field (ReaxFF) for zinc oxide . Surf. Sci. 602, 1020–1031 (2008).
LaBrosse, M. R., Johnson, J. K. & van Duin, A. C. T. Development of a transferable reactive force field for cobalt . J. Phys. Chem. A 114, 5855–5861 (2010).
Goddard, W. A., Mueller, J. E. & van Duin, A. C. T. Development and validation of ReaxFF reactive force field for hydrocarbon chemistry catalyzed by nickel . J. Phys. Chem. C 114, 4939–4949 (2010).
Neyts, E. C. et al. Defect healing and enhanced nucleation of carbon nanotubes by low-energy ion bombardment . Phys. Rev. Lett. 110, 065501 (2013).
Neyts, E. C., van Duin, A. C. T. & Bogaerts, A. Insights in the plasma-assisted growth of carbon nanotubes through atomic scale simulations: effect of electric field . J. Am. Chem. Soc. 134, 1256–1260 (2012).
Neyts, E. C., Shibuta, Y., van Duin, A. C. T. & Bogaerts, A. Catalyzed growth of carbon nanotube with definable chirality by hybrid molecular dynamics-force biased monte carlo simulations . ACS Nano 4, 6665–6672 (2010).
del Alamo, J. A. Nanometre-scale electronics with III-V compound semiconductors . Nature 479, 317–323 (2011).
Frank, M. M. in 2011 Proceedings of the ESSCIRC (ESSCIRC) 50–58 (IEEE, Helsinki, 2011).
Kamata, Y. High-k/Ge MOSFETs for future nanoelectronics . Mater. Today 11, 30–38 (2008).
Houssa, M., Chagarov, E. & Kummel, A. Surface defects and passivation of Ge and III-V interfaces . MRS Bull. 34, 504–513 (2009).
Wallace, R. M., McIntyre, P. C., Kim, J. & Nishi, Y. Atomic layer deposition of dielectrics on ge and iii-v materials for ultrahigh performance transistors . MRS Bull. 34, 493–503 (2009).
Engel-Herbert, R., Hwang, Y. & Stemmer, S. Comparison of methods to quantify interface trap densities at dielectric/III-V semiconductor interfaces . J. Appl. Phys. 108, 124101 (2010).
George, S. M. Atomic layer deposition: an overview . Chem. Rev. 110, 111–131 (2009).
Swaminathan, S., Sun, Y., Pianetta, P. & McIntyre, P. C. Ultrathin ALD-Al2O3 layers for Ge (001) gate stacks: Local composition evolution and dielectric properties . J. Appl. Phys. 110, 094105 (2011).
Lee, J. S. et al. Atomic imaging of nucleation of trimethylaluminum on clean and H2O functionalized Ge(100) surfaces . J. Chem. Phys. 135, 054705 (2011).
Fadida, S. et al. Hf-based high-k dielectrics for p-Ge MOS gate stacks . J. Vac. Sci. Technol. B 32, 03D105 (2014).
Seo, K.-I. et al. Chemical states and electronic structure of a HfO2∕Ge(001) interface . Appl. Phys. Lett. 87, 042902 (2005).
Yoshiki, K., Yuuichi, K., Tsunehiro, I. & Akira, N. Direct Comparison of ZrO2 and HfO2 on Ge Substrate in Terms of the Realization of Ultrathin High-κ Gate Stacks . Jpn J. Appl. Phys. 44, 2323 (2005).
Kutsuki, K., Okamoto, G., Hosoi, T., Shimura, T. & Watanabe, H. Germanium oxynitride gate dielectrics formed by plasma nitridation of ultrathin thermal oxides on Ge(100) . Appl. Phys. Lett. 95, 022102 (2009).
On Chui, C., Ramanathan, S., Triplett, B. B., McIntyre, P. C. & Saraswat, K. C. Germanium MOS capacitors incorporating ultrathin high-/spl kappa/ gate dielectric . IEEE Electr. Device Lett. 23, 473–475 (2002).
Swaminathan, S., Shandalov, M., Oshima, Y. & McIntyre, P. C. Bilayer metal oxide gate insulators for scaled Ge-channel metal-oxide-semiconductor devices . Appl. Phys. Lett. 96, 082904 (2010).
Kana, H. et al. Fabrication of Ge Metal-Oxide-Semiconductor Capacitors with High-Quality Interface by Ultrathin SiO2/GeO2 Bilayer Passivation and Postmetallization Annealing Effect of Al . Jpn J. Appl. Phys. 50, 04DA10 (2011).
Xie, Q. et al. Implementing TiO2 as gate dielectric for Ge-channel complementary metal-oxide-semiconductor devices by using HfO2/GeO2 interlayer . Appl. Phys. Lett. 97, 112905 (2010).
Zheng, Y. X. et al. In situ process control of trilayer gate-stacks on p-germanium with 0.85-nm EOT . IEEE Electr. Device Lett. 36, 881–883 (2015).
Onofrio, N., Guzman, D. & Strachan, A. Atomic origin of ultrafast resistance switching in nanoscale electrometallization cells . Nat. Mater. 14, 440–446 (2015).
Tavazza, F., Senftle, T. P., Zou, C., Becker, C. A. & van Duin, A. C. T. Molecular dynamics investigation of the effects of tip-substrate interactions during nanoindentation . J. Phys. Chem. C 119, 13580–13589 (2015).
Bagri, A. et al. Structural evolution during the reduction of chemically derived graphene oxide . Nat. Chem. 2, 581–587 (2010).
Srinivasan, S. G., van Duin, A. C. T. & Ganesh, P. Development of a ReaxFF potential for carbon condensed phases and its application to the thermal fragmentation of a large fullerene . J. Phys. Chem. A 119, 571–580 (2015).
Goverapet Srinivasan, S. & van Duin, A. C. T. Molecular-dynamics-based study of the collisions of hyperthermal atomic oxygen with graphene using the ReaxFF reactive force field . J. Phys. Chem. A 115, 13269–13280 (2011).
Fantauzzi, D. et al. Development of a ReaxFF potential for Pt-O systems describing the energetics and dynamics of Pt-oxide formation . Phys. Chem. Chem. Phys. 16, 23118–23133 (2014).
Senftle, T. P., Meyer, R. J., Janik, M. J. & van Duin, A. C. T. Development of a ReaxFF potential for Pd/O and application to palladium oxide formation . J. Chem. Phys. 139, 044109–044115 (2013).
Ludwig, J., Vlachos, D. G., van Duin, A. C. T. & Goddard, W. A. Dynamics of the dissociation of hydrogen on stepped platinum surfaces using the reaxff reactive force field . J. Phys. Chem. B 110, 4274–4282 (2006).
Senftle, T. P., Janik, M. J. & van Duin, A. C. T. A ReaxFF investigation of hydride formation in palladium nanoclusters via Monte Carlo and molecular dynamics simulations . J. Phys. Chem. C 118, 4967–4981 (2014).
Zou, C., Duin, A. T. & Sorescu, D. Theoretical investigation of hydrogen adsorption and dissociation on iron and iron carbide surfaces using the ReaxFF reactive force field method . Top. Catal. 55, 391–401 (2012).
Raju, M., van Duin, A. C. T. & Fichthorn, K. A. Mechanisms of oriented attachment of TiO2 nanocrystals in vacuum and humid environments: reactive molecular dynamics . Nano Lett. 14, 1836–1842 (2014).
Rahaman, O., van Duin, A. C. T., Goddard, W. A. & Doren, D. J. Development of a ReaxFF reactive force field for glycine and application to solvent effect and tautomerization . J. Phys. Chem. B 115, 249–261 (2011).
Hatzell, M. C. et al. Effect of strong acid functional groups on electrode rise potential in capacitive mixing by double layer expansion . Environ. Sci. Technol. 48, 14041–14048 (2014).
Senftle, T. P., van Duin, A. C. T. & Janik, M. J. Determining in situ phases of a nanoparticle catalyst via grand canonical Monte Carlo simulations with the ReaxFF potential . Catal. Commun. 52, 72–77 (2014).
Addou, R. et al. Influence of hydroxyls on Pd Atom mobility and clustering on rutile TiO2(011)-2×1 . ACS Nano 8, 6321–6333 (2014).
Spanjers, C. S. et al. Illuminating surface atoms in nanoclusters by differential X-ray absorption spectroscopy . Phys. Chem. Chem. Phys. 16, 26528–26538 (2014).
Raju, M., Ganesh, P., Kent, P. R. C. & van Duin, A. C. T. Reactive force field study of Li/C systems for electrical energy storage . J. Chem. Theor. Comput. 11, 2156–2166 (2015).
Islam, M. M. et al. ReaxFF molecular dynamics simulations on lithiated sulfur cathode materials . Phys. Chem. Chem. Phys. 17, 3383–3393 (2015).
Dereli, G. Stillinger-Weber type potentials in monte carlo simulation of amorphous silicon . Mol. Simulat. 8, 351–360 (1992).
Neyts, E. C. & Bogaerts, A. Numerical study of the size-dependent melting mechanisms of nickel nanoclusters . J. Phys. Chem. C 113, 2771–2776 (2009).
Voter, A. F. Parallel replica method for dynamics of infrequent events . Phys. Rev. B 57, R13985–R13988 (1998).
Joshi, K. L., Raman, S. & van Duin, A. C. T. Connectivity-based parallel replica dynamics for chemically reactive systems: from femtoseconds to microseconds . J. Phys. Chem. Lett. 4, 3792–3797 (2013).
Cheng, T., Jaramillo-Botero, A., Goddard, W. A. & Sun, H. Adaptive accelerated ReaxFF reactive dynamics with validation from simulating hydrogen combustion . J. Am. Chem. Soc. 136, 9434–9442 (2014).
Miron, R. A. & Fichthorn, K. A. Accelerated molecular dynamics with the bond-boost method . J. Chem. Phys. 119, 6210–6216 (2003).
Mortier, W. J., Ghosh, S. K. & Shankar, S. Electronegativity equalization method for the calculation of atomic charges in molecules . J. Am. Chem. Soc. 108, 4315–4320 (1986).
Chen, J. & Martinez, T. J. QTPIE: Charge transfer with polarization current equalization. A fluctuating charge model with correct asymptotics . Chem. Phys. Lett. 438, 315–320 (2007).
Lee Warren, G., Davis, J. E. & Patel, S. Origin and control of superlinear polarizability scaling in chemical potential equalization methods . J. Chem. Phys. 128, 144110 (2008).
Chelli, R., Procacci, P., Righini, R. & Califano, S. Electrical response in chemical potential equalization schemes . J. Chem. Phys. 111, 8569–8575 (1999).
Verstraelen, T. et al. Assessment of atomic charge models for gas-phase computations on polypeptides . J. Chem. Theor. Comput. 8, 661–676 (2012).
Verstraelen, T. et al. The significance of parameters in charge equilibration models . J. Chem. Theor. Comput. 7, 1750–1764 (2011).
Valentini, P., Schwartzentruber, T. E. & Cozmuta, I. Molecular dynamics simulation of O2 sticking on Pt(111) using the ab initio based ReaxFF reactive force field . J. Chem. Phys. 133, 084703/084701–084703/084709 (2010).
Verstraelen, T., Ayers, P. W., Van Speybroeck, V. & Waroquier, M. ACKS2: Atom-condensed Kohn-Sham DFT approximated to second order . J. Chem. Phys. 138, 074108 (2013).
Verstraelen, T., Vandenbrande, S. & Ayers, P. W. Direct computation of parameters for accurate polarizable force fields . J. Chem. Phys. 141, 194114 (2014).
Nistor, R. A., Polihronov, J. G., Müser, M. H. & Mosey, N. J. A generalization of the charge equilibration method for nonmetallic materials . J. Chem. Phys. 125, 094108 (2006).
Su, J. T. & Goddard, W. A. The dynamics of highly excited electronic systems: applications of the electron force field . J. Chem. Phys. 131, 244501 (2009).
Kale, S., Herzfeld, J., Dai, S. & Blank, M. Lewis-inspired representation of dissociable water in clusters and Grotthuss chains . J. Biol. Phys. 38, 49–59 (2012).
Rienstra-Kiracofe, J. C., Tschumper, G. S., Schaefer, H. F., Nandi, S. & Ellison, G. B. Atomic and molecular electron affinities: photoelectron experiments and theoretical computations . Chem. Rev. 102, 231–282 (2002).
Jordan, K. D. & Burrow, P. D. Studies of the temporary anion states of unsaturated hydrocarbons by electron transmission spectroscopy . Acc. Chem. Res. 11, 341–348 (1978).
Bochevarov, A. D. et al. Jaguar: A high-performance quantum chemistry software program with strengths in life and materials sciences . Int. J. Quantum. Chem. 113, 2110–2142 (2013).
van Duin, A. C. T., Baas, J. M. A. & van de Graaf, B. Delft molecular mechanics: A new approach to hydrocarbon force fields . J. Chem. Soc. 90, 2881–2895 (1994).
Iype, E., Hutter, M., Jansen, A. P. J., Nedea, S. V. & Rindt, C. C. M. Parameterization of a Reactive Force Field using a Monte Carlo Algorithm . J. Comput. Chem. 34, 1143–1154 (2013).
Jaramillo-Botero, A., Naserifar, S. & Goddard, W. A. General Multiobjective Force Field Optimization Framework, with Application to Reactive Force Fields for Silicon Carbide . J. Chem. Theor. Comput. 10, 1426–1439 (2014).
Larsson, H., Hartke, B. & van Duin, A. Global optimization of parameters in the reactive force field ReaxFF for SiOH . J. Comput. Chem. 34, 2178–2189 (2013).
Rice, B. M., Larentzos, J. P., Byrd, E. F. C. & Weingarten, N. S. Parameterizing complex reactive force fields using multiple objective evolutionary strategies (MOES): Part 2: transferability of ReaxFF models to C-H-N-O energetic materials . J. Theor. Comput. Chem. 11, 392–405 (2015).
Dittner, M., Muller, J., Aktulga, H. M. & Hartke, B. Efficient global optimization of reactive force-field parameters . J. Comput. Chem. 36, 1550–1561 (2015).
Plimpton, S. J. Fast parallel algorithms for short-range molecular dynamics . J. Comput. Phys. 117, 1–19 (1995).
Budzien, J., Thompson, A. P. & Zybin, S. V. Reactive molecular dynamics simulations of shock through a single crystal of pentaerythritol tetranitrate . J. Phys. Chem. B 113 (2009).
Zybin, S. V., Goddard, W. A. III, Xu, P., van Duin, A. C. T. & Thompson, A. P. Physical mechanism of anisotropic sensitivity in pentaerythritol tetranitrate from compressive-shear reaction dynamics simulations . Appl. Phys. Lett. 96, 081918/081911–081918/081913 (2010).
Li, Y., Kalia, R. J., Nakano, A., Nomura, K. & Vashishta, P. Multistage reaction pathways in detonating high explosives . Appl. Phys. Lett. 105, 204101–204103 (2014).
Nomura, K., Kalia, R. K., Nakano, A., Vashishta, P. & van Duin, A. C. T. Mechanochemistry of nanobubble collapse near silica in water . Appl. Phys. Lett. 101, 073108/073101–073108/073104 (2012).
Vedadi, M. et al. Structure and dynamics of shock-induced nanobubble collapse in water . Phys. Rev. Lett. 105, 014503/014501–014503/014504 (2010).
Chen, H.-P. et al. Embrittlement of metal by solute segregation-induced amorphization . Phys. Rev. Lett. 104, 155502/155501–155502/155504 (2010).
Nakano, A. et al. De novo ultrascale atomistic simulations on high-end parallel supercomputers . Int. J. High Perform. Comput. Appl. 22, 113–128 (2008).
Nomura, K. I. et al. Dynamic transition in the structure of an energetic crystal during chemical reactions at shock front prior to detonation . Phys. Rev. Lett. 99, 148303 (2007).
Nakano, A. et al. A divide-and-conquer/cellular-decomposition framework for million-to-billion atom simulations of chemical reactions . Comput. Mater. Sci. 38, 642–652 (2007).
Zheng, M., Li, X. & Guo, L. Algorithms of GPU-enabled reactive force field (ReaxFF) molecular dynamic . J. Mol. Graph. Model. 41, 1–11 (2013).
Aktulga, H. M., Knight, C., Coffman, P., Shan, T.-R. & Jiang, W. Optimizing the performance of reactive molecular dynamics simulations for multi-core architectures . IEEE Transac. Parallel Distrib. Syst. (2015) (in preparation).
Liu, B., Lusk, M. T., Ely, J. F., van Duin, A. C. T. & Goddard, W. A. Reactive molecular dynamics force field for the dissociation of light hydrocarbons on Ni(111) . Mol. Simulat. 34, 967–972 (2008).
Castro-Marcano, F., Kamat, A. M., Russo, M. F. Jr, van Duin, A. C. T. & Mathews, J. P. Combustion of an Illinois No. 6 coal char simulated using an atomistic char representation and the ReaxFF reactive force field . Combust. Flame 159, 1272–1285 (2012).
Vasenkov, A. et al. Reactive molecular dynamics study of Mo-based alloys under high-pressure, high-temperature conditions . J. Appl. Phys. 112, 013511 (2012).
Agrawalla, S. & van Duin, A. C. T. Development and application of a ReaxFF reactive force field for hydrogen combustion . J. Phys. Chem. A 115, 960–972 (2011).
Manzano, H. et al. Confined water dissociation in microporous defective silicates: mechanism, dipole distribution, and impact on substrate properties . J. Am. Chem. Soc. 134, 2208–2215 (2012).
TPS, MJJ, and ACTvD acknowledge funding from the National Science Foundation, grant CBET-1032979. ACTvD acknowledges support from the Fluid Interfaces Reactions, Structures and Transport (FIRST) center funded by the US Department of Energy, Office of Energy, Office of Basic Energy Sciences. ACTvD and MMI acknowledge support from a grant from the US Army Research Laboratory through the Collaborative Research Alliance (CRA) for Multi Scale Multidisciplinary Modeling of Electronic Materials (MSME). TV is a post-doctoral fellow of the Fund for Scientific Research-Flanders (FWO) and received additional funding from the Research Board of the Ghent University (BOF) and BELSPO in the frame of IAP/7/05.
Correspondence to Adri C T van Duin.

References: V. 
 V. 
 V. 
 V. 
 V. 
 V. 
 V.