# A Posteriori Error Control for a Quasicontinuum Approximation of a Periodic Chain

## Abstract

We consider a 1D periodic atomistic model, for which we formulate and analyze an adaptive variant of a quasicontinuum method. We establish a posteriori error estimates for the energy norm and for the energy, based on a posteriori residual and stability estimates. We formulate adaptive mesh refinement algorithms based on these error estimators. Our numerical experiments indicate optimal convergence rates of these algorithms.

CO was supported by an EPSRC Grant “Analysis of Atomistic-to-Continuum Coupling Methods”. HW was supported by China Scholarship Council and University of Oxford’s “Chinese Ministry of Education-University of Oxford Scholarship” and Sichuan University’s National high Level University Construction Programme

## 1Introduction

Quasicontinuum (QC) methods, or more generally, atomistic-to-continuum coupling (a/c) methods, are a class of multiscale methods for coupling an atomistic model of a solid with a continuum model. These methods have been widely employed in atomistic simulations, where a fully atomistic model would result in prohibitive computational cost but is required in certain regions of interest (typically neighbourhoods of crystal defects) to achieve the desired degree of accuracy. A continuum model is applied to reduce the cost of the computation of the elastic far-fields. We refer to the recent review articles [13] for introductions to QC methods.

In this paper we present an a posteriori error analysis of a simplified variant of the energy-based QC approximation proposed in [26].

Considerable effort has been devoted to the a priori error analysis of QC methods [11]. A posteriori error control has received comparatively little attention: Arndt and Luskin use a goal-oriented approach for the a posteriori error estimates for the QC approximation of a Frenkel-Kontorova model [2]. The estimates are used to optimize the choice of the atomistic region and the finite element mesh in the continuum region. Prudhomme et al. [24] also use the goal-oriented approach to provide a posteriori error control for the original energy-based QC approximation [17], which is inconsistent. Both of these approaches require the use of the solutions of dual problems. In [22] a posteriori error bounds for an energy norm are derived, using a similar approach as in the present work. However, the QC method analyzed in [22] does not contain an approximation of the stored energy and its 2D/3D version is therefore not a computationally efficient method.

In the present work, we use the residual-based approach [28], which is well established in the finite element approximation of partial differential equations, to derive a posteriori error bounds for the energy-norm and for the energy itself. What distinguishes our setting from the classical one, are two particular features: 1. the model approximation (“variational crime”) is fundamentally different than quadrature approximations; 2. it is possible to adapt the mesh to a degree that the residual vanishes exactly in certain regions of interest (the atomistic region). Our work extends [18], which considers a simplified setting.

We do not tackle the question of model and mesh adaptivity for defect nucleation, but focus only on automatically choosing the *size* of the atomistic region and the finite element mesh in the continuum region. In many applications, a defect is inserted into the crystal before the computation, or it is known a priori where a defect will nucleate. An extension of our work to include defect nucleation would be valuable, but cannot be pursued in the simplified 1D setting we are considering here.

### 1.1Outline

In Section 2 we introduce a 1D atomistic model problem, which mimics the behaviour of crystal defects in 2D/3D. Moreover, we review the construction of a QC method to efficiently approximate its solutions, and introduce the notation that will be used throughout the paper.

In Section 3, we derive a residual estimate for the QC method in a discrete negative Sobolev norm. In Section Section 4, we present an a posteriori stability result. In Section 5, we combine the residual estimate and the stability result to give a posteriori error estimates for the deformation gradient and for the energy.

In Section 6, we describe three mesh refinement algorithms based on our a posteriori error analysis and on previous a priori error estimates, and present a numerical example to illustrate the performance of these algorithms.

## 2Model Problem and QC Approximation

### 2.1Atomistic Model

Following previous works [7] we formulate a model problem in a discrete periodic domain containing atoms, where . Let denote a macroscopic stretch and the lattice spacing, both of which we fix throughout. We define the displacement and deformation spaces, respectively, by

In particular, we observe that if and only if for some .

The *stored energy* (per period) of an admissible deformation is given by

