Intent of this issue

The introduction of physically sound interatomic potential energy functions that go beyond simple pair-additive interactions (e.g., many-body potentials) beginning in the 1980s opened tremendous modeling capabilities that continue to shape new directions and create critical breakthroughs in materials research.Reference Vitek and Srolovitz^1–Reference Brenner, Shenderova, Areshkin, Lipkowitz and Boyd⁵ The key to their proven usefulness is a combination of relative accuracy in reproducing important structures (including defects) across a wide range of material types and their overall computational efficiency. This combination of features allows atomic simulations that are large enough to explore phenomena such as correlated dynamics associated with plastic flow in metals and accurate enough to be compared to specific materials and structures. With this capability, continuum concepts can be tested at the atomic scale, experimental results interpreted in new ways, and virtual experiments carried out that are at the forefront of the development of new materials.

Our intentions with this issue of MRS Bulletin are to celebrate the rapid succession of many-body potentials that were introduced in the early to mid-1980s, to review the impact that these potentials have had on research carried out by the materials community in general, and to outline where the field is headed in the next three decades. Contributions in this issue are included from some of the original potential developers, from researchers who have made seminal contributions to materials research using these potentials, and from researchers who are at the forefront of developing and applying the next generation of methods. We expect that this “continuum” of modeling concepts and applications will help guide and inspire the next generation of computational materials scientists and engineers as they expand this capability to new and exciting areas of materials research.

Brief history of atomic simulation

To understand how this field has progressed to its current state and to predict where it is going, it is useful to briefly review the history of atomic simulation using classical trajectories (see the sidebar). The first reported study using classical trajectories to model a chemical process was published in 1936 by Hirschfelder, Eyring, and Topley.Reference Hirschfelder, Eyring and Topley⁶ At this point in the development of atomic simulation, the concept of a continuous potential energy hyper-surface over which atoms moved was still relatively new, and they were interested in understanding how to generate useful potential energy surfaces that define atomic motion for chemical reactions and how to best estimate rates from these surfaces. As part of their studies, they used a single, hand-calculated classical trajectory to study the reaction H + H₂ → H₂ + H constrained to a linear configuration. This reaction was chosen because with only three electrons, it is relatively simple. The linear configuration was used because it both reduces the number of degrees of freedom that had to be followed and corresponds to the lowest energy barrier for reaction (and hence the configuration that most contributes to the overall reaction rate). Their potential energy surface was generated from an analytic approximation to a quantum chemical calculation in which the exchange and Coulomb integrals are approximated by pair functions that depend on distance. The potential energy surface that they used included a potential energy “well” (known later as Lake EyringReference Truhlar and Wyatt⁷) that yielded a stable H₃ molecule (Figure 1). One of the consequences of this well is that the reaction never completed the hydrogen atom exchange during the course of the trajectory that was published; instead the system remained an H₃ intermediate. Nonetheless, this single trajectory was a beautiful illustration of the utility of classical trajectories for understanding complex chemical dynamics.

Figure 1. Three-dimensional representation of the potential energy surface used by Eyring and co-workers to model H + H₂ → H₂ + H.⁶ Each horizontal axis corresponds to the distance between the outer hydrogen atom and the center atom, and the vertical axis corresponds to potential energy. The different shades represent potential energy for the interacting atoms. The darker the color, the lower the potential energy. This is a representation rather than the exact surface used in Reference Reference Hirschfelder, Eyring and Topley6, so units are not given on the axes. “Lake Eyring” is a dip in the potential energy that corresponds to a stable (and non-physical) H₃ molecule.

Since the early calculations of Eyring and co-workers, classical trajectories have proven to be a powerful tool in chemical physics not only for calculating reaction rates, but also for testing and refining theories of chemical kinetics. For example, the central concept of transition state theory (TST) is that the kinetics for a given reaction can be calculated from a single configuration along a reaction path. This configuration is often taken as that at the highest intermediate energy along the reaction path (e.g., the activated complex), but this is often not the best choice due to recrossings back toward reactants. Classical trajectories have been used for a wide range of systems to test the assumptions of TST and to refine the theory by finding more suitable configurations for the transition state.Reference Garrett and Truhlar^8, Reference Truhlar and Wyatt⁹

During the 1940s and 1950s, classical trajectories using pair interactions were utilized to evaluate some of the basic principles of statistical mechanics, including the rate at which equilibrium conditions are established (if at all) for relatively simple systems.Reference Berman and Izrailev^10, Reference Alder and Wainwright¹¹ From a macroscopic viewpoint, the properties of many-body systems are irreversible, as dictated by the second law of thermodynamics. Classical trajectories, however, are in principle deterministic and reversible, implying that routes to systems involving an increase in order (decreased entropy) are possible. Among the conclusions of these studiesReference Berman and Izrailev^10, Reference Alder and Wainwright¹¹ is that some relatively small simplified systems, such as short harmonic chains with weak anharmonic terms, may never fully equilibrate, and therefore traditional statistical mechanics is not generally valid for these systems. On the other hand, numerical approximations for solving classical equations of motion for particles can ensure an approach to equilibrium due to effects such as small errors arising from the truncation of numbers during the calculation, similar to small disturbances (e.g., external temperature fluctuations) to physical systems that help drive a system toward equilibrium.

