There are of course whole classes of problems for which the lack of ability to model bond breaking and forming is a show-stopper. Though in principle quantum calcuations could be done in these cases, the number of atoms typically included in a system usually precludes this. Instead, a lot of work has gone into developing reactive potentials. The three main classes we consider here are the embedded-atom method, bond-order potentials and ReaxFF.
The embedded-atom method (EAM) was originally developed to model metals. [28] The basic idea of EAM is that atoms interact in a pairwise manner with the nearest neighbors but they also interact with a global field of electron density that is explicitly many-body in nature. For a system of \(N\) atoms, the EAM potential is
where
\begin{equation} \rho _i = \sum _{j\ne i}g\left (r_{ij}\right ), \end{equation} \begin{equation} g\left (r\right ) = \exp {\left (-\beta r\right )}, \end{equation} \begin{equation}\label {eq:EAM:phi} \phi \left (r\right ) = V_2\left (r\right ) - 2F\left [g\left (r\right )\right ], \end{equation} \begin{equation} V_2\left (r\right ) = \begin {cases} V_{\rm ZBL}\left (r\right ), & r<r_1, \\ \alpha _0+\alpha _1r+\alpha _2r^2+\alpha _3r^3, & r_1\le r<r_2, \\ -2c\exp {\left (-{{\beta }\over {2}}r\right )}+\Phi _0\exp {\left (-\alpha r\right )}, & r\ge r_2, \end {cases} \end{equation} \begin{equation}\label {eq:ZBL} V_{\rm ZBL}\left (r\right ) = {{Z_1Z_2e^2}\over {4\pi \epsilon _0r}}\sum _{i=1}^{4}c_i\exp \left (-d_i{{r}\over {a}}\right ) \end{equation}
The second term in Eq. 227 accounts for the fact that the electron density into which an atom is embedded does not include electrons from that atom itself. Eq. 229 is the Ziegler-Biersack-Littmark (ZBL) screened nuclear repulsion potential used for modeling high-energy collisions between atoms. The three branches that make up \(V_2\) produce a spline-connection between the ZBL and a Morse-like attractive tail.
Among the many systems simulated using EAM is the sputtering of copper. [29]
Bond-order potentials aim to capture the effect of the nearest-neighbor environment on the behavior of any bond. Such potentials began with the work of Abell [30]. Here I only very lightly gloss over this very deep field of research. In the bond-order formalism, the total potential energy due to covalent bonds is: \begin{equation} \label {eqn:U} U = \sum _i{\sum _{j>i} {\phi _{ij}}}, \end{equation} where the bond energy \(\phi _{ij}\) between atoms \(i\) and \(j\) has repulsive and attractive components: \begin{equation} \label {eqn:phi} \phi _{ij} = V_R(r_{ij}) - \overline {b}_{ij}V_A(r_{ij}) , \end{equation} where \(V_R\) and \(V_A\) are Morse-type pair potentials: \begin{equation} \label {eqn:Vr} V_R(r_{ij}) = f_{ij}(r_{ij})A_{ij}\exp (-\lambda _{ij}r_{ij}) {\rm \hspace {5mm} and} \end{equation} \begin{equation} \label {eqn:Va} V_A(r_{ij}) = f_{ij}(r_{ij})B_{ij}\exp (-\mu _{ij}r_{ij}) \end{equation} Here, \(r_{ij}\) is the scalar separation between atoms \(i\) and \(j\). \(A_{ij}\), \(B_{ij}\), \(\lambda _{ij}\), and \(\mu _{ij}\) are all parameters specific to the two elements participating in the bond. (I use the following \(ij\)-subscript convention: when \(ij\) appears on a variable, \(i\) and \(j\) refer to the individual atom indices; when \(ij\) appears on a parameter or function, \(i\) and \(j\) are specific only to the elements of atoms \(i\) and \(j\).) The cutoff function \(f_{ij}\) decays smoothly from 1 at some “inner” radius to 0 at some “outer” radius.
The bond order \(\overline {b_{ij}}\) models all of the many-body chemistry: \begin{equation} \label {eqn:bijbar} \overline {b_{ij}} = {1\over 2} \left [b_{ij}+b_{ji}+\mbox {corrections}\right ], \end{equation} where \(b_{ij}\) is the contribution of atoms that neighbor atom \(i\) to the bond order of the \(ij\) bond. These involve complicated 3- and 4-body interactions. The corrections arise from the need to expand set of thermochemical data which the potentials are fit against.
Bond-order potentials first applied to metals, but expanded into silicon and hydrocarbons with the work of and Tersoff [31]and Brenner [32], respectively, producing what is called the “reactive empirical bond order potential (REBO). The more recent adaptive intermolecular reactive empirical bond-order (AIREBO) potential combines REBO with Lennard-Jones interactions and specific torsional potentials for better modeling of hydrocarbon chains. [33]
Polarizable versions of reactive potentials have also been developed. The charge-optimized many-body (COMB) potential is an extension of REBO in which the charge on each atom is allowed to change according to energies dictated by input parameters such as atom electronegativity. [34] Adding oxygen to the AIREBO hydrocarbon potential also necessitated including polarizability, leading to the qAIREBO potential. [35]
ReaxFF was introduced in 2001 as a new type of reactive force field for hydrocarbons, and its formalism has now greatly expanded to describe reactive interatomic chemistries for a wide swath of the periodic table [36, 37]. ReaxFF adopts the co-called “central-force” formalism wherein all pairs of atoms interact and there are no switching functions. ReaxFF breaks down the potential into seven distinct additive contributions:
The role of the “overcoordiation penalty” \(U_{\rm over}\) is to monitor atom valencies and penalize deviations from an ideal. This keeps hydrogens from binding to more than one partner, and carbons more the four. The bond, angle, and torsion terms mimic those of class-I potentials but essentially determine parameters on the fly based on geometries. The Coulomb and VdW terms represent electrostatics and dispersion/excluded-volume interactions, while atom charges are forced to equilibrate given atom electronegativities. The “specific” term is the “secret sauce” where a lot of things are put to specialize ReaxFF for a particular system.
ReaxFF has been used for a large variety of reactive systems to date. It is fairly expensive to run and requires a lot of tuning when applying it to systems it has never been applied to before.