Accuracy of Quasicontinuum Approximations Near Instabilities

Accuracy of Quasicontinuum Approximations
Near Instabilities

M. Dobson M. Luskin School of Mathematics, 206 Church St. SE, University of Minnesota, Minneapolis, MN 55455, USA,  and  C. Ortner Mathematical Institute, 24–29 St. Giles’, Oxford OX1 3LB, UK
July 5, 2019

The formation and motion of lattice defects such as cracks, dislocations, or grain boundaries, occurs when the lattice configuration loses stability, that is, when an eigenvalue of the Hessian of the lattice energy functional becomes negative. When the atomistic energy is approximated by a hybrid energy that couples atomistic and continuum models, the accuracy of the approximation can only be guaranteed near deformations where both the atomistic energy as well as the hybrid energy are stable. We propose, therefore, that it is essential for the evaluation of the predictive capability of atomistic-to-continuum coupling methods near instabilities that a theoretical analysis be performed, at least for some representative model problems, that determines whether the hybrid energies remain stable up to the onset of instability of the atomistic energy.

We formulate a one-dimensional model problem with nearest and next-nearest neighbor interactions and use rigorous analysis, asymptotic methods, and numerical experiments to obtain such sharp stability estimates for the basic conservative quasicontinuum (QC) approximations. Our results show that the consistent quasi-nonlocal QC approximation correctly reproduces the stability of the atomistic system, whereas the inconsistent energy-based QC approximation incorrectly predicts instability at a significantly reduced applied load that we describe by an analytic criterion in terms of the derivatives of the atomistic potential.

Key words and phrases:
atomistic-to-continuum coupling, defects, quasicontinuum method, sharp stability estimates
2000 Mathematics Subject Classification:
This work was supported in part by DMS-0757355, DMS-0811039, the Department of Energy under Award Number DE-FG02-05ER25706, the Institute for Mathematics and Its Applications, the University of Minnesota Supercomputing Institute, the University of Minnesota Doctoral Dissertation Fellowship, and the EPSRC critical mass programme “New Frontier in the Mathematics of Solids” (OxMoS). The project was initiated during a visit of ML to OxMoS.
.  To appear in J. Mech. Phys. Solids.

1. Introduction

An important application of atomistic-to-continuum coupling methods is the study of the quasistatic deformation of a crystal in order to model instabilities such as dislocation formation during nanoindentation, crack growth, or the deformation of grain boundaries [18]. In each of these applications, the quasistatic evolution provides an accurate approximation of the crystal deformation until the evolution approaches an unstable configuration. This occurs, for example, when a dislocation forms or moves or when a crack tip advances. The crystal will then typically undergo a dynamic process until it reaches a new stable configuration. In order to guarantee an accurate approximation of the entire quasistatic crystal deformation, up to the formation of an instability, it is crucial that the equilibrium in the atomistic/continuum hybrid method is stable whenever the corresponding atomistic equilibrium is. The purpose of this work is to investigate whether the quasicontinuum (QC) method has this property. In technical terms, this requires sharp estimates on the stability constant in the QC approximation.

The QC method is an atomistic-to-continuum coupling method that models the continuum region by using an energy density that exactly reproduces the lattice-based energy density at uniform strain (the Cauchy-Born rule) [18, 21, 25]. Several variants of the QC approximation have been proposed that differ in how the atomistic and continuum regions are coupled [3, 8, 18, 26]. In this paper, we present sharp stability analyses for the main examples of conservative QC approximations as a means to evaluate their relative predictive properties for defect formation and motion. Our sharp stability analyses compare the loads for which the atomistic energy is stable, that is, those loads where the Hessian of the atomistic energy is positive-definite, with the loads for which the QC energies are stable. It has previously been suggested and then observed in computational experiments that inconsistency at the atomistic-to-continuum interface can reduce the accuracy for computing a critical applied load [25, 18, 19, 8]. In this paper, we give an analytical method to estimate the error in the critical applied load by deriving stability criteria in terms of the derivatives of the atomistic potential.

Although we present our techniques in a precise mathematical format, we believe that these techniques can be utilized in a more informal way by computational scientists to quantitatively evaluate the predictive capability of other atomistic-to-continuum or multiphysics models as they arise. For example, our quantitative approach has the potential to estimate the reduced critical applied load in QC approximations such as the quasi-nonlocal QC approximation (QNL), that are consistent for next-nearest interactions but not for longer range interactions. Since the longer range interactions are generally weak, such an estimate may give an analytical basis to judging that the reduced critical applied load for QNL with finite range interactions is within an acceptable error tolerance.

