Introduction
Many devices are built from materials with very different atomic structures and involve various types of bonding: covalent, ionic, metallic, and secondary bonding. The analysis of such systems at the atomic level is challenging. Electronic-structure methods, of which techniques based on density functional theory (DFT) are the most prevalent, generally provide high materials fidelity.Reference Sholl and Steckel1 However, because of the large computational load, DFT methods are limited to relatively small system sizes, typically a few hundred atoms. Advances in DFT methods, whose computational cost scales linearly with system size, are increasing this number. Currently, however, DFT methods are generally not able to address many important and topical issues in materials science, particularly those associated with complex microstructures.
Alternately, atomic-level simulation methods attempt to describe materials structures using empirical descriptions of the interactions that encapsulate the complex bonding effects of the valence electrons without explicitly describing the electrons themselves. The advantage of these methods is their high computational speed and their ability to probe large systems (millions or even billions of atoms) at finite temperatures and under external driving forces. Of course, the empirical approach usually involves sacrificing some materials fidelity; typically, however, with a given interatomic potential, it is possible to describe a range of phenomena with at least semi-quantitative precision. Because they describe very different bonding types, empirical potentials for metals, ceramics, and semiconductors generally have very different functional forms. Indeed, it has not generally been possible for empirical methods to describe structures in which more than one bonding type is present.
The reactive force field (ReaxFF)Reference van Duin, Dasgupta, Lorant and Goddard2 and charged optimized many-body (COMB)Reference Yu, Sinnott and Phillpot3, Reference Shan, Devine, Hawkins, Asthagiri, Phillpot and Sinnott4 potential methods are two independently developed platforms for describing complex bonding environments found in real materials. Although different in a large number of details, they are built around the same two fundamental concepts: self-consistent charge equilibration and bond order. A fundamental difference between the COMB and ReaxFF formalisms on the one hand, and most other force fields on the other, is that COMB and ReaxFF do not use fixed connectivity for the chemical bonds.
Both COMBReference Yu, Sinnott and Phillpot3, Reference Shan, Devine, Hawkins, Asthagiri, Phillpot and Sinnott4 and ReaxFFReference van Duin, Dasgupta, Lorant and Goddard2 were developed to help bridge the gap between quantum mechanical treatments of atomic forces and the non-reactive interatomic potentials traditionally used in atomic-level simulations (for example, the MM3Reference Allinger, Yuh and Lii5 and AMBERReference Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell and Kollman6 [assisted model building with energy refinement] force fields). Both are reactive force fields and include the effects of many-body interactions in a physically realistic, albeit empirical, manner. While most non-reactive empirical force fields are generally reliable near equilibrium with fixed bond orders and fixed charges, ReaxFF and COMB can describe reactive processes by allowing bond formation and bond breaking during molecular dynamics simulations. They can also be used to characterize the energetics of various reaction intermediates in addition to the reactants and products.Reference Chen, Kalia, Kaxiras, Lu, Nakano, Nomura, van Duin, Vashishta and Yuan7–Reference Nakano, Kalia, Nomura, Sharma, Vashishta, Shimojo, van Duin, Goddard, Biswas and Srivastava9 The force field parameters used in the potential terms are optimized against DFT and/or quantum chemical calculations and/or literature values. There is a price to be paid for this greater flexibility: these methods are significantly slower (approximately two orders of magnitude) than traditional, non-reactive force fields typically used to describe metallic, ionic, and covalent materials. Nevertheless, they are still much faster than DFT calculations, enabling parallelized applications to simulate systems larger than 10Reference Cornell, Cieplak, Bayly, Gould, Merz, Ferguson, Spellmeyer, Fox, Caldwell and Kollman6 atoms at nanosecond time scales.Reference Chen, Kalia, Kaxiras, Lu, Nakano, Nomura, van Duin, Vashishta and Yuan7–Reference Zybin, Goddard, Xu, van Duin and Thompson10 Consequently, COMB and ReaxFF can overcome many of the limitations inherent to nonreactive molecular dynamics (MD) simulation while retaining the desirable scalability. In this article, we provide an integrated overview of these methods. Their capabilities are illustrated, and their similarities and differences are discussed.
Methodology
In this section, we discuss the frameworks underpinning the ReaxFF and COMB methods. Each has bond order–dependent energy contributions as well as long-range non-bonded interactions. The charge equilibration schemes, by which the atoms determine their own charges, are very similar. The ReaxFFReference van Duin, Dasgupta, Lorant and Goddard2, Reference Chenoweth, van Duin and Goddard11 and COMBReference Yu, Sinnott and Phillpot3, Reference Shan, Devine, Hawkins, Asthagiri, Phillpot and Sinnott4 overall system energies are partitioned into several contributions:
The first four terms in each of the energy functions essentially describe the same fundamental physics and chemistry. In particular, the self-energy term, E self, captures the difference in energy of the various charge states of an atom, while the Coulomb term, E Coul, captures the electrostatic interactions. Both potentials also include a description of non-bonding interactions through the van der Waals energy, E vdW. Much of the complex physics and chemistry is captured in the bond-order term, E bond. The concept of bond order is quite simple: the strength of the interaction between two atoms depends not only on their relative positions, but also on the number and location of other atoms with which they interact. The bond-order terms are corrected for over coordination to allow the long-range bond order necessary for accurate reaction barriers.Reference van Duin, Dasgupta, Lorant and Goddard2, Reference Chenoweth, van Duin and Goddard11, Reference Chenoweth, van Duin, Persson, Cheng, Oxgaard and Goddard12 As an example of the ReaxFF reaction barrier performance, Figure 1 shows the pathway for conversion of propane molecules on a V4O10 cluster, indicating good performance for both reaction energies as well as barriers. ReaxFF also includes explicit energy contributions related to valence angle and torsion angle distortion (E valence and E torsion, respectively), energy contributions related to three-body and four-body conjugation (E conjugation), a weak hydrogen bond term (E H-bond) to properly describe X-H-X non-bonded interactions in polar molecules, energy terms related to the number of valence electrons around an atom (E lone-pair, E over, and E under), and several energy terms (E others) to correct for allene-type, terminal triple-bond, and C2 species. All the energy contributions in ReaxFF, excluding the van der Waals and Coulomb terms, are bond order–dependent and disappear upon bond dissociation.
Comparison between a reactive force field (ReaxFF) and quantum mechancis (QM) for the oxidative dehydrogenation conversion of two propane molecules on a V4O10 cluster. The structures that correspond to transition states have been labeled TS. C/H/O/V atoms are colored grey/white/red/blue, respectively. Reproduced with permission from Reference Reference Chenoweth, van Duin, Persson, Cheng, Oxgaard and Goddard12. ©2008, American Chemical Society.
The physics that is associated with the E angle, E over, E under, E torsion, and E conjugation terms in ReaxFF is embedded in the bond order term (E bond) in COMB, which is coupled with the charge-dependent pairwise repulsion, V R(Rij,qi,qj), and attraction, V A(Rij,qi,qj):
![$$E^{{\rm{bond}}} = {1 \over 2}\sum\limits_i {\sum\limits_{j \ne i} {\left[ {V^{\rm{R}} \left({R_{ij} ,q_i ,q_j } \right) - \left({{\rm{b}}^{{\rm{angle}}} + {\rm{b}}^{{\rm{coord}}} + {\rm{b}}^{{\rm{torsion}}} + {\rm{b}}^{{\rm{conjugation n}}} } \right)V^{\rm{A}} \left({R_{ij} ,q_i {\rm{,}}q_j } \right)} \right]} } $$](https://static.cambridge.org/binary/version/id/urn:cambridge.org:id:binary:20151130144030719-0897:S0883769412000954_eqn3.gif?pub-status=live)
Here Rij is the distance between atoms i and j of charge qi and qj. The various bond-order terms, b angle, etc., are complicated functions. E others in COMB is used to correct the specific energies associated with bond angles. The reader is referred to Chenoweth et al.Reference Chenoweth, van Duin, Persson, Cheng, Oxgaard and Goddard12 and Russo et al.Reference Russo and van Duin13 for the details of the ReaxFF energy functions and to Shan et al.Reference Shan, Devine, Hawkins, Asthagiri, Phillpot and Sinnott4 and Devine et al.Reference Devine, Shan, Cheng, McGaughey, Lee, Phillpot and Sinnott14 for details of the COMB energy functions.
The other crucial components of these potentials are the self-energy terms, E self, and Coulomb electrostatic interactions, E Coul. These self-energy terms account for the energy cost associated with the charge state of each atom, which is normally expressed as a Taylor expansion series with respect to its charge:Reference Chenoweth, van Duin and Goddard11, Reference Chenoweth, van Duin, Persson, Cheng, Oxgaard and Goddard12
where χi is defined as the atomic electronegativity, and ηi is associated with the atomic hardness or the self-Coulomb interaction. The Coulombic electrostatic terms include charge-charge (Eij qq) and nuclear charge-charge (Eij qZ) interactions between each pair of atoms i and j. The repulsion between nuclear charges of atoms i and j is normally assumed to be already included in the pairwise repulsion in E bond.
Two temperature model
While replacing explicit electrons with effective chemical bonding interactions makes large-scale molecular simulations feasible, there are properties that are neglected with this approach that are independent of how carefully a potential function is developed. Thermal conductivity of metals near room temperature, for example, is dominated by electronic degrees of freedom. To include these and related effects, including the ability to model lasers interacting with metals, Zhigilei and co-workers introduced a hybrid “two temperature” model in which electron dynamics are treated as a continuum that interacts with the atomic degrees of freedom (C. Schafer, H.M. Urbassek, L.V. Zhigilei, Phys. Rev. B 66, 115404 [2002]). A similar scheme was also developed by Brenner and co-workers that included Joule heating due to current flow (C.W. Padgett, D.W. Brenner, Molecular Simulation 31, 749 [2005]).
The charge equilibration method (EEM [electronegativity equalization method]Reference Mortier, Ghosh and Shankar15 or QEq [charge equilibration]Reference Rappe and Goddard16) used by ReaxFF and COMB, by which the atomic charge is determined, is based on the concepts of electronegativity equalization. Many years ago, in part to overcome the limitations of fixed charges, Streitz and MintmireReference Streitz and Mintmire17 employed the QEq method combined with an embedded-atom method potential in a study of α-alumina and aluminum. Zhou et al.Reference Zhou, Wadley, Filhol and Neurock18 modified this charge model to fix the charge fluctuations by imposing reasonable charge bounds on the metallic and oxide charges. This charge model has been applied to ion bombardment and oxidation of metal alloys.Reference Zhou, Wadley, Filhol and Neurock18
COMB makes use of Streitz’s and Mintmire’s charge density distribution functions to describe the charge of an atom. The Coulomb interactions Eij qq and Eij qZ are calculated through the Coulomb integral operator and the shielded nuclear charge.Reference Streitz and Mintmire17 Since there is no nuclear charge term, the electrostatic interactions in ReaxFF are expressed as a shielded Coulomb potential between charges associated with atoms i and j. In addition to the local charge terms in Equation 4, ReaxFF takes into account the shielding between overlapping electron clouds at short distances. Furthermore, it employs a seventh order tapering function to screen long-range Coulomb interactions. As demonstrated by de Vos Burchart et al.,Reference de Vos Burchart, Verheij, van Bekkum and van de Graaf19, Reference de Vos Burchart20 this tapering treatment provides a close approximation to the formal Ewald sum calculationReference Allen and Tildesley21 of the total Coulombic energy for both crystalline and liquid systems. The local environmental dependence of the electronegativity is reflected in three places in COMB: E self, E Coul, and E bond. The self-energy term in COMB includes terms to adjust the atomic electronegativity and self-Coulomb changes of an atom embedded in a molecule or a crystalline structure. The overall integration of the Coulomb integral over 1s Slater-type orbitals and the shielded nuclear charge integral that is used for the Eij qq and Eij qZ in COMB also results in tapering to account for the electronegativity variation as a function of the local environments of an atom. As indicated in Equation 3, the short-range bond energy terms, V R(Rij, qi, qj) and V A(Rij, qi, qj), also include charge and spatial dependence. With these spatial dependences, COMB is capable of accurately describing the charge transfer in various local charge environments. The spatial dependences in E self and E bond were not included in the initial formalisms proposed by Streitz and Mintmire.Reference Streitz and Mintmire17 As a result, their formalism cannot generate physically reasonable charge distributions in some specific situations.Reference Thomas, Marks and Begg22 The enhancements introduced in COMB enable dynamic charge transfer to be appropriately captured in molecules, interfaces, and crystalline structures.
ReaxFF and COMB share many common features and are conceptually very similar. The intent in both COMB and ReaxFF is that the same element-specific parameters as well as the bond-type specific parameters will be used for all different bond types in a system. Thus the parameter sets in COMB and ReaxFF should be transferable to new systems.
In the development of the ReaxFF framework, specialized force fields have been produced on occasion, which deviate from the completely transferable main branches, thus enabling focused force field development without the burden of the large training sets associated with the main branches of development of the formalism. Examples of these specialized force fields include aminoborane,Reference Weismiller, van Duin, Lee and Yetter23 H2/O2,Reference Agrawalla and van Duin24 and the nitramine/peroxide ReaxFF descriptions.Reference Zhang, Zybin, van Duin and Goddard25–Reference Strachan, van Duin, Chakraborty, Dasgupta and Goddard30 There are currently two main, inter-transferable branches of ReaxFF parameter sets: one based on the ReaxFF water description, which has been extended to proteins,Reference Rahaman, van Duin, Goddard and Doren31, Reference Rahaman, van Duin, Bryantsev, Mueller, Solares, Goddard and Doren32 phosphates,Reference Zhu, Janetzko, Zhang, van Duin, Goddard and Salahub33, Reference Abolfath, van Duin, Biswas and Brabec34 and various metal oxides,Reference Russo, Li, Mench and van Duin35–Reference Aryanpour, van Duin and Kubicki40 and one based on the C/H/O combustion parameters developed by Chenoweth et al.Reference Chenoweth, van Duin and Goddard11 and applied to a range of combustion simulations,Reference Srinivasan and van Duin41–Reference Chenoweth, van Duin, Dasgupta and Goddard47 which has been extended to Mo/V/Bi-oxide chemistry,Reference Chenoweth, van Duin, Persson, Cheng, Oxgaard and Goddard12, Reference Goddard, Mueller, Chenoweth and van Duin48–Reference Goddard, Chenoweth, Pudar, van Duin and Cheng50 metal surface chemistry,Reference Neyts, van Duin and Bogaerts51–Reference Mueller, van Duin and Goddard56 and high-nitrogen compounds.Reference Kamat, van Duin and Yakovlev57 The COMB potentials were parameterized for Si/SiO2,Reference van Duin, Dasgupta, Lorant and Goddard2, Reference Shan, Devine, Hawkins, Asthagiri, Phillpot and Sinnott4 Hf/HfO2,Reference Shan, Devine, Kemper, Sinnott and Phillpot58 Cu/SiO2,Reference Shan, Devine, Sinnott and Phillpot59 and Cu/Cu2O,Reference Devine, Shan, Cheng, McGaughey, Lee, Phillpot and Sinnott14 and are currently being extended to a variety of systems, including Zn/ZnO/Cu/Cu2O, Zr/ZrO2, U/UO2, Al/AlN/Al2O3, and Ti/TiN/TiO2.
Applications of ReaxFF and COMB
As described previously, the COMB and ReaxFF potentials provide a flexible and transferable framework to simulate a wide range of chemical environments and reactions. This is illustrated here with a representative selection of applications. Silica is one of the most widely used materials in the design of electronic devices. While a number of empirical force fields have been developed for silica and related structures, their transferability is not wide enough to provide an adequate description of the physical and chemical properties in all its relevant structural forms, including crystalline, amorphous, liquid, bulk, and clusters. The description of the interactions of molecular systems with condensed phases is illustrated by water interacting with silica and with the interaction of organic molecules with metal surfaces and metal nanoparticles.
Phase order for silica
Silica is an ideal system in which to demonstrate the capabilities of these methods because of the wide range of bond lengths and bond angles found in its many polymorphs. The ReaxFF and COMB potentials for Si/SiO2 were developed with the aim of providing a consistent description for crystalline Si and all forms of silicon oxides. The COMB parameters for Si, which have been used in numerous studies, are the same as those in the original Tersoff article.Reference Tersoff60 The ReaxFF Si potentialReference van Duin, Strachan, Stewman, Zhang, Xu and Goddard61 closely matches quantum chemistry data and successfully reproduces the high-pressure transition from a four-coordinated tetragonal diamond phase to a six-coordinated β-Si phase. Figure 2 shows the volume versus energy equation of state for three relevant Si-phases, including the zero-pressure ground state (alpha-diamond), the six-coordinated beta-phase, and the symmetric six-coordinate simple cubic phase. COMB reproduces the correct phase order of the crystalline phases of SiO2 and a good description of the amorphous phase.Reference van Duin, Dasgupta, Lorant and Goddard2 ReaxFF accurately reproduces the heat of formation determined from DFT and experiments for a wide range of SiO2 phases, including four- and six-coordinated (stishovite) structures, as well as low-density zeolite structures.Reference van Duin, Strachan, Stewman, Zhang, Xu and Goddard61 For α-quartz, ReaxFF assigns charges of +1.346 to Si atoms and –0.673 to O atoms, which are very similar to the QEq charge distribution (+1.3180 for Si and –0.6590 for O) and closely matches the Mulliken charges from small basis set calculations (typically 6–31G*).Reference Mueller, van Duin and Goddard55 COMB assigns charges of +2.92 to Si atoms and –1.46 to O atoms, which are fit to the Bader charge analysis.Reference Henkelman, Arnaldsson and Jónsson62 As such, ReaxFF describes silica as a significantly covalent material, while COMB provides a description more in terms of formal ionic charges. The original ReaxFF Si/O descriptionReference van Duin, Strachan, Stewman, Zhang, Xu and Goddard61 of tridymite stated it was the global energy minimum, since the DFT data (the “training” set) used in fitting the parameters of the potential described tridymite as such. It has since been re-formulated by Fogarty et al.Reference Fogarty, Aktulga, Grama, van Duin and Pandit39 using the same condensed-phase training set, to be transferable with development of the ReaxFF water branch.
Equation of state (compression and expansion) for three phases of Si: α is tetrahedral (diamond), β is the six-coordinate, and cubic is a simple cubic six-coordinate structure, calculated using the (a) reactive force field (ReaxFF) and (b) quantum-chemical (QC) methods. Reproduced with permission from Reference Reference Mueller, van Duin and Goddard55. ©2003, American Chemical Society.
Covalent-ionic bonding at Si/SiO2 interfaces
While reproducing the correct phase order in SiO2 is a non-trivial demonstration of their capabilities, ReaxFF and COMB were specifically designed to describe systems with mixed types of atomic bonding. Thus both have been used to model various Si/SiO2 interfaces, at which covalently and ionically bonded Si atoms meet and interact. Such Si/SiO2 interfaces are prevalent in silicon-based materials. The COMB potential was used to model a Si (001)/α-quartz (010) interface.Reference Shan, Devine, Hawkins, Asthagiri, Phillpot and Sinnott4 The initial structure, constructed using simple crystallography, contained Si coordination defects: half of the Si atoms were two-fold coordinated, while the rest had the correct four-fold coordination. Upon annealing at high temperature, the structure evolved to a more physically realistic structure, in which all of the Si atoms became four-fold coordinated.
The ReaxFF method was used to model a slightly different problem:Reference van Duin, Strachan, Stewman, Zhang, Xu and Goddard61 the Si (001)/amorphous silica interface. An α-quartz slab was sandwiched within a Si slab and heated to melting. The simulation box was then cooled and equilibrated, forming the interface. After annealing, analysis of the Si–O–Si bond angles showed that most were between 135°–170°. However, about 10% of Si–O–Si bond angles were close to 90°, and some of the O atoms at the interface were three-fold coordinated. The distribution of Si–O–Si bond angles obtained using ReaxFF and COMB were similar and consistent with experiment: COMB predicted the distribution of Si–O–Si bond angles to be 110°–160° (135.4° average), while experiments yielded 130°–160°, with an average of 144°.Reference Coombs, Denatale, Hood, McElfresh, Wortman and Shackelford63
Chemisorption of oxygen on silicon surfaces
As discussed previously, one of the key features of these approaches is their ability to describe bonding and unbonding processes. Building on the expertise developed in describing the structural aspects of the Si-O system, this ability to describe chemistry has been explored using a very simple case: that of the adsorption of oxygen on Si, a first step in the oxidation of Si. In particular, the Si (001) {2×1} reconstructed surface and its reaction with oxygen has been studied with both ReaxFF and COMB. Three different adsorption sites were studied for a 0.5 monolayer surface coverage. COMB matched the energy trend in adsorption energies from DFT calculations within 6%. In a complementary study, the impact of atomic and molecular oxygen at various energies with the Si surface was studied using ReaxFF.Reference Neyts, Khalilov, Portois and van Duin64 It was found that the extent of penetration of O2 and O atoms (Figure 3) into the silicon surface depends on the impact energy, from chemisorption on the top layer at low energy, to the knock-on displacement of Si at high energy. Because the momentum of molecular oxygen is distributed into two oxygen atoms after dissociation due to impact, penetration of these oxygen atoms into the lattice is weaker than that of hyperthermal atomic oxygen. These simulations indicated that at low impact energies (less than 1.0 eV), the oxygen atoms remain confined to the surface Si layers, and the damage to the bulk is very limited and purely induced by the chemisorption process. At higher energies (higher than 1 eV), O atoms can penetrate the crystal to a depth of up to 20 Å. At these higher energies, the damage to the Si lattice is much more pronounced, with primary and secondary knock-on displacements of Si atoms.
A ReaxFF (reactive force field) simulation on the mechanism of defect creation in a silicon slab by oxygen atom chemisorption at low impact energy below 10 eV. The red atom is oxygen, all other atoms (beige, blue, and green) are silicon atoms. The green and blue atom coloring enables visualization of the silicon atom displacement due to O-chemisorption. (a) Initial, pre-impact configuration; (b) post-impact configuration, with re-configured Si-slab and oxygen incorporated in the Si-subsurface. Reproduced with permission from Reference Reference Neyts, Khalilov, Portois and van Duin64. ©2011, American Chemical Society.
Metal–ionic bonding at Cu/SiO2 interfaces
The interface between Cu and SiO2 is prevalent in microelectronics, particularly in the context of interconnects, and is a prototypical metal/ceramic interface. Such interfaces have traditionally been very difficult to model due to the differences in bonding types on the two sides. To assess the capabilities of the COMB potential, a comparison between three different Cu/SiO2 interfaces was performed: the Cu (001)/α-cristobalite (001), Cu (001)/α-quartz (010), and Cu (001)/amorphous silica interfaces.Reference Shan, Devine, Sinnott and Phillpot59 The number of oxygen atoms at the interface is expected to depend on the oxygen partial pressure; therefore, three different oxygen terminations were explored. The three resulting interfaces for Cu (001)/α-cristobalite (001) are shown in Figure 4. The DFT and MD simulations with the COMB potentials give qualitatively and quantitatively similar results.Reference Shan, Devine, Sinnott and Phillpot59 For all the interfaces studied, the adhesion energy between the interfaces depends on the oxygen content, with oxygen-rich interfaces having larger adhesion energies than oxygen-lean interfaces. Contact between Cu and amorphous silica shows that oxygen atoms play an important role in determining the bonding of the interface and the formation of a thin Cu2O oxide layer. The Cu–O bonds have a significant role in providing stronger adhesion between copper and silica substrates, in part due to enhanced charge transfer across the interface.
Minimized structures of the fcc Cu (001)/α-cristobalite (001) interfaces, OO terminated, O terminated, and Si terminated, left to right, with molecular dynamics simulations using charged optimized many-body potentials viewing along (a) [100] and (b) [010] directions. The Cu layer is on the top, while the SiO2 layer is at the bottom. Atoms are color coded by the charge values in units of the magnitude of the electron charge, |e|, indicated by the color bar (navy is negatively charged, cyan is charge neutral, and red is positively charged). Reproduced with permission from Reference Reference Shan, Devine, Sinnott and Phillpot59. ©2011, American Physical Society.
Hafnia-silicon interfaces in high-κ structures
The ever-decreasing feature dimensions of microelectronic structures required ever thinner thicknesses of the dielectric. As the thickness of the conventional silica layer decreased below 2 nm, the leakage current increased significantly. Thus, in the 45 nm technology introduced by Intel in 2007, the silica gate oxide was replaced by hafnia.65 Hafnia, HfO2, has a significantly higher dielectric constant than silica (20–25 versus 4),Reference Rignanese66, Reference Zhao and Vanderbilt67 thereby allowing thicker dielectric layers and reducing the leakage current. In addition, the high thermal stability of hafnia in contact with the semiconductor and the metal gate is another appealing factor.Reference Gusev, Cabral, Copel, D’Emic and Gribelyuk68 It was found that the COMB formalism reproduces the structure and properties of metallic Hf quite well.Reference Shan, Devine, Kemper, Sinnott and Phillpot58 The predicted c/a ratio of the hcp phase is comparable to the experimental value, indicating the capability of COMB to capture directional bonds in non-ideal hcp structures. The elastic constants, stacking fault energy, and coefficient of linear thermal expansion are in good agreement with experimental and/or ab initio values. Most particularly, it reproduced the correct phase order: hcp structure in the ground state and fcc and bcc structures having somewhat higher energies. In describing the oxide phase, COMB was parameterized to yield the correct phase order. Most importantly, a monoclinic HfO2 structure has the lowest energy. Moreover, on heating the monoclinic structure, the system transformed to the tetragonal phase at around 2300 K, in agreement with experiments.Reference Weast69 It was also found using COMB that the monoclinic-to-tetragonal phase transition nucleated from the free surface, which corresponds well with the notion that defects and surfaces are common nucleation sites. The phase transition is particularly important because DFT predicts, and COMB reproduces, three other phases (α-PbO2 structure, orthorhombic I, and rutile) as having lower energy than the tetragonal and cubic phases.Reference Shan, Devine, Kemper, Sinnott and Phillpot58 However, they do not appear in the experimental heating curves or in the heating curves from simulation using the COMB potentials, indicating significantly higher phase transition barriers from the monoclinic phase to this phase; this is consistent with them not being present in the ambient pressure phase diagram. In describing the defect energetics of monoclinic hafnia, the oxygen interstitial is predicted to be the most stable point defect, in agreement with ab initio calculations.Reference Foster, Gejo, Shluger and Nieminen70 With these strong indications of the high materials fidelity of the Hf/HfO2 COMB parameterization, the structure of the Si/HfO2 interface was characterized. Hafnia deposited on Si via atomic layer deposition, pulsed laser deposition, or chemical vapor deposition methods is typically amorphous as depositedReference Ritala, Leskela, Niinisto, Prohaska, Friedbacher and Grasserbauer71–Reference Fang, Zhang, Wang, Modreanu, O’Sullivan, Hurley, Leedham, Hywel, Audier, Jimenez, Senateur and Boyd73 and polycrystalline after annealing.Reference Luo, Demkov, Triyoso, Fejes, Gregory and Zollner74 A recent study showed that the cubic phase of hafnia is a promising polymorph for the gate oxide because of its higher dielectric constant and because the dielectric constant is isotropic.Reference Ramana, Noor-A-Alam, Gengler and Jones75Figure 5 shows the structure of a coherent cubic-HfO2/Si interface equilibrated with the COMB formalism. As characterized by their yellow and green colors, the Si atoms at the interface take on a considerable charge, indicating that they have oxidized. These Si atoms also form four-fold coordination with O atoms, characteristic of ionic bonding at the interface as well as the presence of a silica layer. The presence of such a silica layer is undesirable from the point of view of the high-κ behavior: capacitances in series are added inversely, thereby lowering the overall dielectric constant of the gate stack. However, the strong ionic bonding contributes to a high work of adhesion at the interface; such strong adhesion should aid in template growth of HfO2 on Si and in the minimization of harmful interfacial traps.
Structure of the hafnium oxide/silicon interface. The atoms are colored by charge. The larger atoms (shades of red) are the hafnium atoms, the small atoms (blue) are the oxygen atoms, and the medium-sized atoms (light green, yellow, and blue) are the silicon atoms. Reproduced with permission from Reference Reference Shan, Devine, Kemper, Sinnott and Phillpot58. ©2011, American Physical Society.
Metal surface catalysis: Complex hydrocarbon conversion on Ni-nanoparticles
Metal surfaces and metal nanoparticles are key components of many successful heterogeneous catalysts. These catalysts often have highly complex surface structures, including reactant atoms diffusing among surface, subsurface, and bulk, as well as various surface facets, step sites, and grain boundaries. In order to understand, and optimize, the structure of these catalysts, a detailed understanding of the dynamics of these surfaces is of vital importance. The ReaxFF potential has been used to study these dynamics.
In order to study the long-time (nanosecond) structure and catalytic activity of Ni-surfaces and nanoparticles, Mueller et al.Reference Mueller, van Duin and Goddard54–Reference Mueller, van Duin and Goddard56 developed a ReaxFF potential for hydrocarbon binding and conversion on Ni-metal surfaces. This potential was trained against a substantial DFT-based training set, containing CxHy binding energies to Ni[111] and Ni[100] surfaces (x = 0–2; y = 0–6), hydrocarbon dehydrogenation barriers, and equations of state for Ni-carbides, Ni-metal, and Ni-hydrides.
Figure 6 shows an application of this Ni/C/H description. In this simulation, a Ni-nanoparticle was exposed to a propene gas phase and gradually heated from a temperature of T = 500 K to T = 2500 K. At low temperatures, we observed the adsorption of propene molecules, making a relatively strong bond through the C=C double bond interacting with Ni-surface atoms. At around 1000 K, the first dehydrogenation steps are observed, leading to the formation of a wide range of C–C cracking and dehydrogenation steps and H2 release from the surface, eventually resulting in the formation of a Ni-carbide nanoparticle, covered with a complex mixture of hydrocarbon fragments.
A ReaxFF (reactive force field) (a) initial and (b) final structures for propene adsorption and decomposition on a Ni468-nanoparticle during heating from 500 K to 2500 K, resulting in H2 release and formation of surface and subsurface CHx-species. C/H/Ni atoms are colored green/red/blue, respectively. Reproduced with permission from Reference Reference Mueller, van Duin and Goddard55. ©2010, American Chemical Society.
This same ReaxFF Ni/C/H potential was subsequently used to study the formation of carbon nanotubes on Ni-nanoparticles and was recently successfully applied in the elucidation of the key steps in the formation of chiral carbon nanostructures.Reference Neyts, van Duin and Bogaerts51, Reference Neyts, Shibuta, van Duin and Bogaerts53
Water/silica and DMMP/silica interfaces
The interaction of silica with water is a prototype for a wide range of geological, biological, and technological issues. Many ab initio calculations of the reaction pathway for breaking the surface Si–O bond have been performed, yielding estimates of the activation energy for dissolution.Reference Nangia and Garrison76–Reference Xiao and Lasaga78 Moreover, many reactive classical molecular dynamics studies have been reported at silica-water and quartz-water interfacesReference Argyris, Cole and Striolo79–Reference Vanbeest, Kramer and van Santen83 and have shown proton transfer at the water/silica interface.Reference Lockwood and Garofalini84–Reference Mahadevan and Garofalini86 To describe the proton-transfer reactions at the water/amorphous silica interface, Fogarty et al.Reference Fogarty, Aktulga, Grama, van Duin and Pandit39 modified ReaxFF by re-optimizing the Si/O/H bond and angle parameters to specifically fine-tune the interfacial structure and to more accurately reproduce the properties of water. This ReaxFF description successfully reproduced the surface chemistry, most importantly the structure of water and amorphous silica at the interface. In addition, the analysis of orientation of water molecules at the interface indicated that the water molecules are strongly oriented, causing a large charge distribution along the interface and an electric potential difference between the bulk water and the interface. Proton transfer into the bulk silica was also observed. ReaxFF has also been used to characterize the interactions of organophosphates with silica surfaces.Reference Quenneville, Taylor and van Duin87
Outlook
The ultimate success of a potential description relies on it being widely adopted, thereby allowing its capabilities and limitations with regard to specific materials systems to be assessed. For example, it can be argued that silicon as modeled by the Stillinger-Weber potentialReference Stillnger and Weber88 has become so ubiquitous that it can be regarded as a virtual material with well-characterized physical properties. The variable-charged potentials used here are quite complicated to implement. Therefore, to promote wide adoption of these approaches and to enable the community as a whole to develop them further, both have been implemented into the LAMMPS (large-scale atomic/molecular massively parallel simulation) program.89, Reference Plimpton90 As discussed in the article by Plimpton and Thompson in this issue, the current implementations of both of these methods are much slower than conventional potentials. Some of this slowness is intrinsic, a result of the computational cost of the charge equilibration. However, neither is fully optimized from an algorithmic view yet; moreover, there may also be opportunities for making approximations that reflect the underlying physics and further improve computational efficiency.
In conclusion, the COMB (charge optimized many-body) and ReaxFF (reactive force field) methodologies have already demonstrated significant power to describe complex physical and chemical behavior in a wide range of materials. They encapsulate in slightly different forms two essential aspects of chemical bonding: the ability of an atom to autonomously determine its charge state by interrogating its environment, and the modulation of the interaction of an atom with a neighbor, also by identifying other proximate atoms and their positions. Many of the specific differences in the current formalisms arise from the different materials problems to which they have been applied. Moreover, as they are asked to describe different kinds of phenomena, it can be expected that additional physical and chemical processes will be encoded in additional terms. We can thus expect both of them to evolve, sometimes in a convergent manner, sometimes in a divergent manner, and in unexpected ways. Moreover, it is quite possible that at some point, they will both be subsumed in a more powerful, more flexible formalism or that the concepts underlying these methods will provide the foundation for the next-generation of transferable, efficient empirical bond order potentials.
Acknowledgments
The development of COMB (T.-R.S., S.R.P., and S.B.S.) has been supported by the National Science Foundation under DMR-0426870 and DMR-1005779. The work of T.L. was supported as part of the Center for Atomic Level Catalyst Design, an Energy Frontier Research Center funded by the US Department of Energy, Office of Science, Office of Basic Energy Sciences under Award Number DE-SC0001058, and by NSF CHE-0809376. The work of M.J.N. was supported by the Consortium for Advanced Simulation of Light Water Reactors, an Energy Innovation Hub for Modeling and Simulation of Nuclear Reactors under US Department of Energy Contract No. DE-AC05–000R22725. The development of ReaxFF (A.C.T.v.D. and Y.K.S.) was originally funded by the British Royal Society and is currently funded by NSF Grants OCT 1047857 and CBET 1032979, NETL/RUA contract 662.884.001, the FIRST center (an EFRC/DoE funded center), AFOSR Grant FA9550–10–1-0563, DoE grant DE-FE0005867, and AFRL/SBIR contract FA8650–11-C-2185.