where and we note that is a second-neighbour difference, and where is a Lennard-Jones type interaction potential. We assume throughout that there exists such that is convex in and concave in .

Given a periodic dead load , we define the *external energy* (per period) by

where . Thus, the *total energy* (per period) under a deformation is given by

We wish to compute

where denotes the set of local minimizers.

1. Since the internal energy is translation invariant, our choice (instead of the more common constraint ) does not alter the problem but simplifies the treatment of external forces in Sections Section 3.2 and Section 3.4.

2. The external energy can be thought of as the linearisation of a nonlinear external energy about the state .

3. We have included the external forces primarily to render the problem non-trivial. In realistic applications in 2D/3D, defects or forces applied on a boundary are the cause of deformation of the crystal, however, in 1D such complex behaviour cannot be described directly.

### 2.2Notation for lattice functions

Throughout, we identify a lattice function with its continuous and piecewise affine interpolant with nodal values . Vice-versa, if then we always identify it with its vectors of nodal values .

Let be a subset of . For a vector , we define

If the label is omitted, then we understand this to mean .

We define the first discrete derivatives . It is then straightforward to see that, if at the same time we identify with the pointwise derivative, then .

We equip the space with the discrete Sobolev norm (recall that )

The norm on the dual is defined by

### 2.3Finite element notation

To construct the QC approximation in the next section, we first define some convenient notation. Let denote the computational cell. We choose an atomistic region , where atomistic accuracy is required, and we define the continuum region by . We will also use the periodic extension of , denoted by .

We assume throughout that is an open interval with . All our results (with the exception of Section 3.4) can be extended without difficulty to the case when consists of a finite union of open intervals.

Let be a partition of into closed intervals with , where are the nodes of the partition. We assume, without loss of generality, that is the left-most node and the right-most nodes in the interval .

The length of an element is denoted by . The space of continuous piecewise affine functions with respect to the partition is denoted by .

We assume throughout that the partition has the following properties:

is periodic: there exists such that for all .

has atomistic resolution in : .

The a/c interface points are finite element nodes: . In particular, each element belongs entirely to either the atomistic or continuum region.

If then .

Property (T4) is not strictly required, but simplifies the analysis and is not a significant restriction. Note also that we have not required (except in the atomistic region) that finite element nodes must be positioned on atomic sites. Although not necessary in 1D, it somewhat simplifies mesh generation, and to some extend mimics the fact that element edges or faces in 2D/3D cannot normally be aligned with the underlying crystal lattice.

The finite element displacement and deformation spaces are defined, respectively, by

For , we define the interpolation operator by

We note that .

For future reference, we also define the micro-elements for . Analogously, we define to be the nodal interpolant with respect to the atomistic grid. We will require this interpolant since the mesh nodes do not necessarily coincide with lattice sites.

### 2.4QC Approximation

The QC approximation we analyze in this paper is the 1D variant of the ACC method described in [26]. (An earlier variant of the idea was described in [18] and a similar construction in [10]. We focus on the formulation proposed in [26] since it can be readily generalised to 2D.)

The idea of the ACC method is based on the splitting of interaction bonds. A bond is an open interval for and (since we consider only first and second neighbour interactions). Since our computational domain is the set of bonds over which the atomistic energy is defined is given by

For each bond we define .

For any open interval (e.g., for a bond) with length we define the finite difference operator

Note that, with this notation, we can rewrite .

For each bond and deformation field we define its atomistic and continuum energy contributions to the stored energy, respectively, by

If then is ill-defined; in that case we define . The QC energy (the ACC energy of [26]) of a deformation field is now defined by

The key property of this definition is that it satisfies the patch test [26],

The *external energy* (per period) is given by

Note that, in the atomistic region, this reduces to the same form as . The *total energy* (per period) of a deformation is then given by

In the QC approximation we seek

It is initially not obvious why the formulation should reduce the complexity of the computation of over that of , since is still written as a summation over all bonds. However, one can readily check (see [26] for the details), that

where is the Cauchy–Born stored energy function. This formulation requires only a sum over all bonds within the atomistic region, and a standard finite element assembly procedure to compute the energy contribution from the continuum region. Thus, the evaluation of has only a computational complexity proportional to .