The accuracy of various QC approximations and other atomistic-to-continuum coupling methods is currently being investigated by both computational experiments and numerical analysis  [1, 2, 4, 5, 7, 9, 11, 12, 14, 15, 16, 20, 23, 24]. The main issue that has been studied to date in the mathematical analyses is the rate of convergence with respect to the smoothness of the continuum solution (however, see [2, 7, 23] for analyses of the error of the QC solutions with respect to the atomistic solution, possibly containing defects). Some error estimates have been obtained that give theoretical justification for the accuracy of a QC approximation for all loads up to the critical atomistic load where the atomistic model loses stability [5, 23], but other error estimates that have been presented do not hold near the atomistic limit loads. It is important to understand whether the break-down of these error estimates is an artifact of the analysis, or whether the particular QC approximation actually does incorrectly predict an instability before the applied load has reached the correct limit load of the atomistic model.

Two key ingredients in any approximation error analysis are the consistency and stability of the approximation scheme. For energy minimization problems, consistency means that the truncation error for the equilibrium equations is small in a suitably chosen norm, and stability is usually understood as the positivity of the Hessian of the functional. For the highly non-convex problems we consider here, stability must necessarily be a local property: The configuration space can be divided into stable and unstable regions, and the question we ask is whether the stability regions of different QC approximations approximate the stability region of the full atomistic model in a way that can be controlled in the setup of the method (for example, by a judicious choice of the atomistic region).

In this work, we initiate such a systematic study of the stability of QC approximations. In the present paper, we investigate conservative QC approximations, that is, QC approximations which are formulated in terms of the minimization of an energy functional. In a companion paper [6], we study the stability of a force-based approach to atomistic-to-continuum coupling that is nonconservative.

In computational experiments, one often studies the evolution of a system under incremental loading. There, the critical load at which the system “jumps” from one energy well to another is often the goal of the computation. Thus, we will also study the effect of the “stability error” on the error in the critical load.

We will formulate a simple model problem, a one dimensional periodic atomistic chain with pairwise next-nearest neighbour interactions of Lennard-Jones type, for which we can analyze the issues layed out in the previous paragraphs. It is well known that the uniform configuration is stable only up to a critical value of the tensile strain (fracture). We use analytic, asymptotic, and numerical approaches to obtain sharp results for the stability of different QC approximations when applied to this simple model.

In Section 2, we describe the model and the various QC approximations that we will analyze. In Section 4, we study the stability of the atomistic model as well as two consistent QC approximations: the local QC approximation (QCL) and the quasi-nonlocal QC approximation (QNL). We prove that the critical applied strains for both of these approximations are equal to the critical applied strain for the atomistic model, up to second-order in the atomistic spacing.

A similar analysis for the inconsistent QCE approximation is more difficult because the uniform configuration is not an equilibrium. Thus, in Section 5, we construct a first-order correction of the uniform configuration to approximate an equilibrium configuration, and we study the positive-definiteness of the Hessian for the linearization about this configuration. We explicitly construct a test function with strain concentrated in the atomistic-continuum interface that is unstable for applied strains bounded well away from the atomistic critical applied strain.

In Section 6, we analyze the accuracy in predicting the critical strain for onset of instability. For the QCL and QNL approximations, this involves comparing the effect of the difference between their modified stability criteria and that of the atomistic model. For QCE, since the solution to the nonlinear equilibrium equations are non-trivial, we provide computational results in addition to an analysis of the critical QCE strain predicted by the approximations derived in Section 5.

2. The atomistic and quasicontinuum models

2.1. The atomistic model problem

Figure 1. Lennard-Jones type interaction potential. The bond length is the turning point between the convex and concave regions of .

Suppose that the infinite lattice is deformed uniformly into the lattice where is the macroscopic deformation gradient and where scales the reference atomic spacing, that is,

We admit -periodic perturbations from the uniformly deformed lattice More precisely, for fixed , we admit deformations from the space

where is the space of -periodic displacements with zero mean,

