All-atom Potential Energy Functions

8 All-atom Potential Energy Functions

The modeling of molecular structure and inter/intra-molecular interactions is the job of the potential energy function $\mathscr {U}(\rb ^N)$. Modern potential energy functions are actively maintained and carefully curated and optimized by several devoted groups around the world, and for the most part they are made freely available to the research community. In this section, we’ll consider a few of the more popular potentials out there. There are many good reviews out there about all-atom potential energy functions; most of what I present here is adapted from the recent review by Harrison et al. [17]

First, a note on terminology. A potential energy function, or just “potential”, is conceptually just a function that can compute potential energy from all atomic positions. Pairwise Lennard-Jones is an example we have used extensively already. When used specifically by MD simulations, potentials are often referred to as “force-fields”, since it is really their gradients that are used in MD; potential energy itself is more of a diagnostic in most standard MD simulations (though it is extremely important in some advanced free-energy methods). Because of this, potentials that are used in MD simulations have to be differentiable, for the most part. (There are a few MD methods that can use so-called “discontinuous forces” but no large-scale MD codes can.)

Second, a note on physical reality vs what potentials actually model. Interactions among atoms are quantum-mechanical in nature. Modeling them accurately involves in-depth quantum chemical calculations that can be very expensive. Potentials used in most molecular simulations are empirical functions that generally only very roughly approximate true interatomic interactions. Nuclei are treated as point masses and electrons aren’t treated at all. Because of this simplification, empirical potentials always have parameters that must be tuned against more accurate quantum chemical calculations and experimental observations. Parameter tuning in potentials leads to specialization that reflects where such potentials are most in demand. This will become evident as we start discussing important examples.

8.1 Class-I Potentials

Class-I potentials break down interatomic interactions between bonded and non-bonded interactions. Bonded interactions include bonds, valence angles, and both proper and improper dihedrals (Fig. 29). Non-bonded interactions include electrostatic and Lennard-Jones interactions. Class-I potentials apply for systems of fixed bonded topology; i.e., bonds between atoms are permanent, never breaking or forming. A simulation using a Class-I potential is therefore unable to model “chemistry”. Examples of Class-I potentials include CHARMM [18, 19, 20], AMBER [21], Gromos [22], OPLS [23, 24], and TraPPE [25].

Figure 29: Schematic representation of (A) the $i$-$j$ bond vector $\rb _{ij}$, (B) the $j$-$i$-$k$ angle $\theta _{jik}$, (C) the $i$-$j$-$k$-$l$ torsion $\phi _{ijkl}$ around the $j$-$k$ bond, and (D) the $i$-$j$-$k$-$l$ improper $\phi _{ijkl}$, comprising bonds $i$-$j$, $i$-$k$, and $i$-$l$ (solid), with vectors $\rb _{ij}$, $\rb _{jk}$, $\rb _{kl}$ (dashed) defining the dihedral angle. In (C), the angle $\phi _{ijkl}$ is acute and positive by IUPAC/IUB convention; if we sight down the $j$-$k$ bond with $j$ nearer to the eye, the $i$-$j$ bond must rotate in a clockwise direction by an angle less than $\pi $ to eclipse the $k$-$l$ bond.

The general form of a Class-I potential is

\begin{align} \label {eq:class-i-potential} \mathscr {U} = & \sum _{\rm bonds} K_b(l-l_0)^2 + \sum _{\rm angles} K_\theta (\theta -\theta _0)^2\\ & + \sum _{\rm dihedrals} K_{\phi ,n}\left [1+\cos (n\phi - \delta _n)\right ] \\ & + \sum _{\rm impropers} K_\phi (\phi -\phi _0)^2 \\ & + \sum _{i<j} \left \{ \frac {q_iq_j}{r_{ij}} + 4\epsilon _{ij}\left [\left (\frac {\sigma }{r_{ij}}\right )^{12}-\left (\frac {\sigma }{r_{ij}}\right )^{6}\right ]\right \} \end{align}