While making important contributions to chemistry and physics, classical atomic trajectories did not, in general, have the same early impact on the field of materials science and engineering due mainly to limited computational power relative to what is needed to model the structures of traditional interest to materials scientists. One exception to this was simulations by Vineyard and co-workers at Brookhaven National Laboratory, first published in 1960, who used classical dynamics to examine irradiation damage in metals (Figure 2).Reference Gibson, Goland, Milgram and Vineyard¹² These calculations were performed in three dimensions with continuous pair additive potentials (as opposed, for example, to hard spheres). Soon after the Brookhaven calculations, Rahman at Argonne National Laboratory used classical trajectories with simple pair additive potentials to characterize the structure and atom motion in a model liquid.Reference Rahman¹³ Like the work of Eyring and co-workers on chemical dynamics, Rahman’s simulations demonstrated the utility of classical simulations for understanding dynamics in fluids and for generating quasi-experimental data that could be used to test and refine existing theories of liquid structure and correlated dynamics. As computing power continued to increase, Rahman and others extended this initial work to more complicated fluids, including water,Reference Stillinger and Rahman¹⁴ as well as biological systems.Reference McCammon, Gelin and Karplus¹⁵

Figure 2. Illustration of atomic trajectories from a 500 atom simulation of radiation damage in copper carried out by Vineyard (based on figures in Reference Reference Gibson, Goland, Milgram and Vineyard12). Spheres represent initial atomic positions in a crystal; the lines and dots trace atomic motion during a collision cascade started by motion of the atom represented by the red sphere.

There are several common themes that emerge from these early atomic simulations. First, molecular simulation was emerging as a powerful tool in physics and chemistry that complemented both theory and experiment. Not only could it provide data but it could be used to test and refine theories under exceedingly well-controlled conditions. Second, the development of atomic simulation as a tool was driven in large part by the availability of increasingly powerful computing resources, such as those at US national laboratories. Finally, it became clear that demands placed on a potential energy function depended very strongly on the application. For chemistry, well-refined potential energy surfaces are needed to provide qualitative data, but even approximate surfaces can yield useful data for exploring general phenomena. Similarly, in statistical mechanics, simple potentials can yield exceedingly useful information.

Over the same time period that simulations were emerging as important tools in chemistry and physics, new and sophisticated dislocation-based theories of materials deformation were being developed. Why did atomic simulation not have the same impact on materials science over the same time period, in particular in understanding material deformation? There are two likely leading reasons. First, stress fields, which play a critical role in dislocation dynamics and plasticity, require relatively large systems to be completely contained, and computing resources were not sufficient to allow simulation of necessary system sizes. Second, even if sufficient computing resources were available to model relatively large systems, the pair potentials that were being used have well-established limitations for quantitatively reproducing the properties of crystalline solids that are required to satisfactorily describe dislocations. These limitations include elastic constants that are forced to obey the Cauchy relations (C ₁₂ = C ₄₄, where C ₁₂ and C ₄₄ are elastic constants) in cubic systems, generally inaccurate surface and vacancy formation energies and relaxations, and a general inability to describe as the lowest energy state all but close-packed lattices.

There was a convergence of capabilities in the 1980s that led to a rapid rise in the use of atomic simulation in materials science. First, there were new theories of interatomic bonding,Reference Norskov and Lang¹⁶ out of which came many-body potential energy expressions that did not have the inherent limitations of pair potentials for quantitatively describing materials properties but had computing demands that were slightly higher than pair potentials.Reference Carlsson, Ehrenreich and Turnbill¹⁷ Second, major computing resources that had driven much of the early work at the US national laboratories (and at places such as IBM) were becoming readily available to academic and other institutions. Third, there was a pressing need for new materials with tailored properties down to the atomic scale that could not be adequately calculated with the existing theory and modeling tools. Finally, researchers such as the developers of the embedded-atom method (EAM)Reference Daw, Foiles and Baskes¹⁸ not only published papers, they also freely released their computer codes to researchers (Figure 3). Further details about the EAM are given in the contribution by Foiles and Baskes in this issue.

Figure 3. Atomic configurations from a molecular dynamics simulation of dislocations created by high strain rate compression cutting an AlCu precipitate in an Al matrix. The precipitate is about 5 nm thick. Only atoms with local symmetry different from the bulk are shown. The compression comes from the left, which drives the dislocations from left to right. The simulations were carried out using the embedded-atom method potentials and the large-scale atomic/molecular massively parallel simulator (LAMMPS) code.