Analogous results hold in 2D [26]; a more involved assembly procedure is required to make the ACC method efficient in 3D [25].

## 3Residual Analysis

### 3.1The residual

Let be a solution of the atomistic problem . It is straightforward to see that, if , then is differentiable at and hence the first order optimality condition for is satisfied:

Similarly, if is a solution of the QC problem with on , then it satisfies the corresponding first order optimality condition

Let be a solution of , then we define its residual by

Using we can rewrite this as

where we call the internal residual and the external residual. We will bound them separately in the next two sections.

### 3.2Estimate for the internal residual

In this section, we analyze the internal residual . To state the main theorem, we define the index set of all nodes in the continuum region and a/c interface (recall that is closed),

* Let such that and , then *

1. The expressions for are reminiscent of the flux (or stress) jump terms that occur in the classical residual error analysis for partial differential equations. The origin in our case, is somewhat different however, and results only from the model approximation and not the finite element coarsening.

2. With some additional work, the form of the estimates can be turned into element contributions and further simplified by computing more explicit representations.

Let . From and we obtain

We subtract and add the terms

to split into three components,

We will show that and estimate .

For this follows simply from the fact that in and hence for all such that .

For , following the computations in [26], we obtain

where . Since for all and since is constant on each element it follows that .

Finally, we turn to the analysis of , where we define

Using the fact that we obtain

If , then and , and hence . Similarly, if , then and and hence . Thus, we observe that only bonds crossing continuum element boundaries, or the atomistic/continuum interface, contribute to the residual. These are precisely the points contained in . In particular, we obtain

where we used the fact that no bond can cross more than one point due to our assumption that all elements have at least length .

From the definition of , and applying two Cauchy–Schwarz inequalities, it is straightforward to estimate

and after applying another Cauchy–Schwarz inequality,

where is defined in .

Combing our foregoing estimates, we arrive at

and we are only left to estimate the sums involving the test function.

To that end, we simply note that, due to (T4), for any fixed point , the maximal number of bonds appearing in the sum on the left-hand side below and crossing is three; hence,

This concludes the proof.

### 3.3Estimate of the external residual

We now turn to the estimate of the residual of the external energy, which was defined as

We outline the key points of the argument for estimating , before stating the result.

We define a slightly extended continuum region,

then in , and therefore

The three terms can be estimated using standard interpolation error results, hence we only give a brief outline of the proof of the resulting bound.

* Let , then where the error due to external forces and the “quadrature error” are, respectively, defined as follows: for , *

1. Note that there is an error contribution from the atomistic region, due to the fact that in the elements touching the a/c interface, the “quadrature” approximation of the external forces is not exact. For the purpose of mesh refinement, we count this error towards the neighbouring elements in the continuum region.

2. An alternative residual estimate that does not use , but only the discrete setting, is presented in [29]. This requires a much more involved argument.

From we obtain

Applying standard interpolation error estimates (see, e.g., [5]) on elements , we obtain

Summing over all elements, applying the Cauchy–Schwarz inequality, and defining for ,

We now use the estimates (which exploit the fact that is the -orthogonal projection of onto piecewise constants)

to deduce that

The result follows by splitting the norms inside the brackets over elements.

### 3.4External residual estimate for singular forces

In our numerical experiments in Section 6 we shall employ an external force that behaves like near (we use the “singularity” in the force to mimic a defect). Let us suppose that we also have and near the origin. We now give a formal motivation why the quadrature estimates employed in Proposition ? are inadequate in this situation.

Applying the quadrature estimates to such a force field, we obtain

We notice that, for near the origin, . Moreover, the quadrature estimate is and cannot be controlled. By contrast,

from which we conclude that is dominated by , but that is itself dominated by .

The origin of this undesirable effect is the (ab-)use of the Poincaré inequality in the proof of Proposition ?. In the remainder of this section, we shall remedy the situation by replacing the standard Poincaré inequality with a weighted variant. This approach is inspired by [20].

* Let be defined as in Section Section 3.3, and let , then *

We begin by noting that

Hence, we can estimate

Since