The first four terms are the “bonded” potential. Each bond is treated like a harmonic spring, with a specific spring constant $2K_b$ and equilibrium length $l_0$. Similarly, every angle is also harmonic. Dihedrals (also called torsions) are periodic and therefore modeled with a trigonometric expansion; each dihedral may have a specific $K_{\phi ,n}$ for one or more values of $n$. Finally, impropers (improper dihedrals) refer to groupings of four atoms meant to be coplanar, such as three atoms bound to a nitrogen atom in a trigonal-planar configuration.

Fig. 29 depicts representative bonded features with atoms labeled to permit presentation of formulas for computing bond lengths, valence angles, dihedrals, and impropers. The measures are all defined conventionally using vector arithmetic and the right-hand rule for vector cross-products. The bond vector $\rb _{ij}$ is defined

\begin{equation} \rb _{ij} = \rb _i - \rb _j \end{equation}

with magnitude

\begin{align} r_{ij} & = \left |\rb _{ij}\right | = \left |\rb _i-\rb _j\right |\\ & = \left [(r_{i,x}-r_{j,x})^2+(r_{i,y}-r_{j,y})^2+(r_{i,z}-r_{j,z})^2\right ]^\frac 12. \end{align}

For atoms $j$ and $k$ both bonded to $i$, the valence angle $\theta _{jik}$ is found via

\begin{equation} \cos \theta _{jik} = \frac {\rb _{ij}\cdot \rb _{ik}}{r_{ij}r_{jk}} \end{equation}

For atoms $j$ and $k$ bonded to each other, the dihedral angle defined by atom $i$ bonded to $j$ and atom $l$ bonded to $k$ is found by computing the angle of intersection between the two planes defined by the $i$-$j$/$j$-$k$ bonds and the $j$-$k$/$k$-$l$ bonds.

\begin{equation} \cos \phi _{ijkl} = \frac {(\rb _{ji}\times \rb _{kj})\cdot (\rb _{kj}\times \rb _{lk})}{\left |\rb _{ji}\times \rb _{kj}\right |\left |\rb _{kj}\times \rb _{lk}\right |} \end{equation}

Now, for improper dihedrals, this formula is usually sufficient because (a) $\phi _0$ is zero for coplanarity and thus $\phi $ should be close to zero, and (b) $\arccos ()$ always returns an angle on $[0,\pi ]$. However, for proper dihedrals (torsions), since these are fully periodic and not guaranteed to have symmetry of reflection about 0 or $\pi $, we must define $\phi $ on $[-\pi ,\pi ]$, and the $\arccos ()$ is not enough. So, the convention used to determine the sign of the torsion angle is to compute its sine:

\begin{equation} \sin \phi _{ijkl} = \frac {[(\rb _{ji}\times \rb _{kj})\cdot (\rb _{kj}\times \rb _{lk})]\cdot \rb _{kj}}{\left |\rb _{ji}\times \rb _{kj}\right |\left |\rb _{kj}\times \rb _{lk}\right |r_{kj}} \end{equation}

The sign of $\phi _{ijkl}$ is determined by its sine thusly: \begin{equation} \phi _{ijkl} = \begin {cases} -\cos ^{-1}\cos \phi _{ijkl}, & \text {if $\sin _{ijkl}<0$}\\ \cos ^{-1}\cos \phi _{ijkl}, & \text {otherwise} \end {cases} \end{equation}

Most Class-I potentials are based on the concept of context-specific “atom types”. For example, a carbon atom in an aliphatic chain is of a different type than a carbon atom in an aromatic ring, even though both are “carbon”. Then, bonds, angles, dihedrals, and impropers are classified by the types of their member atoms. This results in huge databases of potential parameters and the necessity to fully specify all bonded parameters for any system that is to be simulated. Luckily, most MD packages have helper routines to do just that.

Determining these parameters for a given force field is context-specific and therefore a bit of an art. Bonded potential parameters are mostly found by fitting the given analytical forms to series of single-point energy calculations peformed using high-level quantum-mechanical simulations. For example, parameters for the C1-C2-C3-C4 torsion in olefins was parameterized by sweeping the torsion angle of a quantum mechanical model of butane. Because most molecules involve a superpositition of such bonded potentials, their use confers a vibrational mode spectrum to molecules that can also be matched with experimental spectra for parameter adjustment.

