# Mapping the stochastic response of nanostructures

See allHide authors and affiliations

Edited by L. B. Freund, University of Illinois at Urbana–Champaign, Urbana, IL, and approved March 7, 2014 (received for review February 1, 2014)

## Significance

A challenge in nanotechnology is that nanoscale structures exhibit a rate-dependent and apparently random (stochastic) load-deflection behavior when repeated experimental measurements are compared. The variability in mechanical response is due to the extreme nonconvexity of the nanostructure’s energy versus structure function. The many minima of this function correspond to equilibrium states, accessible during loading, that create many possible pathways. Traditional atomistic simulation techniques cannot systematically address this complexity. Here, we explore a new method based on constructing an “equilibrium map” (akin to a phase diagram) that efficiently and systematically explores a nanostructure's stochastic, rate-dependent response to specified loading conditions. We demonstrate the method's capabilities and its surprisingly complex results for the case of a nanoslab of nickel under uniaxial compression.

## Abstract

Nanostructures are technological devices constructed on a nanometer length scale more than a thousand times thinner than a human hair. Due to the unique properties of matter at this scale, such devices offer great potential for creating novel materials and behaviors that can be leveraged to benefit mankind. This paper addresses a particular challenge involved in the design of nanostructures—their stochastic or apparently random response to external loading. This is because fundamentally the function that relates the energy of a nanostructure to the arrangement of its atoms is extremely nonconvex, with each minimum corresponding to a possible equilibrium state that may be visited as the system responds to loading. Traditional atomistic simulation techniques are not capable of systematically addressing this complexity. Instead, we construct an equilibrium map (EM) for the nanostructure, analogous to a phase diagram for bulk materials, which fully characterizes its response. Using the EM, definitive predictions can be made in limiting cases and the spectrum of responses at any desired loading rate can be obtained. The latter is important because standard atomistic methods are fundamentally limited, by computational feasibility, to simulations of loading rates that are many orders of magnitude faster than reality. In contrast, the EM-based approach makes possible the direct simulation of nanostructure experiments. We demonstrate the method’s capabilities and its surprisingly complex results for the case of a nanoslab of nickel under compression.

The current revolution in nanotechnology promises to bring great benefits to mankind. Much of the excitement focuses on the creation of nanostructures and nanodevices that, due to their small size, have unique properties relative to macroscopic bulk materials. However, a challenge in harnessing such properties is that nanostructures can behave in a stochastic or apparently random manner when loaded. Thermal and mechanical fluctuations can cause nanostructures to morph from one deformed structure to another and nominally identical loadings can lead to many different response sequences. An example of this is the uniaxial compression of nanopillars that exhibit stochastic behavior when internal defects are nucleated (1). This indeterminacy in response stems from the fact that for a given state of external loading there are a great many ways for the atoms comprising the nanostructure to arrange internally. Mathematically, this can be understood by considering the dependence of the system’s potential energy on the positions of the atoms ** u** and a scalar

*λ*representing the applied loading, i.e., . The position of every atom has three coordinates, therefore

**has components (or “degrees of freedom”) less those fixed by imposed constraints. For a given value of**

*u**λ*, one can imagine plotting this function in an -dimensional space (this is of course not practically possible). [At most one can visualize this for , where takes the form of a topographical map at some fixed value of λ.] One then obtains an

*n*-dimensional hypersurface called the “potential energy surface” (PES), which is a conceptually useful tool to understand the behavior of physical systems. In particular, minima on this surface [i.e., values of

**for which is a local minimum] correspond to metastable equilibrium states of the nanostructure. Physically these minima correspond to all of the internal defect structures, crystal phases, and the like, that can occur due to the applied loading. Saddle points on the PES correspond to (unstable) transition states between equilibrium states. The stochastic response that is observed experimentally and computationally is due to the system hopping from one minimum to another through transition states.**

*u*Atomistic and multiscale methods aim to predict a system’s response to a changing load directly from the atomic interactions represented by the PES. In static simulations, the loading parameter is incremented in stepwise fashion, , and the equilibrium configuration in the new state is obtained by minimizing using as an initial guess for ** u** the equilibrium configuration obtained for . (Because the PES contains multiple minima, i.e., it is nonconvex, the choice of initial guess is critical.) Alternatively, a molecular dynamics (MD) simulation can be performed in which a time-dependent trajectory of the system is obtained by numerically integrating the equations of motion, . Here the forces acting on the degrees of freedom are the negative of the PES gradient, . The end result for both static and dynamic simulations is the response trajectory, , i.e., the arrangement of the atoms as a function of the applied load.