We set throughout so that the reference length of the periodic domain is fixed. Even though the energies and forces we will introduce are well-defined for all -periodic displacements, we require that they have zero mean in order to obtain locally unique solutions to the equilibrium equations. These zero mean constraints are an artifact of our periodic boundary conditions and are similarly used in the analysis of continuum problems with periodic boundary conditions.

We assume that the stored energy per period of a deformation is given by a next-nearest neighbour pair interaction model,


where is the backward difference

and where is a Lennard-Jones type interaction potential satisfying (see also Figure 1)

  1. ,

  2. there exists such that is convex in and concave in , and

  3. rapidly as , for .

We have used the scaled interaction potential, in the definition of the stored energy, to obtain a continuum limit as Assumptions (i) and (ii) are used throughout our analysis, while assumption (iii) serves primarily to motivate that next-nearest neighbour interaction terms are typically dominated by nearest-neighbour terms. Note, however, that even with assumption (iii), the relative size of next-nearest and nearest neighbour interactions is comparable when strains approach .

We denote the first variation of the energy functional, at a deformation by

for . In the absence of external forces, the uniformly deformed lattice is an equilibrium of the atomistic energy under perturbations from , that is,


We identify the stability of with linear stability under perturbations from the space . To make this precise, we denote the second variation of the energy functional, evaluated at a deformation , by


for all . The matrix is the Hessian for the energy functional. We say that the equilibrium is stable for the atomistic model if this Hessian, evaluated at , is positive definite on the subspace of zero mean displacements, or equivalently, if


In Section 3, Definition 3, we extend this definition of stability to the various QC approximation and their equilibria.

Note that if , then and for all . Therefore, upon defining the quantities

we can rewrite (3) as follows


(We will use later.) The quantities and will play a prominent role in the analysis of the stability of the atomistic model and its QC approximations and describe the strength of the nearest neighbor and next-nearest neighbor interactions, respectively. We similarly define the quantities for all and for all . For most realistic interaction potentials the second-nearest neighbour coefficient is non-positive, except in the case of extreme compression (see Figure 1). Therefore, in order to avoid having to distinguish several cases, we will assume throughout our analysis that . In this case, property (ii) of the interaction potential shows that .

We also note that, for , both and are understood as -periodic chains, that is, where the centered second difference is defined by

For , we also define the weighted -norms

as well as the weighted -inner product

2.2. The local QC approximation (QCL)

Before we introduce different flavors of QC approximations, we note that we can rewrite the atomistic energy as a sum over the contributions from each atom,

If is “smooth,” i.e., varies slowly, then where

and where is the so-called Cauchy-Born stored energy density. In this case, we may expect that the atomistic model is accurately represented by the local QC (or continuum) model


The main feature of this continuum model is that the next-nearest neighbour interactions have been replaced by nearest neighbour interactions, thus yielding a model with more locality. Such a model can subsequently be coarse-grained (i.e., degrees of freedom are removed) which yields efficient numerical methods.

2.3. The energy-based QC approximation (QCE)

If is “smooth” in the majority of the computational domain, but not in a small neighbourhood, say, , where , then we can obtain sufficient accuracy and efficiency by coupling the atomistic model to the local QC model by simply choosing energy contributions in the atomistic region and in the continuum region . This approximation of the atomistic energy is often called the energy based QC approximation [21] and yields the energy functional

It is now well-understood [3, 4, 5, 8, 25] that the QCE approximation exhibits an inconsistency (“ghost force”) near the interface, which is displayed in the fact that The first remedy of this lack of consistency was the ghost force correction scheme [25] which eventually led to the derivation of the force-based QC approximation [3] and which we analyze in [6] and [7].

2.4. Quasi-nonlocal coupling (QNL)

An alternative approach was suggested in [26], which requires a modification of the energy at the interface. This idea is best understood in terms of interactions rather than energy contributions of individual atoms (see also [8] where this has been extended to longer range interactions). The nearest neighbour interactions are left unchanged. A next-nearest neighbour interaction is left unchanged if at least one of the atoms belong to the atomistic region and is replaced by a Cauchy–Born approximation,

if both atoms belong to the continuum region. This idea leads to the energy functional

where and are modified atomistic and continuum regions. The QNL approximation is consistent, that is, is an equilibrium of the QNL energy functional. The label QNL comes from the original intuition of considering interfacial atoms as quasi-nonlocal, i.e., they interact by different rules with atoms in the atomistic and continuum regions.