Non-bonded parameters like charges and LJ parameters are adjusted so that errors in thermodynamics property calculations are minimized. This includes things like vacancy energies in solids, heats of vaporization in liquids, and volumetric equations of state. Most force-field developers strive to make their parameters sets as transferrable as possible, but one must always be aware of the context a force-field is most appropriate for. CHARMM and AMBER are used mainly to model biological molecules, while TRAPPE is pretty much strictly used for fluids. OPLS and GROMOS are used in a variety of settings.

It is important to note that all potentials considered so far assume each atom has a fixed partial charge (which may be zero). In fact, for class-I potentials, “charge” is an atom attribute wholly separate from its type, and they are usually derived from QM calculations on small fragments. Potentials for which the charge on each atom is a degree of freedom rather than a fixed quantity are described as “polarizable”. Polarizable potentials are not widely used, but there are good arguments for their requirement in situations where atoms with high polarizabilities might lead to local dipole moments that influence interatomic interactions.

One important polarizable model is the classical Drude oscillator, in which each atom is assigned a “ghost” particle tethered to its center that carries some of the charge. This allows each atom to behave like a literal “fluctuating dipole”, but it requires an extended Hamiltonian approach where the Drude particles’ equations of motion are solved along with those of the atom nuclei, making Drude-enabled simulations substantially more expensive than fixed-charge simulations. For example, the polarizable version of CHARMM is called CHARMM-Drude [26, 27].

8.2 Reactive Potentials

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.

8.2.1 Embedded Atom Method (EAM)

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

\begin{equation} \mathscr {U}(\rb ^N) = \sum _{i=1}^{N}\left [F_i(\rho _{h,i}) + \frac {1}{2}\sum _{j\ne i} \phi _{ij}(r_{ij})\right ] \end{equation}

where

\begin{eqnarray} F\left (\rho \right ) & = & -\sqrt {\rho }, \\ \rho _i & = & \sum _{j\ne i}g\left (r_{ij}\right ), \\ g\left (r\right ) & = & \exp {\left (-\beta r\right )}, \\ \label {eq:EAM:phi}\phi \left (r\right ) & = & V_2\left (r\right ) - 2F\left [g\left (r\right )\right ], \\ V_2\left (r\right ) & = & \left \{ \begin {array}{ll} V_{\rm ZBL}\left (r\right ), & \mbox {$r<r_1$}, \\ \alpha _0+\alpha _1r+\alpha _2r^2+\alpha _3r^3, & \mbox {$r_1\le r<r_2$}, \\ -2c\exp {\left (-{{\beta }\over {2}}r\right )}+\Phi _0\exp {\left (-\alpha r\right )}, & \mbox {$r\ge r_2$},\\ \end {array} \right .\\ \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{eqnarray}

The second term in Eq. ?? accounts for the fact that the electron density into which an atom is embedded does not include electrons from that atom itself. Eq. ?? 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]

8.2.2 Bond-Order Potentials

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]

8.2.3 ReaxFF

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:

\begin{equation} \mathscr {U}_{\rm ReaxFF} = U_{\rm bond} + U_{\rm over} + U_{\rm angle} + U_{\rm torsion} + U_{\rm VdW} + U_{\rm Coulomb} + U_{\rm specific} \end{equation}

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.

8.3 Case Study: Stillinger-Weber Silicon

As a precursor to the Tersoff formalism, perhaps the earliest attempt to use molecular simulation with a reactive potential to study a realistic atomic-scale model of silicon was due to Stillinger and Weber [38]. The C-code that implements the Stillinger-Weber (SW) potential for silicon is mdswsi.c. This code computes in reduced units as well; $\sigma $ = 0.20951 nm, $\epsilon $ = 2.1678 eV, and $m$ = 28.085 amu, which are appropriate for a system of pure silicon. One reduced unit of temperature, $\epsilon /k_B$, corresponds to 25156.74 K.

The main reason to introduce Stillinger-Weber silicon here is to give you a historical example of an implementation of a simple reactive force-field. Silicon forms 4-coordinated tetrahedral bonded structures. The SW potential includes three-body interactions that enforce this symmetry as well as permit breaking and reforming of bonds to produce defects, for example.

The total potential is expressed as two sums, one for unique pair interactions, and another for unique triplet interactions:

\begin{equation} \mathscr {U} = \sum _{i<j}v_2(r_{ij}) + \sum _{i<j<k}v_3({\bf r}_i, {\bf r}_j, {\bf r}_k) \end{equation}

The two-body term models the bonds: \begin{equation} v_2 \left (r\right ) = \left \{\begin {array}{ll} \epsilon A \left (Br^{-p}-r^{-q}\right )\exp \left [\left (r-a\right )^{-1}\right ], & \mbox {$r < a$}\\ 0 & \mbox {$r\ge a$} \end {array} \right . \end{equation} It is very much like a Lennard-Jones potential, only with different exponents and a “smooth cutoff” as the interatom separation distance, $r$, approaches some cutoff, $a$, given by the factor $\exp \left [\left (r-a\right )^{-1}\right ]$.

The three body models the angles, and is the sum of functions of each of the three angles of a triplet, $ijk$:

\begin{equation} v_3\left ({\bf r}_i, {\bf r}_j, {\bf r}_k\right ) = h_{jik} + h_{ijk} + h_{ikj} \end{equation}

Here I have employed the shorthand notation $\displaystyle h_{jik} \equiv h\left (r_{ij}, r_{ik}, \theta _{jik}\right )$. Note that, in the notation of this potential, $\theta _{jik}$ is subtended at ${\bf r}_i$, and $\cos \theta _{jik}^c = -\frac {1}{3}$:

\begin{equation} h\left (r_{ij},r_{ik},\theta _{jik}\right ) \equiv \begin {cases} \epsilon \lambda \exp \left [\frac {\gamma }{r_{ij}-a} + \frac {\gamma }{r_{ik}-a}\right ] \left (\cos \theta _{jik}-\cos \theta _{jik}^c\right )^2 & \mbox {if $r_{ij} < a$, and} \cr 0 & \mbox {otherwise} \end {cases} \end{equation}

One computes the angle-$j$ term, $h_{ijk}$, and the angle-$k$ term, $h_{ikj}$, by permuting the indices appropriately.

The parameters used in the original study by Stillinger and Weber are:

\begin{eqnarray} A & = & 7.049556277\\ B & = & 0.6022245584\\ p & = & 4\\ q & = & 0\\ a & = & 1.80\\ \lambda & = & 21.0\\ \gamma & = & 1.20 \end{eqnarray}

As in any MD simulation, one computes the force on any particle $i$ from the negative gradient of the total potential:

\begin{eqnarray} {\bf f}_i & = & -\nabla _{{\bf r}_i}\mathscr {U}\\ & = & -\sum _{j\ne i} \nabla _{{\bf r}_i}v_2\left (r_{ij}\right ) - \sum _{j\ne i}\sum _{k\ne i,j} \nabla _{{\bf r}_i}v_3\left ({\bf r}_i,{\bf r}_j,{\bf r}_k\right ) \end{eqnarray}

The two-body term for the $i$-$j$ interaction is only slightly more complicated than Lennard-Jones, due to the smooth cutoff. Here, assuming $i$ and $j$ are within interaction range ($r_{ij} < a$), we have for the force on $i$ due to $j$:

\begin{eqnarray} -\nabla _{{\bf r}_i} v_2\left (r_{ij}\right ) & = & -\frac {\partial }{\partial {\bf r}_i} \left \{\epsilon A \left (Br_{ij}^{-p}-r_{ij}^{-q}\right )\exp \left [\left (r_{ij}-a\right )^{-1}\right ]\right \}\\ & = & -\epsilon A \frac {{\bf r}_{ij}}{r_{ij}}\left (\left .\left [\frac {\partial }{\partial r} \left (Br^{-p}-r^{-q}\right )\right ]\right |_{r_{ij}}\exp \left [\left (r_{ij}-a\right )^{-1}\right ]\right .\\ & + & \left .\left (Br_{ij}^{-p}-r_{ij}^{-q}\right )\left .\left \{\frac {\partial }{\partial r}\exp \left [\left (r-a\right )^{-1}\right ]\right \}\right |_{r_{ij}}\right )\\ & = & -\epsilon A \frac {{\bf r}_{ij}}{r_{ij}}\left \{\left (-pBr_{ij}^{-p-1}+qr_{ij}^{-q-1}\right )\exp \left [\left (r_{ij}-a\right )^{-1}\right ]\right .\\ & - & \left .\left (Br_{ij}^{-p}-r_{ij}^{-q}\right )\exp \left [\left (r_{ij}-a\right )^{-1}\right ]\left (r_{ij}-a\right )^{-2}\right \}\\ \label {eq:swfij} & = & v_2\frac {{\bf r}_{ij}}{r_{ij}}\left [\frac {pBr_{ij}^{-p-1}-qr_{ij}^{-q-1}}{Br_{ij}^{-p}-r_{ij}^{-q}} + \left (r_{ij}-a\right )^{-2}\right ]\\ & \equiv & {\bf f}_{ij} \end{eqnarray}

Note that it is still true that ${\bf f}_{ij} = -{\bf f}_{ji}$. As an exercise, you should be able to show that the right-hand-side of Eq. ?? is well-behaved (that is, it does not diverge, and in fact vanishes) as $r_{ij}\rightarrow a^-$.

It is comparatively much more tedious to evaluate the three-body gradients:

\begin{eqnarray} -\nabla _{{\bf r}_i}v_3\left ({\bf r}_i,{\bf r}_j,{\bf r}_k\right ) & = & -\nabla _{{\bf r}_i}\left (h_{jik} + h_{ijk} + h_{ikj}\right )\\ \label {eq:sw3_partials} & = & -\left (\frac {\partial h_{jik}}{\partial {\bf r}_i} +\frac {\partial h_{ijk}}{\partial {\bf r}_i} +\frac {\partial h_{ikj}}{\partial {\bf r}_i}\right )\\ & \equiv & {\bf f}_{i\leftarrow jk} \end{eqnarray}

The triplet forces on the other members of the triplet, ${\bf f}_{j\leftarrow ik}$ and ${\bf f}_{k\leftarrow ij}$, are defined analogously.

Each of the partials in Eq. ?? is unique: \begin{equation} \label {eq:sw3jik} \begin {array}{lll} \mbox {\raisebox {-0.5cm}{\includegraphics [width=1cm]{figpng/hjik.png}}\ \ } \displaystyle {\frac {\partial h_{jik}}{\partial {\bf r}_i}} & \displaystyle = &\displaystyle - \gamma h_{jik}\left [ {{{\bf r}_{ij}}\over {r_{ij}}}{{1}\over {\left (r_{ij}-a\right )^2}} + \frac {{\bf r}_{ik}}{r_{ik}}\frac {1}{\left (r_{ik}-a\right )^2} \right ] \\ & \displaystyle + & \displaystyle 2\lambda \exp \left ( \frac {\gamma }{r_{ij}-a} + \frac {\gamma }{r_{ik}-a} \right ) \left (\cos \theta _{jik}-\cos \theta _{jik}^c\right ) \\ & \displaystyle \times & \displaystyle \left [ \frac {{\bf r}_{ij}}{r_{ij}}\frac {1}{r_{ik}} +\frac {{\bf r}_{ik}}{r_{ik}}\frac {1}{r_{ij}} -\left (\frac {{\bf r}_{ij}}{r_{ij}}\frac {1}{r_{ik}} +\frac {{\bf r}_{ik}}{r_{ik}}\frac {1}{r_{ij}}\right ) \cos \theta _{jik}\right ] \end {array} \end{equation} \begin{equation} \label {eq:sw3ijk} \begin {array}{lll} \mbox {\raisebox {-0.5cm}{\includegraphics [width=1cm]{figpng/hijk.png}}\ \ } \displaystyle \frac {\partial h_{ijk}}{\partial {\bf r}_i} & \displaystyle = & \displaystyle - \gamma _j h_{ijk}\left [ \frac {{\bf r}_{ij}}{r_{ij}}\frac {1}{\left (r_{ij}-a\right )^2} \right ] \\ & \displaystyle + & \displaystyle 2\lambda _j\exp \left ( \frac {\gamma _j}{r_{ij}-a} + \frac {\gamma _j}{r_{jk}-a} \right ) \left (\cos \theta _{ijk}-\cos \theta _{ijk}^c\right ) \\ & \displaystyle \times & \displaystyle \left [ \frac {{\bf r}_{jk}}{r_{jk}}\frac {1}{r_{ij}} +\frac {{\bf r}_{ij}}{r_{ij}}\frac {1}{r_{ij}} \cos \theta _{ijk}\right ] \end {array} \end{equation} \begin{equation} \label {eq:sw3ikj} \begin {array}{lll} \mbox {\raisebox {-0.5cm}{\includegraphics [width=1cm]{figpng/hikj.png}}\ \ } \displaystyle \frac {\partial h_{ikj}}{\partial {\bf r}_i} & \displaystyle = & - \displaystyle \gamma _k h_{ikj}\left [ \frac {{\bf r}_{ik}}{r_{ik}}\frac {1}{\left (r_{ik}-a\right )^2} \right ] \\ & \displaystyle + & \displaystyle 2\lambda _k\exp \left ( \frac {\gamma _k}{r_{ik}-a} + \frac {\gamma _k}{r_{jk}-a} \right ) \left (\cos \theta _{ikj}-\cos \theta _{ikj}^c\right ) \\ & \displaystyle \times & \displaystyle \left [ -\frac {{\bf r}_{jk}}{r_{jk}}\frac {1}{r_{ik}} +\frac {{\bf r}_{ik}}{r_{ik}}\frac {1}{r_{ik}} \cos \theta _{ikj}\right ] \end {array} \end{equation} Note that the right hand sides of each of Eqns. 204, 205, and 206 vanish when the appropriate $r$’s are greater than $a$’s, as in Eq. 212. As an exercise, it is easy to show that ${\bf f}_{i\leftarrow jk} + {\bf f}_{j\leftarrow ik} + {\bf f}_{k\leftarrow ij} = 0$.

Let us now run two MD simulations for 10,000 time steps; one with and and the other without three body forces, at reduced initial temperature $T$ = 0.12 and reduced density $\rho $ = 0.46 for a small system of 216 atoms. The code mdswsi.c can initialize atoms on a diamond-cubic lattice; the DC unit cell has 8 atoms, so the number of atoms one should specify for a perfect crystal is 8 times any product of three integers. Using the -nc #,#,# switch we indicate how may DC unit cells we want in the $x$, $y$, and $z$, directions, respectively. The code also has a switch -do-three-body which is by default on, but can be turned off by passing a zero. Running this twice generates two logs and two trajectories:

$ ./mdswsi -nc 3,3,3 -ns 10000 -fs 10 -prog 10 \
  -do-three-body 1 -traj yes3.xyz > yes3.log
$ ./mdswsi -nc 3,3,3 -ns 10000 -fs 10 -prog 10 \
  -do-three-body 0 -traj no3.xyz > no3.log

Fig. 30 shows two 3-D system representations. The first is a perfect DC lattice with 216 atoms, where a dot is drawn every 10 time steps for each atom, and each atom is colored uniquely. The second shows the same view of a system for which the three-body forces are turned off; notice that it appears liquid-like. Three-body forces are required to stabilize the DC lattice since each atom only has four nearest neighbors.

Figure 30: Two views of a 10,000-time-step NVE MD simulation of a system of 216 silicon atoms initialized on a diamond-cubic lattice; (left) with three-body force and (right) without three-body forces. Each point is an atom position, and atom positions are rendered every 10 time-steps, and each atom is colored uniquely.

Fig. 31 shows a plot of two-body and three-body potential energy vs. time from the first simulation, and a plot of instantaneous temperature vs. time from both simulations on the right. Notice that the three-body forces act to keep the system oscillating in its crystalline state, and the lack of three-body forces results in system melting. This latter occurrence is because the DC lattice is a fairly unfavorable configuration for a system with only two-body forces active.

Figure 31: Results of a 216-atom NVE MD simulation of SW silicon initialized on a DC lattice. (Left) Two-body and three-body contributions to the potential energy vs time. (Right) Temperature vs. time for systems with and without three-body forces active.

The code mdswsi.c is very slow when three-body forces are turned on because it uses a brute-force $N^3$ loop. In practice, any code that implements two- and three-body potentials with short-range cutoffs will use data structures like the neighbor list and the link-cell algorithm to make the pair and triplet loops more efficient. And although SW silicon is historically important, there are much better (and more complicated) potentials out there for silicon-based systems, such as COMB [34].

← Previous Contents Next →