We have mentioned that the PES is nonconvex. In fact, it is extremely nonconvex. This poses practical and conceptual problems to the types of simulations described above which are often overlooked. Stillinger and Weber (2, 3) estimate that a system of *N* atoms of a single species exhibits distinct structural minima where *ν* is a positive number. Based on this estimate (4), one gram of argon would have of the order of distinct structural minima! Such is the complexity of the problem. Unfortunately, the simulation procedures described above provide only one of the myriad of possible equilibrium pathways. Even more troubling is the fact that the obtained pathway will depend on the details of the numerical solver used. For dynamic trajectories there is the added problem that simulated loading rates are typically many orders of magnitude larger than the experimental rate due to the short time scale of atomic vibrations that limit the numerical integration. This can lead to qualitatively different behavior in cases where the response is rate dependent.

To address these limitations and explore the implications of PES nonconvexity, we propose an approach that is completely different from current-day static minimization techniques and thermodynamic sampling methods (such as MD and Monte Carlo). Rather than statically or dynamically generating a single response trajectory, we use branch-following and bifurcation (BFB) algorithms to generate a diagram we call the equilibrium map (EM), which contains information on all possible response trajectories. Specifically, by definition, the EM includes representatives of all stable and unstable equilibrium states as a function of the applied loading. [The term “bifurcation diagram” is commonly used for such a graph. We introduce the term “EM” to reflect the qualitative difference in the size of the maps we generate and their application.] Similar to the way in which a phase diagram shows the available bulk states of matter as a function of pressure and temperature, the EM describes the possible metastable states associated with a given value of the applied loading. For example, for nanopillar compression, the EM can be represented as a plot of stress versus strain where the metastable states associated with different internal defect arrangements (microstructures) appear as disconnected segments. Here *λ* is either the stress or the strain, depending on whether the load or displacement is controlled. (An example of an EM for a similar problem is given later in Fig. 3.)

Given a system’s EM, we can approach the question of its response to loading in a more systematic fashion. Rather than limiting the analysis to a single stochastic (static or dynamic) response trajectory that may be an artifact of the solution algorithm, we classify the response into three regimes based on the loading rate: (*i*) quasistatic process (QP); (*ii*) quenched dynamics (QD); and (*iii*) driven dynamics (DD). As the name suggests, QP corresponds to the thermodynamic limit wherein for an infinitesimal increase in *λ*, the system is given infinite time to sample all possible configurations and to settle into the global minimum of the PES. Physically this corresponds to a system loaded infinitely slowly and in contact with a heat bath (at ). In the second regime, QD corresponds to the case where kinetic energy is instantaneously removed from the system. At a point of instability, the system can only roll down the PES to a local minimum without the ability to overcome energy barriers. This is an idealization of a system with very high damping and very high thermal conductivity such that any mechanical vibrations are instantly damped out and the heat released is transferred to the surroundings at a rate that is very high compared with the rate of loading. QP and QD correspond to two well-defined physical paths through the set of available stable segments of the EM. (This is in contrast with static energy minimization schemes that would essentially choose a random path through the EM.) In addition to the idealized cases described above, the last case of DD refers to the real-world experimental case where a system is loaded at a finite rate. This is explained further in *Results and Discussion*.

To illustrate this novel EM-based simulation strategy, a model problem of a nanoslab subjected to uniaxial compressive displacement is presented. The set of possible equilibrium states found is much more complex than first expected and is a vivid illustration of the rich behavior of which nanoscale systems are capable. We discuss the results within the context of the loading regimes identified above.

## BFB Investigation: Theoretical Background

Starting from a physically observable (generally an unloaded and thus a stable) structure, equilibria are obtained by solving the force equilibrium equation,The implicit function theorem (5) states that, under certain continuity conditions, solutions to Eq. **1** form locally unique continuous curves, . These curves are called “equilibrium paths” or “branches,” where *α* is analogous to an arc length or distance along the curve. Solution points are locally stable if all of the associated eigenvalues of the Hessian (second derivative matrix), , are positive, and unstable otherwise. Stable points correspond to local minima on the PES (at a given value of loading *λ*), whereas unstable points correspond to local maxima and saddle points. For more details see, for example, refs. 6⇓⇓⇓⇓⇓⇓⇓⇓–15. Points on the path where the Hessian becomes singular are termed critical points. These points are often associated with a transition between stable and unstable equilibrium. If two or more distinct paths meet at a critical point it is termed a bifurcation point. A turning point is a point on the path where the tangent is perpendicular to the *λ* axis.