3. Stability of Quasicontinuum Approximations: Summary

In this section, we briefly summarize the our main results.

We begin by giving a careful definition of a notion of stability. Our condition is slightly stronger than local minimality, which is the natural concept of stability in statics. However, an analysis of local minimality alone is usually not tractable. Moreover, for the deformations that we consider, our definition is in fact sufficiently general.

Definition 1 (Stable Equilibrium). Let . We say that is a stable equilibrium of if is twice differentiable at and the following conditions hold:

  • for all ,

  • for all .

If only (i) holds, then we call a critical point of .

We focus on the deformation and ask for which macroscopic strains this deformation is a stable equilibrium. We know from (2) that is a critical point of the atomistic energy , and it is easy to see that is also a critical point of the QCL energy and of the QNL energy . Our analysis in Section 4 gives the following conditions under which is stable:

  1. is a stable equilibrium of if and only if ;

  2. is a stable equilibrium of if and only if ;

  3. is a stable equilibrium of if and only if

where we recall that is the continuum elastic modulus for the Cauchy–Born stored energy function . Points 1., 2., and 3. are established, respectively, in Propositions 4.1, 4.2, and 4.3.

If we envision a quasistatic process in which is slowly increased, then we may wish to find the critical strain at which is no longer a stable equilibrium (fracture instability). If we denote the critical strains in the atomistic, QCL, and QNL models, respectively, by , and , then 1.–3. imply that (cf. Section 6)

For the QCE approximation defined in Section 2.3, the situation is more complicated. The occurrence of a “ghost force” in the QCE model implies that is not a critical point of , and consequently, we will need to analyze the stability of the second variation where is an appropriately chosen equilibrium of . Since solves a nonlinear equation, we will replace it by an approximate equilibrium in our analysis in Section 5 where we obtain the following (simplified) result:

  1. For to be a stable equilibrium of it is necessary that

    where is assumed to be small.

We remark that 4. gives only a necessary but not a sufficient condition for stability of the QCE equilibrium , which, moreover, depend on assumptions on the parameter . We refer to Remark 5 for a careful discussion of the role of .

If we let denote the critical strain at which 4. fails (ignoring the term), then we obtain

which suggests that the QCE method is unable to predict the onset of fracture instability accurately. In Section 6, we confirm this asymptotic prediction with numerical experiments.

We have shown in [17] that the stability properties of the ghost force correction scheme (GFC) can be understood for uniaxial tensile loading by considering the stability of the QC energy


We note that so is an equilibrium of the energy under perturbations from We can therefore analyze the stability of at by studying the Hessian We show in Remark 5 in Section 2.3 that

  1. is a stable equilibrium of if and only if where

We give analytic and computation results in Sections 5 and 6 showing that the ghost force correction scheme can be expected to improve the accuracy of the computation of the critical strain by the QCE method, that is,

where is the critical strain at which is no longer positive definite, but the error in computing the critical strain by the GFC scheme is still

4. Sharp Stability Analysis of Consistent QC Approximations

In this section, we analyze the stability of the atomistic model and two consistent QC approximations: the local QC approximation and the quasi-nonlocal QC approximation. In each case, we will give precise conditions on under which is stable in the respective approximation. The inconsistent energy-based QC approximation (QCE) is analyzed in Section 5. The corresponding result for QCE is less exact than for QCL and QNL, but shows that there is a much more significant loss of stability.

4.1. Atomistic model

Recalling the representation of from (5) and noting that


we obtain


To quantify the influence of the strain gradient term, we define

Since is periodic, it follows that has zero mean. In this case, the eigenvalue is known to be attained by the eigenfunction and is given by [27, Exercise 13.9]


Since as , it follows that as . Thus, we obtain the following stability result for the atomistic model.

Proposition 2. Suppose . Then is stable in the atomistic model if and only if , where is the eigenvalue defined in (10).


By the definition of , and using (9), we have

4.2. The Local QC approximation

The equilibrium system, in variational form, for the QCL approximation is

Since has zero mean, it follows that is a critical point of for all . The second variation of the local QC energy, evaluated at , is given by

Thus, recalling our definition of stability from Section 2.1, we obtain the following result.

Proposition 3. The deformation is a stable equilibrium of the local QC approximation if and only if .