An interesting aspect of this research convergence in the 1980s was that the potential functions that enabled simulations targeted at particular materials came almost simultaneously from multiple groups in different countries. The Finnis-Sinclair potentials from the United Kingdom,Reference Finnis and Sinclair¹⁹ the EAM potentials from the United States,Reference Daw and Baskes²⁰ the glue potentials from Italy,Reference Ercolessi, Tosatti and Parrinello²¹ the effective medium potentials from Denmark,Reference Jacobsen, Norskov and Puska²² and the corrected effective medium theory potentials from the United StatesReference Kress and DePristo²³ all looked very similar in form, despite different theoretical derivations and intended applications. Even the bond-order potentials developed around that same time,Reference Tersoff^24–Reference Khor and Das Sarma²⁷ despite initially looking very different in form from functions such as the EAM,Reference Brenner²⁸ came from similar moments expansion of the local electronic density of states.Reference Brenner, Shenderova, Areshkin, Lipkowitz and Boyd⁵ One exception to this common form (see the sidebar) is the Stillinger-Weber potential that was introduced to model silicon.Reference Stillinger and Weber²⁹

Today’s challenges and tomorrow’s opportunities

The topics outlined in this issue were chosen to illustrate the research convergence of the 1980s in this area: the theoretical formalisms that led to new many-body potentials, the transition of these theories to potential functions that are useful for large-scale simulations, the codes that spread these potentials throughout the research community, and the materials and applications of atomic modeling that have been opened by these capabilities. In the first article, Finnis discusses some of the quantum mechanical concepts that have led to the development of effective many-body potentials as well as some of the functional forms and emerging methods being used to fit parameters for these functions. This is followed by a contribution from Foiles and Baskes that focuses on the widely used EAM and related potentials, including some contemporary applications in materials science. Next is a contribution from Pastewka et al., in which the bond-order potentials and their application to covalently bound ceramic materials are discussed in detail. This is followed by a contribution from Shin and co-workers that discusses two relatively new dynamic charge-transfer methods: the reactive force field (ReaxFF) and charge optimized many-body (COMB) potentials. The focus is on atomic-scale modeling of systems with non-uniform bonding environments, including surface chemistry between covalently bound molecules and metallic or metal-oxide surfaces. In the final contribution, Plimpton and Thompson examine the performance and application of multiple potentials, including those discussed in this issue of MRS Bulletin, that are incorporated within the highly parallelizable, open-source large-scale atomic/molecular massively parallel simulator (LAMMPS) molecular dynamics software.

These articles illustrate the capabilities of state-of-the-art many-body potentials and their application to fundamental problems in materials science. As was the case in the 1980s, the trifecta of requirements of many-body reactive potentials by the computational community remains (1) incorporation of appropriately high levels of theory, (2) implementation in a manner that enhances practical application, and (3) utilization in software that takes advantage of evolving computing capabilities. As discussed in the article by Plimpton and Thompson, the computational cost of many-body potentials developed over the last three decades varies by several orders of magnitude. This provides a wide array of new opportunities for the simulation community to take advantage of the associated advances in computing power to study ever larger materials systems (∼10¹¹ atoms) with well-established and less expensive potentials, or to examine smaller systems (∼10⁴ atoms) with recent, more computationally intensive potentials that are intended for use on more complex systems.

New directions and capabilities have recently emerged that are expected to greatly influence atomistic materials science and engineering over the next three decades. One example of a new direction is the combination of atomic simulations with traditional engineering approaches, such as finite element analysis that is beginning to blur the distinction between atomic and continuum concepts.Reference Miller and Tadmor³⁰ Among the emerging capabilities is the continuing ad hoc incorporation of electronic effects into analytic potentials that is going significantly beyond the moments expansions of the 1980s that is discussed in the following contributions. Examples include charge equilibration approaches and electronic force fields that can simulate excited electronic states in complex, heterogeneous systems.Reference Su and Goddard³¹ Given the importance of excited states on the dynamic properties of complex material systems, such as photovoltaic devices, it is expected that the latter capabilities will find increasing utilization by the computational materials science community.

Enabling all of these advances are continuing enhancements in computing capabilities; these capabilities include faster and more numerous parallel processors, larger data storage with faster access speeds, specialized visualization and data analysis hardware, and the emerging availability of these resources through cloud computing. Graphics processing units (the processors traditionally used in video cards), for example, are currently providing individual researchers access to a high degree of parallel processing at a cost per processor that is substantially less than traditional computer processing units. At the same time, high-end computing is looking to transition from the peta- to the exascale (i.e., from 10¹⁵ to 10¹⁸ floating point operations per second) with access to these resources through the computing cloud. The next three decades will be an exciting time for atomic simulations in materials science and engineering.