The set of all stable and unstable equilibrium paths constitutes the EM. The EM is thus a complex interconnected network of solution curve segments that join together at bifurcation points. The solution curve passing through the initial configuration is termed the principal equilibrium path and the bifurcation points along this path are called principal bifurcation points. The paths that branch out of these points are first generation paths and those that branch out of those are second generation paths, and so on.

As a simple example to help clarify these ideas, consider the schematic one-dimensional model shown in Fig. 1. The surface represents the dependence of an energy function on a single degree of freedom *u* and loading parameter *λ*. A surface of this type is typical of systems undergoing buckling, such as the classic case of a thin-walled axially compressed cylindrical shell. Using the physics of such shells as a motivating example, we understand the figure as follows. The shell exists in a stable cylindrical configuration at low values of axial compression (low values of *λ*) and in a buckled configuration at high compression. The curves in Fig. 1 (*Insets*) show the PES at different compressive loads. The transition from the single-well cylindrical configuration to two symmetry-related variants of the buckled configuration is clear. The diagram at the bottom of the figure shows the EM for the system. The solid lines correspond to stable (local energy minimizing) equilibrium configurations and the dashed lines to unstable (local maximum or saddle point) equilibrium configurations. As the compressive load is continuously increased, the well associated with the cylindrical configuration becomes more shallow and eventually transforms into a maximum losing its stability at high compression. At the point of change in topography from a minimum to a maximum, the Hessian becomes singular at a bifurcation point. The two dashed golden-brown segments bifurcating from the primary path represent the first generation paths. Points along this path correspond to the maxima separating the cylindrical and buckled configurations and are thus unstable. These paths connect to the stable buckled configuration minima at turning points.

Given the extreme nonconvexity of the PES, the construction of the EM for realistic systems is a daunting task that is only now becoming possible due to recent developments in BFB algorithms for massively parallel supercomputers. In practice, an EM is generated by tracing a single solution curve and identifying all bifurcation points, then returning to the bifurcation points, tracing out the connecting curves and their bifurcation points, and so on in a recursive fashion. This process is ideal for parallelization because the calculation of each curve is independent of the others.

There is a rich literature dealing with the characterization of local bifurcations. For example, see refs. 6⇓⇓⇓⇓⇓⇓⇓⇓–15. There is also a mature literature on numerical branch-following (continuation) algorithms, including refs. 5, 16⇓⇓–19. These methodologies have been applied in a variety of settings including problems in crystal and atomic-scale stability (20⇓⇓⇓⇓–25).

The exploitation of symmetry in the construction of the EM is often a critical step. In fact, the existence of bifurcation points is usually due to the presence of symmetry (as exemplified in Fig. 1). Once the symmetry group of a PES function is identified, the theory of equivariant systems (17, 26) can be used to reduce the dimension of the equilibrium equations and significantly improve the computational efficiency of the branch-following procedure. This reduction of dimension also regularizes the equilibrium equations near bifurcation points allowing for the accurate numerical determination of these points. Further, the use of equivariant bifurcation theory (14, 17, 26⇓–28) shows that once one equilibrium path has been computed all symmetry-related paths may be constructed without further computation. These theoretical and computational techniques must play a central role in any automated approach for the generation of the EM.

## Compression of a Face-Centered Cubic Nanoslab of Nickel

As an example of the EM-based simulation methodology, we focus on a simple atomistic application as a model problem. A face-centered cubic (fcc) nickel nanoslab of finite width and height and infinite thickness in the out-of-plane direction is uniformly compressed by a controlled displacement of its top and bottom surfaces as shown in Fig. 2. [Note that the EM can be used for the interpretation of both displacement and load control problems. However, in this publication we restrict our discussion to the displacement control case.] The nanoslab contains 104 unconstrained atoms (per unit depth) for a total of 312 degrees of freedom. For this problem, *λ* is the compressive displacement of the nanoslab and ** u** is the column vector of nuclear coordinates. The atomic interactions are modeled by an embedded atom method potential for nickel due to Angelo et al. (29) and Baskes et al. (30) modified to have smooth first and second derivatives.

