On a convergence analysis of a quasicontinuum method^{†}^{†}thanks: This work is supported by a grant of the Deutsche Forschungsgemeinschaft (DFG) SCHL 1706/21.
Abstract
In this article, we investigate a quasicontinuum method by means of analytical tools. More precisely, we compare a discretetocontinuum analysis of an atomistic onedimensional model problem with a corresponding quasicontinuum model. We consider next and nexttonearest neighbour interactions of LennardJones type and focus on the socalled quasinonlocal quasicontinuum approximation. Our analysis, which applies convergence techniques, shows that, in an elastic setting, minimizers and the minimal energies of the fully atomistic problem and its related quasicontinuum approximation have the same limiting behaviour as the number of atoms tends to infinity. In case of fracture this is in general not true. It turns out that the choice of representative atoms in the quasicontinuum approximation has an impact on the fracture energy and on the location of fracture. We give sufficient conditions for the choice of representative atoms such that, also in case of fracture, the minimal energies of the fully atomistic energy and its quasicontinuum approximation coincide in the limit and such that the crack is located in the atomistic region of the quasicontinuum model as desired.
Key words. Quasicontinuum method; atomistictocontinuum; –convergence; fracture.
AMS subject classifications. 49J45, 74R10, 74G15, 74G10, 74G65, 70C20
1 Introduction
The quasicontinuum (QC) method was introduced by Tadmor, Ortiz and Phillips [35] as a computational tool for atomistic simulations of crystalline solids at zero temperature. The key idea is to split the computational domain into regions where a very detailed (atomistic, nonlocal) description is needed and regions where a coarser (continuum, local) description is sufficient. The QCmethod and improvements of it are successfully used to study crystal defects such as dislocations, nanoindentations or cracks and their impact on the overall behaviour of the material, see e.g. [25].
There are various types of QCmethods: Some are formulated in an energy based framework, some in a force based framework; further, different couplings between the atomistic and continuum parts and different models in the continuum region are considered. In the previous decade, many articles related to the numerical analysis of such coupling methods were published. We refer to [15, 23] for recent overviews, in particular on the large literature including work on error analysis.
In this article, we consider a onedimensional problem and focus on the socalled quasinonlocal quasicontinuum (QNL) method, first proposed in [33]. The QNLmethod and further generalizations of it (see e.g. [16, 30]) are energybased QCmethods and are constructed to overcome asymmetries (so called ghostforces) at the atomistic/continuum interface which arise in the classical energy based QCmethod.
We are interested in an analytical approach in order to verify the QNLmethod as an appropriate mechanical model by means of a discretetocontinuum limit. This is embedded into the general aim of deriving continuum theories from atomistic models, see e.g. [3, Section 4.1], where also the need of a rigorous justification of QCmethods is addressed.
Our approach, announced in [34], is based on convergence, which is a notion for the convergence of variational problems, see e.g. [6]. We start with a onedimensional fully atomistic model problem which takes nearest and nexttonearest neighbour interactions into account. The limiting behaviour of the corresponding discrete model was analyzed by means of convergence techniques in [31] for a large number of atoms. In particular the limit and the first order limit are derived there, which take into account boundary layer effects.
From the fully atomistic model problem, we construct an approximation based on the QNLmethod. In particular, we keep the nearest and nexttonearest neighbour interactions in the atomistic (nonlocal) region and approximate the nexttonearest neighbour interactions in the continuum (local) region by certain nearest neighbour interactions as outlined below. Furthermore, we reduce the degree of freedom of the energy by fixing certain representative atoms and let the deformation of all atoms depend only on the deformation of these representative atoms.
It turns out that the choice of the representative atoms has a considerable impact on the validity of the QCmethod, see Theorem LABEL:5:cor, which is the main result of this work. This theorem asserts that the QCmethod is valid if the representative atoms are chosen in such a way that there is at least one nonrepresentative atom between two neighbouring representative atoms in the local region and in particular at the interface between the local and nonlocal regions. In Proposition LABEL:lemma:fracture, we prove that the mentioned sufficient condition on the choice of the representative atoms is indeed sharp by showing that in cases where the condition is not satisfied the limiting energy functional of the QCmethod does not have the same minima as the limiting energy of the fully atomistic model and thus should not be considered an appropriate approximation. This implies by means of analytical tools that in numerical simulations of fracture one has to make sure to pick a sufficiently large mesh in the continuum region and at the interface.
The outline of this article is as follows. In Section 2 we present the two discrete models, namely the fully atomistic and the quasicontinuum model, in detail.
In Sections 3 and 4 we investigate the limiting behaviour of the quasicontinuum energy functional by deriving the limits of zeroth and first order. It turns out that the limit of zeroth order of the fully atomistic and the quasicontinuum model coincide (Theorem LABEL:theorem:zero). If the boundary conditions are such that the specimen behaves elastically, we prove that both models also have the same limit of first order (Theorem LABEL:theorem:elastic).
If the boundary conditions are such that fracture occurs, the quasicontinuum approximation leads to a limit of first order (Theorem LABEL:theorem:fracture) that is in general different from the one obtained earlier for the fully atomistic model ([31], cf. Theorem LABEL:th:fracturennn). To compare the fully atomistic and the quasicontinuum model also in this regime, we analyze the limits of first order further in Section 5. As mentioned above, it turns out that if we use a sufficiently coarse mesh in the continuum region, the minimal energies of the two first order limits coincide (Theorem LABEL:5:cor). In fact we are able to show that in our particular model problem it is sufficient that the mesh size in the continuum region is at least twice the atomistic lattice distance. With this choice, fracture occurs always in the atomistic region as desired.
Furthermore, the convergence results imply, under suitable assumptions, a rate of convergence of the minimal energy of the quasicontinuum model to the minimal energy of the fully atomistic model (Theorem LABEL:th:limmin). Finally, we show that the condition on the mesh size is sharp. In Proposition LABEL:lemma:fracture, we provide examples where the corresponding limit has a different minimal energy and different minimizers than the fully atomistic system, which is due to poorly chosen meshes. This yields an analytical understanding of why meshes have to be chosen coarse enough in the continuum region.
Similar models as the one we consider here, were investigated previously in terms of numerical analysis. We refer especially to [14, 21, 26, 28, 29] where the QNL method is studied in one dimension. By proving notions of consistency and stability, those authors perform an error analysis in terms of the lattice spacing. To our knowledge, most of the results do not hold for “fractured” deformations. However, in [27] a Galerkin approximation of a discrete system is considered and error bounds are proven also for states with a single crack of which the position is prescribed. Recently, a different approach based on bifurcation theory is used in [22] to study the QCapproximation in the context of crack growth.
In [4], a different onedimensional atomisticcontinuum coupling method is investigated. Similar as in the QCmethod the domain is splitted in a discrete and a continuum region. In the discrete part the energy is given by nearest neighbour LennardJones interaction and in the continuum part by an integral functional with LennardJones energy density. It is shown that fracture is more favourable in the continuum than in the discrete region. To overcome this, the energy density of the continuum model is modified by introducing an additional term which depends on the lattice distance in the discrete region. Furthermore, in [5, p. 420] it is remarked that if the continuum model is replaced by a typical discretized version, the fracture is favourable in the discrete region. As mentioned above, we here treat a similar issue in the QNLmethod, see in particular Theorem LABEL:5:cor, Proposition LABEL:lemma:fracture.
The techniques of our analysis of the QNL method are related to earlier approaches based on convergence for the passage from discrete to continuum models in one dimension, see [8, 9, 10, 11, 12, 31, 32]; see also [18, 19] for a treatment of two dimensional models. Recently, convergence was used in [17] to study a QC approximation. In [17] a different atomistic model, namely a harmonic and defectfree crystal, is considered. Under general conditions it is shown that a quasicontinuum approximation based on summation rules has the same continuum limit as the fully atomistic system.
Common in all those works based on convergence is that primarily information about the global minimum and minimizers are obtained. Since atomistic solutions are not necessary global minimizers, it would be of interest to obtain also results for local minimizers, for instance in the lines of [7, 9]. In this article, we treat systems with nearest and nexttonearest neighbour interaction. A natural question is how the sufficient conditions on the choice of representative atoms change if we consider also interacting neighbours, . Therefore the corresponding fully atomistic model has first to be studied, which is part of ongoing research.
2 Setting of the Problem
First we describe our atomistic model problem which is the same as in [31]. We consider a onedimensional lattice given by with and interpret this as a chain of atoms. We denote by the deformation of the atoms from the reference configuration and write as shorthand. We identify such functions with their piecewise affine interpolations and define
The energy of a deformation is given by
where and are potentials of LennardJones type which will be specified in [LJ1]–[LJ4] below. Moreover, we impose boundary conditions on the first and last two atoms. For given we set
\hb@xt@.01(2.1) 
To consider only deformations which satisfy (LABEL:def:boundarycond), we define the functional
\hb@xt@.01(2.2) 
The goal is to solve the minimization problem
which we consider as our atomistic problem.
The idea of energy based quasicontinuum approximations is to replace the above minimization problem by a simpler one of which minimizers and minimal energies are good approximations of the ones for . Typically this new problem is obtained in two steps:

Define an energy where the long range (in our case nexttonearest neighbour) interactions are replaced by certain nearest neighbour interactions in some regions.

Reduce the degree of freedom by choosing a smaller set of admissible functions.
To obtain (a), the nexttonearest neighbour interactions are approximated as
see e.g. [28]. While this approximation turns out to be appropriate in the bulk, this is not the case close to surfaces, where the second neighbour interactions create boundary layers. This motivates to construct a quasinonlocal quasicontinuum model accordingly: For given let with . For we define the energy by using the above approximation for , cf. Fig. LABEL:fig:qcchain and keeping the atomistic descriptions elsewhere
Analogously to we define the functional
A crucial step for the following analysis is to rewrite the energy in a proper way. By defining
\hb@xt@.01(2.3) 
and , sometimes called CauchyBorn energy density (see [28]), we can write
\hb@xt@.01(2.4) 
for satisfying (LABEL:def:boundarycond). To emphasize the local structure of the continuum approximation, we rewrite the summation over the terms with in (LABEL:hnlqc) as an integral. To this end we use the fact that is constant on for and thus
Then
\hb@xt@.01(2.5) 
for satisfying (LABEL:def:boundarycond).
To obtain (b) we consider instead of the deformation of all atoms just the deformation of a possibly much smaller set of so called representative atoms (repatoms). We denote the set of repatoms by with and define
\hb@xt@.01(2.6) 
Since we are interested in the energy for deformations , we define
\hb@xt@.01(2.7) 
In the following sections we study as tends to infinity. Therefore, we will assume that is such that
\hb@xt@.01(2.8) 
Hence, in particular . The above assumption corresponds to the case that the size of the atomistic region becomes unbounded on a microscopic scale (i), but shrinks to a point on a macroscopic scale (ii). While assumption (i) is crucial, see also Remark LABEL:rem:elastic, the assumption (ii) can be easily replaced by , and . In this case the analysis is essentially the same, but in the case of fracture, see Theorem LABEL:theorem:fracture, one has to distinguish more cases. We assume (LABEL:ass:kn) (ii) here because it is the canonical case from a conceptual point of view. Otherwise the atomistic region and continuum region would be on the same macroscopic scale.
3 ZeroOrder Limit
In this section we derive the limit of the discrete energy (LABEL:def:hnlqct), which we refer to as zeroorder limit. This limit involves the convex and lower semicontinuous envelope of the effective potential energy which is already introduced in [11] defined by
\hb@xt@.01(3.1) 
We state the assumptions on the functions , and under which the following results are obtained.

] (strict convexity)

] (uniqueness of minimal energy configurations) For every such that we have where is defined as
\hb@xt@.01(3.2) This implies
\hb@xt@.01(3.3) 
] (regularity and behaviour at , ). be in on their domains such that on its domain. Let and . Moreover, we assume the following limiting behaviour
\hb@xt@.01(3.4) 
] (structure of , and ). , are such that there exists a convex function
\hb@xt@.01(3.5) and there exist constants such that
\hb@xt@.01(3.6) Further, there exist such that
\hb@xt@.01(3.7) and is strictly convex in on its domain for . Moreover, it holds and for all .
Remark 3.1
(a) The main examples we think of are LennardJones interactions, defined classically as
\hb@xt@.01(3.8) 
and . The calculations in [31, Remark 4.1] show that defined as above satisfy [LJ1]–[LJ4]. Another example of interatomic potentials which satisfy the above assumptions, see [31, Remark 4.1], are Morsepotentials, defined for as
\hb@xt@.01(3.9) 
(b) The assumptions [LJ1]–[LJ4] imply that . In particular, we have
\hb@xt@.01(3.10) 
(c) Note that [LJ4] and (LABEL:lim:j012) imply that either or that there exists such that or for . In [LJ3], we assume for simplicity. However, this could be dropped making suitable assumptions on in the following statements.
To define appropriate function spaces, we use a similar notation as in [8] and [31]. Let be a function with bounded variation. Then we say that if satisfies the Dirichlet boundary conditions and . To allow jumps in respectively , the boundary conditions are replaced by respectively in this case. Analogously, we define for special functions with bounded variations and the above boundary conditions. Let (or in ), then we denote by the jump set of in , and for we set . Moreover we denote by the singular part of the measure with respect to the Lebesgue measure.
Let us now state and prove the zerothorder limit of the functional . It turns out that the limiting functional is equal to the limit of the functional , cf. [31].
Theorem 3.2
Suppose [LJ1]–[LJ4] are satisfied and let . Let satisfy (LABEL:ass:kn) and let with be such that
\hb@xt@.01(3.11) 
Then the limit of defined in (LABEL:def:hnl) and of defined in (LABEL:def:hnlqct) with respect to the –topology is the functional defined by
on .
Proof. The result for follows from [31, Theorem 3.1]. Thus we prove the result for . The following compactness property and lower bound follow from [10, Theorem 3.7] and [11, Theorem 3.1]. For the readers convenience, we present direct proofs here.
Compactness. Let be a sequence with equibounded energy . The definition of and the properties of , imply that . Define the set . Next, we make use of the fact that are bounded from below and that the energy is equibounded. Moreover, we apply (LABEL:ass:psi2) and Jensen’s inequality to obtain
for some independent of . By (LABEL:ass:psi1), we have that for some constant independent of . Moreover, by using the boundary conditions, we obtain
Since , we obtain by the Poincaréinequality that is equibounded. Thus, we can extract a subsequence of which converges weakly to some , see [2, Theorem 3.23]. As argued in [31, Theorem 3.1], we have .
Liminf inequality. Let and be a sequence with equibounded energy which converges to in . The above compactness property and [2, Proposition 3.13] imply that converges to weakly in . By using [LJ3], [LJ4], we obtain for the recession function
with arbitrary. For every there exists such that for every . For large enough, we deduce from (LABEL:hnlqc:integral) by the definition of and [LJ4]
Note that by it follows for all , thus there exists such that
The last inequality is a direct implication of [2, Theorem 2.34], using that weakly converges to . By using that the righthand side above is finite only if , we obtain the liminf inequality from the arbitrariness of .
Limsup inequality. To show the existence of a recovery sequence, we first do not take the boundary conditions into account. Therefore, we define the functional by
For every we show existence of a sequence converging to in such that
\hb@xt@.01(3.12) 
As outlined in the proof of [10, Theorem 3.5] it is enough to show the above inequality for linear and for with a single jump: by density, this proves the statement for and the general estimate follows by relaxation arguments. Firstly, we consider functions with a single jump. Let with , and . By (LABEL:T_n) there exists with and such that for . We define now a sequence by
\hb@xt@.01(3.13) 
Obviously we have in . The functions are defined such that for and for all with . Using , (LABEL:T_n), [LJ3] and [LJ4] this implies
Now let for some . For every sequence satisfying (LABEL:T_n) we find a sequence of natural numbers such that
We define for every a set with , where such that there exist which satisfy
From we deduce and thus . Let us define such that and
By using and , we obtain
and thus in . Indeed, by , and , the last term tends to zero as . For the limsup inequality we argue similarly as in the case of a jump before. By definition, we have