Comparing Proposition 4.2 with Proposition 4.1 we see a first discrepancy, albeit small, between the stability of the full atomistic model and the local QC approximation (or the Cauchy–Born approximation). In Section 6 we will show that this leads to a negligible error in the computed critical load.

4.3. Quasi-nonlocal coupling

By the construction of the QNL coupling rule at the interface, the deformation is an equilibrium of [26]. The second variation of evaluated at is given by

We use (8) to rewrite the second group on the right-hand side (the nonlocal interactions) in the form

to obtain

Except in the case , it now follows immediately that is stable in the QNL approximation if and only if .

Proposition 4. Suppose that and that , then is stable in the QNL approximation if and only if .

5. Stability Analysis of the Energy-based QC approximation

We will explain in Remark 5 that there exists such that

However, is not a critical point of so we must be careful in extending the previous definition of stability to the QCE approximation. We cannot simply consider the positive-definiteness of Instead, we analyze the second variation where solves the QCE equilibrium equations

or equivalently


We will see that, when the second-neighbour interactions are small compared with the first neighbour interactions (which we make precise in Lemma 5), there is a locally unique solution of the equilibrium equations, which is the correct QCE counterpart of . We will then derive a stability criterion for the equilibrium deformation

In Proposition 5 below, we derive an upper bound for the coercivity of with respect to the norm Even though the derivation of this upper bound is only rigorous for strains bounded away from the atomistic critical strain, it clearly identifies a source of instability that cannot be found by analyzing, for example, . Moreover, we will present numerical experiments in Section 6 showing that the critical strain predicted in our following analysis gives a remarkably accurate approximation to actual QCE critical strain.

In Section 6, we consider the critical strain for each approximation, namely the point at which the appropriate equilibrium deformation (either or ) becomes unstable. We will see later in this section, as well as in Section 6 that predicting the loss of stability for the QCE approximation using greatly underestimates the error in approximating the atomistic critical strain by the QCE critical strain.

Due to the nonlinearity and nonlocality of the interaction law, we cannot compute explicitly. Instead, we will construct an approximation which is accurate whenever second-neighbour terms are dominated by first-neighbour terms. In the following paragraphs, we first present a semi-heuristic construction, motivated by the analysis in [4], and then a rigorous approximation result, the proof of which is given in Appendix B.

In (26) in the appendix, we provide an explicit representation of . Inserting , we obtain a variational representation of the atomistic-to-continuum interfacial truncation error terms that are often dubbed “ghost forces,”




We note that (12) makes our claim precise that is not a critical point of .

Motivated by property (iii) of the interaction potential , we will assume that the parameters

are small, and construct an approximation for which is asymptotically of second order as . Although such an approximation will not be valid near the critical strain for the QCE approximation, it will give us a rough impression how the inconsistency affects the stability of the system. We note that is scale invariant since we used a scaled interaction potential, in our definition of the stored energy (1).

A non-dimensionalization of (12) shows that . If is small, then we can linearize (11) about and find the first-order correction , which is given by


We note that this linear system is precisely the one analyzed in detail in [4]. However, instead of using the implicit representation of obtained there, we use the assumption that is small to simplify (14) further and obtain a more explicit approximation.

Writing out the bilinear form explicitly (using (28) as a starting point) gives


where we have only displayed the terms in the right half of the domain and indicated the terms in the left half by dots. Ignoring all terms involving , which are of order relative to the remaining terms, we arrive at the following approximation of (14):

the solution of which is given by

The following lemma makes this approximation rigorous. A complete proof is given in Appendix B.

Lemma 5. If and are sufficiently small, then there exists a (locally unique) solution of (11) such that

where may depend on (and its derivatives) and on , but is independent of .

From now on, we will also assume that is small, and combine the three small parameters into a single parameter

We will neglect all terms which are of order . A careful discussion of the parameter and the validity of the asymptotic analysis is given in Remark 5.

In the following, we will again only show terms appearing on the right half of the domain. Our goal in the remainder of this section is to obtain an estimate for the smallest eigenvalue of . Using (28), we can represent as

We expand all terms containing and neglect all terms which are of order relative to which is the order of magnitude of the coefficient of the diagonal term of For example, we have, for some ,

as . Thus, the perturbation of a second-neighbour term will not affect our final result. On the other hand, expanding a nearest neighbour term gives

as . Proceeding in the same fashion for the remaining terms, we arrive at