## Results and Discussion

Fig. 3 shows the EM of the nickel nanoslab as a plot of the nominal compressive stress (compressive force *F* per unit reference area *A*_{0}) versus the applied compressive strain. Stable segments are shown in blue and unstable segments in golden brown. The loading starts at the point “S” where the applied displacement *λ* is zero and proceeds up and to the right with increased compression. [Note that at point S the nominal compressive stress is negative (tensile). This is due to the constraint of the applied displacements which is preventing surface relaxation.] Fig. 3 (*Inset*) shows in more detail the point of initial loss of stability on the primary path. This EM contains 8,026 solution curves (after removing symmetry-related curves), of which 1,306 are stable with a total of 5.8 million distinct solution points. It was generated on 256 processors of a parallel cluster over about 2 wk of wall time and to our knowledge constitutes the largest BFB simulation performed to date. With the data from the EM, we present selection rules of the stable paths representing the different loading regimes discussed earlier.

We first set out to identify the QP path where at each value of *λ* the stable segment in the EM with the lowest energy is selected. This corresponds to the lower envelope on the energy–strain plot in Fig. 4 where only stable segments of the EM are shown. To differentiate adjacent equilibrium paths, alternate segments are colored red and green and are joined by black lines if discontinuous in this projection. Fig. 4 (*Insets*) show some of the details of the lower envelope structure. The corresponding path represented in (compressive) stress–strain space is shown in Fig. 5. Note that the transition between stable segments occurs before each individual segment loses stability at its end point. The transition is dictated by the changing ground state (lowest-energy) structure as the load is increased.

Each stable segment traversed by the QP path in Fig. 5 corresponds to a globally stable configuration of the nanoslab. The configurations at the start and end of each QP segment are shown in Fig. 6. The black segments cannot be associated with any specific reaction mechanisms, but rather are complex global minimization processes connecting the initial and final structures in the presence of thermal vibrations over long time scales. The color coding in Fig. 6 is based on a local crystal structure indicator (31) (see the legend at the bottom of the figure). The nanoslab starts out in the fcc phase (cyan) and transitions through a series of defective states involving the passage of partial dislocations that leave behind stacking faults (locally hexagonal close-packed structures shown in red). The nanoslab reverts back to a barrel-shaped fcc structure at higher loading. The barreling effect is due to the fixed boundary conditions imposed at the top and bottom ends of the nanoslab.

The next case is the QD path describing an overdamped system. Starting from a stable initial configuration in the EM (point S in Fig. 3), we increase *λ* and follow the stable segment until a loss of stability is encountered. At this critical value, , the system is perturbed in the direction corresponding to the negative curvature of the PES. [There are actually two directions corresponding to perturbations in the positive and negative directions of the unstable eigenvector. However, for our system, these are associated with symmetry-related states and therefore only one direction needs to be followed.] Starting from the perturbed state, a steepest descent path (SDP) is followed on the PES, , until a local minimum on a new stable segment is reached. This new segment is then followed until it loses stability and so on. [Note that the SDP is defined to follow the negative of the local gradient at each point. It is traced out using a fourth-order Runge–Kutta method to integrate the system of first-order differential equations given by , where is held fixed.]

Fig. 7 presents the QD path of the nanoslab system in energy–strain space. Only stable segments are shown. The red segments correspond to the evolution of the system with loading when forced to transition from one local minimum to another by the quenched dynamics. In the figure, black vertical segments represent rolling down the PES at a point of instability (black dot) through an SDP to another minimum. This figure highlights the fact that unlike the QP path, which follows the lower envelope of the potential energy (Fig. 4), the QD path remains on a stable segment until a loss of stability is encountered. The QD path in (compressive) stress–strain space is shown in Fig. 8.

The sequence of configurations corresponding to the start and end of each of the stable segments constituting the QD path is presented in Fig. 9. The color coding of the atoms is the same as in Fig. 6. A close study of Fig. 9 reveals that the start and end configurations corresponding to each stable segment (red lines in Fig. 8) share the same symmetry group, whereas following the SDP to a minimum (black lines in Fig. 8) is associated with either symmetry breaking (observed in the initial stages of the QD path) or addition of symmetry (observed in the final stages of the QD path). This can be clearly observed in the first segment of the path. The structures representing the start and end of the first segment belong to the complete symmetry group (*D*_{2h} point group). Although not visible in the figure, a state of incipient slip exists in the ending structure of the first segment (second structure in Fig. 9) in the form of two X-shaped shear bands on the top and bottom of the nanoslab. At this point a small perturbation in the direction of the unstable eigenmode of the Hessian nucleates the X shear band at the top of the nanoslab at the expense of the other, thereby breaking the symmetry.

From the computed EM and the associated QP and QD paths (shown in Figs. 6 and 9) an unexpected but simple (“obvious” in hindsight) description of the system’s behavior emerges in terms of crystal rotation. Specifically, portions of the fcc crystal rotate by 90° to convert the original free surfaces (on the sides of the nanoslab) into free surfaces that are 12% lower in energy (29). This is confirmed by the energy difference computed between the unloaded initial configuration and the unloaded global minimum energy configuration . Examining Fig. 9 with this in mind, the QD path can be interpreted as the nucleation of a rotated grain followed by its growth until it encompasses the entire specimen (except for the boundaries where the applied constraints do not allow it). In contrast, the QP path in Fig. 6 can be interpreted as one where the system finds the most homogeneous mixture of the two crystal orientations consistent with the boundary conditions at a given value of *λ*. In particular, the system jumps back and forth between a contiguous single grain of the rotated crystal (such as is the case for the final two configurations shown in the third row of Fig. 6) and two grains of the rotated crystal separated by either a stacking fault or a small grain of the original unrotated crystal (such as is the case for the first two configurations in the fourth row of Fig. 6). It is this sort of systematic understanding of system response that becomes possible due to the EM-based approach. Contrast this with a single static or dynamic trajectory (as obtained from standard atomistic methods) that provides no insight into the overall response landscape.

Finally, we consider the DD regime in which the system is loaded at a finite rate as it would be in a real-world experiment. Standard MD simulations apply extremely high loading rates to overcome time scale limitations. However, many physical processes exhibit a sensitivity to loading rate which can lead to a situation where MD incorrectly predicts the system’s behavior. A recent example, relevant to the simulations presented here, are experiments showing strong rate dependence in the compression of copper nanopillars (32). Such rate dependence can be readily understood within the context of the EM as resulting from the availability of multiple solutions and the limited amount of time available for the system to sample the PES as it changes due to the evolving loading.

Consider the case where the loading parameter is a function of time, . Because at each value of loading, the EM includes a complete description of all available metastable states (basins) and the transition states connecting them (saddle points), the EM contains all information needed to predict the system response. [This is true as long as the loading rate is not too fast relative to the equilibration time of the system within each available metastable basin.] An example is given in Fig. 10 showing all available distinct (nonsymmetry-related) metastable states for the nickel nanoslab at . (Symmetry-related states can be automatically generated from these representatives and added to the set of available states.) To determine the connections between these states, we perform two steepest descent minimizations from each saddle point in the EM at , one in the positive direction of the unstable eigenvector and one in the negative direction. From these calculations we obtain the identities of the two basins which are connected by a transition pathway through the saddle point.

At sufficiently low temperature the transition rate from one basin to another across a saddle point can be estimated using harmonic transition state theory (hTST) (33), which locally approximates the PES near a basin or saddle point using a second-order approximation. The hTST transition rate has the form , where *ν* is an attempt frequency related to the normal mode frequencies at the basin and saddle, is the energy difference between the saddle and basin, *k* is Boltzmann’s constant, and *T* is the temperature. [ is directly computable from the EM, and the information needed to compute *ν* is available during the construction of the EM.] Once the transition rates are computed for all of the connected states at a given level of strain, kinetic Monte Carlo (KMC) (34) can be used to simulate the system’s dynamical trajectory and follow its sequence of transitions from state to state. To do so, one constructs a rate table ** R**, such that

*R*

_{ij}is the transition rate from state

*i*to state

*j*(for unconnected states

*i*and

*j*, ). The KMC algorithm proceeds as follows:

1) Given state

*i*and time*t*, compute a time to escape from state*i*according to the exponential distribution for a Poisson process (34), , where is the total escape rate out of state*i*, and*r*is a uniformly distributed random number in the range .2) Randomly choose one of the possible transitions based on their relative probability, .

3) Increment the clock, , move the system to the new state selected in step 2, and go to step 1.

At a constant level of strain, the KMC procedure follows the system as it repeatedly transitions between metastable states with probabilities commensurate with the transition rates between them. However, the standard KMC algorithm described above does not account for the fact that the transition rates change in time for the DD regime. [The transition rates are time dependent because the PES is evolving with the change in *λ* according to the applied loading rate.] Prados et al. (35) have extended the KMC algorithm to account for time-dependent transition rates, and here we adopt their methodology. Our procedure to compute the system’s behavior at any loading rate is as follows. First, we precompute the hTST transition rate tables at load intervals using the information stored in the EM, i.e., , , , …. Second, we choose a fixed loading rate , which determines the time dependence of the transition rates through their dependence on the load *λ*. Third, we apply the Prados et al. method to obtain a KMC trajectory at the specified loading rate:

1) Given state

*i*at time*t*, compute an escape time consistent with the evolving rate table by solving , where is the total escape rate out of state*i*at load*λ*, and*r*is a uniformly distributed random number in the range .2) Increment the clock, , and compute the rate table at the new time

*t*, , by interpolating between the precomputed rate tables.3) Follow the standard KMC procedure (step 2) to select one of the available states at time

*t*based on their relative probabilities. Move the system to the new state and go to step 1.

The time-dependent KMC procedure described above is repeated multiple times to obtain an ensemble of KMC trajectories. These trajectories are averaged to obtain the predicted load–deformation behavior for the system at the specified loading rate.

To illustrate this technique, the procedure was applied to the nickel nanoslab for two different loading rates. Fig. 11 shows, in blue, the average stress versus strain at a low strain rate of 1 s^{−1}. (The strain rate is related to the loading rate by , where is the undeformed height of the nanoslab.) The gray vertical line segments represent the standard deviation for 400 sampled KMC trajectories about their mean value. For comparison, the QP path is shown in red. This strain rate corresponds to loading the nanoslab by holding one end fixed and displacing the other end at a speed of 0.352 nm/s. At low strain rates, the PES evolves slowly compared with the transition rate between different states in the EM. This implies that the system is able to visit most of the available states at any given value of *λ*. However, it spends most of its time in low-energy states so that the resulting average behavior closely approximates the QP limit. For the finite, but low, loading rate shown in Fig. 11, the behavior follows that of the QP limit except near loads where the system transitions from one state to another. Here, the finite loading rate causes the system to transition between states at somewhat lower load levels than would otherwise occur.

Fig. 12 shows the mean evolution of the system at a high strain rate of 10^{8} s^{−1}, corresponding to a displacement speed of 3.52 m/s. The blue, gray, and red lines have the same meaning as in Fig. 11. The response is clearly different from that observed at the low strain rate. At this high strain rate, the PES evolves at a rate comparable to the rate of transition between states. As a result the dynamics is frustrated because the system does not have time to fully sample the available states before their properties significantly change due to the evolving load. This gives rise to a rate dependence in the system response. The system becomes “stuck” in a particular state despite the existence of other lower-energy states. As the strain rate approaches infinity the behavior is expected to approach that of the QD path. For the strain rate of 10^{8} s^{−1}, the stress–strain response in Fig. 12 shows a mixture of the QP and QD behavior.

Once the EM is constructed (which can be very costly but need only be done once), performing the KMC simulations is inexpensive. This means that many KMC simulations can be performed to obtain good statistics. Furthermore, this approach makes it possible to simulate processes at experimental loading rates which are inaccessible to dynamical methods such as MD. (Even for the tiny system studied here the low strain rate simulation is far beyond the capability of current-day computing.) It is also clear from the results above, that—together with our EM—KMC for time-dependent transition rates (35) is capable of predicting rate-dependent behavior. However, many aspects of the approach remain to be fully explored. An important issue is the completeness of the EM. We have limited the KMC calculations to a strain of because above this value the EM we generated is missing some key equilibrium segments and hence some important states and transitions. In principle, this issue is easily resolved by extending the computational effort used for the construction of the EM. In practice, further investigation and study are required to devise and develop a practical and robust computational scheme that will automatically generate important components of the full EM.

## Conclusions

A new technique is proposed for atomistic simulations using a BFB-generated EM of the system. The EM fully characterizes all possible responses of the system to applied loading. Three EM-based path selection and system evolution strategies, denoted QP, QD, and DD, are described that correspond to well-defined loading and environment conditions. This method provides an alternative to standard molecular statics and dynamics simulation strategies that consider only a single ill-defined trajectory over the PES. In situations where dynamics does not play an important role, the QP and QD paths provide the thermodynamic and quenched dynamics limits of the system response. Rate-dependent effects are included through the application of time-dependent KMC simulations as demonstrated in this work.

Similar to optimization techniques aiming to find the global minimum of a nonconvex function, the present method has its limitations. Clearly, the accuracy of the QP, QD, and DD paths obtained from the EM is dependent on its completeness and the possible existence of unexplored disconnected equilibrium segments. Determining the extent to which the EM provides a good representation of important deformation mechanisms for a given problem in nanoscale mechanics is thus an important issue that needs further study. Results from similar studies in structural mechanics are encouraging (14, 26, 28, 36).

Another area for future research is the application of this approach to spatial multiscale methods, such as the quasicontinuum (QC) method (37⇓⇓–40). An EM-based QC method would enable simulations at arbitrary loading rates of systems over length scales that are unavailable to full-resolution atomic calculations. This will make it possible to simultaneously span multiple length- and time-scales.

## Materials and Methods

The algorithms and pseudocodes implemented in the current version of the BFB code used in the construction of the EM are described in ref. 41. Within the BFB framework, a typical continuation–quenching step requires the system’s potential energy and its first and second derivatives with respect to the nuclear coordinates. This is computed by calling subroutines from the QC Method Package, wherein all of the atoms in the nanoslab are declared as nonlocal “representative atoms.” The basic theory of QC is outlined in refs. 39, 40 and the code is freely available for download at www.qcmethod.org (42). Computation of the equilibrium paths is automated using Perl scripts on a parallel cluster system taking advantage of the inherent parallelism in bifurcation problems. EM-based simulations are being made compatible with the newly developed application programming interface standard developed as part of the Knowledgebase of Interatomic Models (KIM) (http://openKIM.org) project (43).

## Acknowledgments

The authors acknowledge Tsvetanka Sendova for her contributions to the initial stages of this work. This work was partly supported by the National Science Foundation (NSF) Cyber-Enabled Discovery and Innovation (CDI) Grant Phy-0941493 (to R.S.E. and E.B.T.), NSF CAREER Grant CMMI-0746628 (to R.S.E.), Department of Energy Grant DE-SC0002085 (to E.B.T.), the Office of the Vice President for Research, University of Minnesota, and the Minnesota Supercomputing Institute.

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. E-mail: tadmor{at}aem.umn.edu.

Author contributions: R.S.E. and E.B.T. designed research; S.P. performed research; S.P., R.S.E., and E.B.T. analyzed data; and S.P., R.S.E., and E.B.T. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

## References

- ↵
- ↵
- ↵
- ↵
- Stillinger FH,
- Weber TA

- ↵
- Keller HB

- ↵
- Thompson JMT,
- Hunt GW

- ↵
- Thompson JMT

- ↵
- Thompson JMT,
- Hunt GW

- ↵
- Golubitsky M,
- Schaeffer DG

*Singularities and Groups in Bifurcation Theory*. Volume 1, Applied Mathematical Sciences (Springer, Berlin), 1st Ed, Vol. 51. - ↵
- Golubitsky M,
- Stewart I,
- Schaeffer DG

*Singularities and Groups in Bifurcation Theory*. Volume 2, Applied Mathematical Sciences (Springer, Berlin), 1st Ed., Vol. 69. - ↵
- Iooss G,
- Joseph DD

- ↵
- Govaerts WJ

- ↵
- Buffoni B,
- Toland J

- ↵
- Ikeda K,
- Murota K

- ↵
- Seydel R

- ↵
- ↵
- ↵
- Allgower EL,
- Georg K

- ↵
- Dankowicz H,
- Schilder F

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Wohlever JC

- ↵
- ↵
- ↵
- Stukowski A

*Model Simul Mater Sci Eng*18(1):015012. - ↵
- ↵
- ↵
- Sickafus KE,
- Kotomin EA

- Voter AF

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵
- Tadmor EB,
- Miller RE

- ↵
- Jusuf V

- ↵
- Tadmor EB,
- Miller RE

- ↵
- Tadmor EB,
- Elliott RS,
- Sethna JP,
- Miller RE,
- Becker CA

## Citation Manager Formats

## Article Classifications

- Physical Sciences
- Engineering