ERROR ESTIMATION AND ATOMISTIC-CONTINUUM ADAPTIVITY
FOR THE QUASICONTINUUM APPROXIMATION
OF A FRENKEL-KONTOROVA MODEL
MARCEL ARNDT AND MITCHELL LUSKIN
Abstract. We propose and analyze a goal-oriented a posteriori error estimator for the atomistic-
continuum modeling error in the quasicontinuum method. Based on this error estimator, we develop
an algorithm which adaptively determines the atomistic and continuum regions to compute a quan-
tity of interest to within a given tolerance. We apply the algorithm to the computation of the
structure of a crystallographic defect described by a Frenkel-Kontorova model and present the re-
sults of numerical experiments. The numerical results show that our method gives an efficient
estimate of the error and a nearly optimal atomistic-continuum modeling strategy.
1. Introduction
The quasicontinuum (QC) method [22, 23, 24] has been successfully used to efficiently cou-
ple atomistic and continuum models for crystalline solids and offers the possibility of computing
mesoscale or macroscale properties by a nearly minimal number of degrees of freedom. Accurate
modeling requires that an atomistic model be used in regions with highly non-uniform deformations
such as around dislocations, whereas a continuum model can be used in regions with nearly uniform
deformations to reduce the number of degrees of freedom.
It is usually not known a priori which regions of some specimen undergo uniform deformations
and which do not, so a posteriori error estimation is important for the design of efficient numerical
approximations by the quasicontinuum method. Since the purpose of a computation is often to
obtain the value of a (usually local) quantity of interest to a desired error tolerance rather than to
obtain a solution to a desired error tolerance for a global norm, there has been great interest in
the development of goal-oriented error estimators for many problems. They are based on duality
techniques and have been developed and used to adaptively refine finite element approximations of
continuum problems [1, 3] and to study and control modeling error [19].
In this paper, we extend this approach to develop an a posteriori error estimator for the qua-
sicontinuum method which quantifies the atomistic-continuum modeling error for a goal function
and allows for an adaptive decision about which regions can be accurately modeled as a continuum
and which regions need to be modeled atomistically. Methods to determine the optimal mesh size
within the continuum region will be studied in a forthcoming paper.
Crystallographic defects [5] provide a challenge to validate atomistic-continuum error estimators
and adaptivity. No such error estimators and adaptive methods currently exist for fully three-
dimensional crystals. As a step in this direction, we develop a rigorous theory for a simple one-
dimensional atomistic model for a defect that is a modification of the Frenkel-Kontorova model [15].
We add next-nearest-neighbor harmonic interactions between the atoms to the nearest-neighbor
harmonic interactions between the atoms in the classical Frenkel-Kontorova model.
A priori analyses for various quasicontinuum approximations have been given in [4, 8, 9, 10,
12, 13, 14, 21]. An a posteriori analysis for a slightly different one-dimensional quasicontinuum
http://arxiv.org/abs/0704.1924v1
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 2
approximation is given in [20]. The development and application of a goal-oriented error estimator
for mesh coarsening in a two-dimensional quasicontinuum method is reported in [17, 18].
Let us mention that the continuum model used in the QC method, which coincides with the
model obtained by the classical thermodynamic limit, is by far not the only reasonable continuum
model to use. A method to derive continuum models which approximate atomistic models up to
an arbitrarily high order has been proposed in [2].
The paper is organized as follows. In Section 2, we give a general formulation of the one-
dimensional quasicontinuum approximation [23] that includes not only two-body and three-body
potentials, but also many body potentials such as the embedded atom potential [6, 7]. In Section 3,
we describe our extension of the Frenkel-Kontorova model and its quasicontinuum approximation.
In Section 4, we introduce the primal and dual problems for our model and formulate our approach
to goal-oriented error estimation.
Next, in Section 5 we extend the approach in [16] to develop an error estimator for atomistic-
continuum modeling. This first error estimator does not allow a decomposition among the atoms
that can be used for atomistic-continuum adaptivity, so we propose and analyze a less accurate
second error estimator that does allow such a decomposition.
Finally, in Section 6 we propose an adaptive atomistic-continuum modeling algorithm and show
that it gives an efficient estimate of the modeling error and a nearly optimal atomistic-continuum
modeling strategy for the computation of defect structure.
2. Quasicontinuum Approximation
The departure point for the QC approximation is the potential energy of the atomistic system.
The potential energy that is utilized fully models the properties of the system. The local minima
of the potential energy model the metastable states of the system, and the potential energy can be
used in Newton’s equations of motion to model the dynamical behavior.
The QC method approximates the potential energy of the atomistic system in two steps. First,
we develop a continuum potential energy that will be used in the adaptively determined continuum
region, and we then show how to reduce the degrees of freedom in the continuum region.
2.1. The Atomistic System. We assume that the atomistic system has 2M atoms with deforma-
tion given by ya = (ya−M+1, . . . y
M) ∈ R
2M . Without loss of generality, we assume that the atoms
are ordered so that their positions satisfy yai < y
i+1. Furthermore, we assume that the atomistic
total potential energy, Ea(ya), can be written as a sum over potential energies associated with each
atom, Eai (y
a), so that
Ea(ya) =
i=−M+1
Eai (y
a). (2.1)
This decomposition can be found for most empirical potentials, including embedded atom potential
energies [6, 7]. For example, if the atomistic total potential energy Ea(ya) is given by
Ea(ya) =
ψ(yaj − y
i ), (2.2)
where ψ(r) is an empirical two-body potential energy, then we can obtain the decomposition (2.1)
by taking
Eai (y
j 6=i
ψ(yaj − y
i ). (2.3)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 3
We note that Eai (y
a) can also contain contributions from external forces, such as for the Frenkel-
Kontorova model described in Section 3, and can thus depend on i.
2.2. The Atomistic-Continuum Energy. For any deformation ya ∈ R2M , we let Li,i+1ya ∈ RZ
denote the linear extrapolation of the atomistic positions yai and y
i+1 given by
(Li,i+1ya)k = (k − i)y
i+1 + (i+ 1− k)y
i for k = −∞, . . . ,∞. (2.4)
The continuum potential energy Eci (y
a) of atom i is obtained from the average of the atomistic
potential energy Eai evaluated at the extrapolations L
i−1,iya and Li,i+1ya by
Eci (y
a) := 1
Eai (L
i−1,i
a) + 1
Eai (L
i,i+1
a), (2.5)
where we note that the domain of Eai has been expanded to the infinite periodic atomistic systems
in the range of Li−1,i and Li,i+1. We assume that Eai is finite for infinite periodic atomistic systems,
which is true for (2.3) when the two-body potential ψ(r) decays fast enough so that
k=1 ψ(kr) is
finite for r 6= 0. At the endpoints of the chain, the extrapolation can be done only to one side, so
we neglect the undefined part and define
Ec−M+1(y
a) := 1
Ea−M+1(L
−M+1,−M+2
a) and EcM (y
a) := 1
EaM (L
M−1,M
a). (2.6)
We then decide for each atom i whether to model its energy atomistically by Eai (y
a) or as a
continuum by Eci (y
a). We thus obtain for the whole chain the atomistic-continuum energy
Eac(ya) :=
i=−M+1
δai E
i=−M+1
δci E
i=−M+1
δai E
a) + 1
i=−M+2
i−1,i
a) + 1
i=−M+1
δci E
i,i+1
(2.7)
where
δai =
1 if atom i is modeled atomistically,
0 if atom i is modeled as continuum,
and δci = 1− δ
i . (2.8)
This approximation allows for a slightly faster evaluation of the energy and its derivatives, especially
if Eai is long-ranged. However, it reveals its full strength only after the quasicontinuum coarsening to
be described next. We note that sometimes atomistic degrees of freedom and energies are referred
to as nonlocal and continuum degrees of freedom and energies are referred to as local [23].
2.3. Repatoms: Reduction of Degrees of Freedom. The quasicontinuum method allows a
reduction of the number of degrees of freedom in the continuum region. To this end, we choose
so-called representative atoms, or more briefly called repatoms. The repatoms are a subset of the
original atoms. The quasicontinuum approximation of the energy is defined completely in terms of
the repatoms.
We choose the repatoms by defining indices ℓj for j = −N + 1, . . . , N where
−M + 1 = ℓ−N+1 < · · · < ℓj < ℓj+1 < · · · < ℓN =M.
The atoms at yai for i = ℓ−N+1, ℓ−N+2, . . . , ℓN are repatoms, and all of the remaining atoms are
non-repatoms. We have that
νj = ℓj+1 − ℓj (2.9)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 4
gives the number of atomistic intervals between the repatoms ℓj and ℓj+1.We require that the chain
not be coarsened in the atomistic regions, which precisely means that δc
= . . . = δc
whenever νj > 1.
Finally, the interactions of the atomistic energy only partially reach into the continuum part if
the atomistic potential has a finite cutoff radius. To allow for an exact calculation of this energy
without atomistic interpolation, we require that these regions are not coarsened as well. As we
will see in the next subsection, the atomistic next-nearest-neighbor interactions from the Frenkel-
Kontorova model studied in this paper reach two atoms into the continuum part. Hence, we require
that νj−2 = νj−1 = νj = νj+1 = 1 whenever δ
= 1. Other potential energies in general require
similar conditions that depend on their cut-off radius.
We denote the position of the j-th repatom by y
j = y
and the vector of all repatoms by
yqc ∈ R2N .
2.4. The Quasicontinuum Energy. Now we define the quasicontinuum energy. To this end, the
missing non-repatoms are implicitly reconstructed. We will see later that this helps to set up the
QC model, but needs not be done for the actual computation.
The reconstruction is done by a linear interpolation between the nearest repatom to the right
and to the left. That is, the vector of all atomistic positions is computed from the vector yqc of
repatom positions by
I : R2N → R2M , (Iyqc)ℓj+m :=
νj −m
j+1, m = 0, . . . , νj . (2.10)
We note that
j = (Iy
qc)ℓj . (2.11)
The underlying idea is that in regions where the lattice spacing of the atoms is nearly constant,
this interpolation is very close to the actual atomistic positions and therefore leads to a good
approximation of the total energy. Only a few repatoms are needed in these regions. This exactly
corresponds to mesh coarsening in classical finite element approximations of continuum models. On
the other hand, in regions where the lattice spacing is non-uniform, such as around a dislocation,
all atoms must be chosen to be repatoms to obtain sufficient accuracy. This guarantees that the
full resolution of the atomistic model in the critical regions is retained and corresponds to a high
refinement in classical finite element continuum models.
We define the QC approximation of the total energy to be
Eqc(yqc) := Eac(Iyqc). (2.12)
Now (2.12) has to be reformulated such that it can be computed efficiently, without the overhead
of the interpolation. Most atomistic potentials are invariant to translations, a property that al-
lows us to simplify (2.12) considerably. For any translationally invariant energy Eai , we have that
Eai (L
i,i+1ya) = φi(y
i+1 − y
i ) and E
i−1,iya) = φi(y
i − y
i−1) for some function φi. If these
functions φi coincide, that is, φi = φj for all i and j, we can write
Eai (L
i−1,i
a) = φ
yai − y
and Eai (L
i,i+1
a) = φ
yai+1 − y
(2.13)
for some function φ : R→ R. Here φ plays the role of a continuum energy density and is given for
the two-body potential (2.2) by
φ(r) =
ψ(kr).
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 5
Equations (2.7), (2.12), and (2.13) lead to
Eqc(yqc) =
i=−M+1
δai E
i (Iy
qc) + 1
i=−M+2
i−1,iIyqc) + 1
i=−M+1
δci E
i,i+1Iyqc)
i=−M+1
δai E
i (Iy
qc) + 1
i=−M+2
δciφ((Iy
qc)i − (Iy
qc)i−1) (2.14)
i=−M+1
δciφ((Iy
qc)i+1 − (Iy
qc)i).
Because Iyqc is the linear interpolation between two repatoms y
j and y
j+1, we have
(Iyqc)i+1 − (Iy
qc)i =
j+1 − y
, i = ℓj, . . . , ℓj+1 − 1. (2.15)
Hence,
Eqc(yqc) =
i=−M+1
δai E
i (Iy
qc) +
j=−N+1
j+1 − y
(2.16)
with weight factors
δcℓj + δ
0 if both y
j and y
j+1 are atomistic,
if exactly one of y
j and y
j+1 is continuum,
νj if both y
and y
j+1 are continuum.
(2.17)
The first sum corresponds to the atomistic region which will be a small region and is thus com-
putationally inexpensive. The second sum only involves at most 2N terms which is a considerable
reduction when N ≪M.
Note that the second term in formula (2.16) coincides with an integral over the energy density φ
as it occurs in finite element discretizations of classical continuum mechanical models. Hence the
apparently unmotivated definitions (2.7) and (2.5) of the continuum energy here result in what is
commonly understood as a continuum energy. The linear interpolation operator I resembles the
Cauchy-Born hypothesis.
3. Frenkel-Kontorova Model
Dislocations are lines in crystals which represent a defect in the lattice structure [15], see Fig-
ure 1. Typically, there is a core of small radius surrounding the dislocation line where the lattice
structure is highly deformed, but the lattice structure is nearly uniform outside the core. A sim-
ple one-dimensional model for a defect such as a dislocation is given by the Frenkel-Kontorova
model [15]. Here, the elastic energy is modeled by harmonic interactions between the atoms in the
one-dimensional chain and the misfit energy of the slip plane is modeled by a periodic potential. A
more accurate model of the same form is given by the Peierls-Nabarro model [11].
3.1. Atomistic Frenkel-Kontorova Model. We study a single defect in the middle of the chain
of 2M atoms. To achieve a symmetric description in terms of bonds, we number the atoms from
−M + 1 to M . The defect is situated between the atoms numbered 0 and 1 (Figure 2).
Recall that the atomistic positions are denoted by ya = (ya−M+1, . . . , y
M ) ∈ R
2M . The total
potential energy for this atomistic system is then a function Ea : R2M → R of the atomistic
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 6
Figure 1: Cross-section through a dislocation in a three-dimensional cubic lattice. The displayed plane
repeats periodically in the three-dimensional crystal. Vertical bonds are shown by lines to emphasize the
topological defect.
(−K−1)a
(K+1)a
= atomistic = continuum
Figure 2: Numbering of the atoms. The dislocation is situated in the middle of the chain between atoms ya
and ya
positions. For the Frenkel-Kontorova model, the energy, Ea = Ea,e + Ea,m, consists of two parts,
namely the part which models the elastic energy of the defect, Ea,e, and the part which models the
misfit energy on the slip plane, Ea,m.
The elastic energy is modeled by Hookean (harmonic) springs between nearest-neighbors (NN)
and next-nearest neighbors (NNN), and the total elastic energy is given by
Ea,e(ya) =
i=−M+1
i+1 − y
i − a0)
i=−M+2
i+1 − y
i−1 − 2a0)
2, (3.1)
where the moduli k1 > 0 and k2 > 0 describe the strength of the elastic interactions, and where
a0 ∈ R denotes the equilibrium distance.
We note that the asymptotic expansion to second order of any nonlinear NN/NNN potential
energy
E(ya) =
i=−M+1
ψ(yai+1 − y
i ) +
i=−M+2
ψ(yai+1 − y
i−1) (3.2)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 7
substrate
modeled
layer
misfit
energy
Figure 3: Frenkel-Kontorova model. The wells depict the misfit energy (3.4).
about aa = [(−M + 1)a0, (−M + 2)a0, · · · , (M − 1)a0,Ma0]
∈ R2M has the form
E(ya) ≈ E(aa) +
ψ′(a0) + 2ψ
′(2a0)
i=−M+1
(yai+1 − y
i − a0)
− ψ′(2a0)(y
M − y
M−1 − a0)− ψ
′(2a0)(y
−M+2 − y
−M+1 − a0)
ψ′′(a0)
i=−M+1
(yai+1 − y
i − a0)
2 + 1
ψ′′(2a0)
i=−M+2
(yai+1 − y
i−1 − 2a0)
(3.3)
We thus see that the elastic energy (3.1) with k1 = ψ
′′(a0) and k2 = ψ
′′(2a0) approximates the
energy (3.2) to second order if ψ′(a0) + 2ψ
′(2a0) = 0 and if we ignore the boundary terms in the
second line of (3.3).
The misfit energy of the slip plane is modeled by a periodic potential (Figure 3). We model this
misfit energy by
Ea,m(ya) =
i=−M+1
yai − a0
, (3.4)
where ⌊x⌋ denotes the largest integer smaller than or equal to x, and where the constant k0 > 0
determines the strength of the misfit energy.
Altogether, the total potential energy of the atomistic system is given by
Ea(ya) = Ea,e(ya) + Ea,m(ya)
i=−M+1
(yai+1 − y
i − a0)
2 + 1
i=−M+2
(yai+1 − y
i−1 − 2a0)
i=−M+1
yai − a0
(3.5)
We restrict ourselves to configurations in which the M leftmost atoms yai for −M + 1 ≤ i ≤ 0 are
situated in the interval
, whereas the M rightmost atoms yai for 1 ≤ i ≤ M
are situated in the interval
. The defect is situated between atoms ya0 and
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 8
ya1 . In this case, the total energy simplifies to
Ea(ya) = 1
i=−M+1
(yai+1 − y
i − a0)
2 + 1
i=−M+2
(yai+1 − y
i−1 − 2a0)
i=−M+1
(yai − (i− 1)a0)
(yai − ia0)
(3.6)
3.2. Quasicontinuum Approximation of the Frenkel-Kontorova Model. We now apply
the quasicontinuum method to the dislocation model described in Section 3.1. The total energy
(3.6) is split up into atom-wise contributions, separately for the elastic interactions and the misfit
interactions:
a) = 1
i − y
i−1 − a0)
2 + 1
i+1 − y
i − a0)
i − y
i−2 − 2a0)
2 + 1
i+2 − y
i − 2a0)
k0 (y
i − (i− 1)a0)
, i = −M + 1, . . . , 0,
k0 (y
i − ia0)
, i = 1, . . . ,M.
(3.7)
To simplify notation, we use the convention that the undefined terms at the endpoints of the chain
are neglected. We thus have that
Ea(ya) = Ea,m(ya) + Ea,e(ya) =
i=−M+1
a) + E
a)] . (3.8)
Since the largest displacement of the atoms is to be expected near the defect, we deem the atoms
−K +1, . . . ,K atomistic and the remaining atoms −M +1, . . . ,−K and K +1, . . . ,M continuum.
Here K < M is some constant whose optimal value will be determined by the algorithm given in
Section 5.
The optimal choice of the repatoms for coarsening is investigated in the second paper of this
series, so we work with a general formulation which holds for any values of ℓj for now. However,
there are two restrictions on the coarsening. Since the atomistic region must not be coarsened and
since we need full refinement in the vicinity of two atoms around the atomistic region due to the
NNN interactions, we have that
ℓj = j, j = −K − 1, . . . ,K + 2. (3.9)
Second, we require that
ℓ−N+1 = −M + 1, ℓ−N+2 = −M + 2, ℓN−1 =M − 1, and ℓN =M (3.10)
to incorporate the boundary conditions later.
The elastic part E
i is translationally invariant, so we perform its QC approximation as described
in the previous section. This leads to the continuum energy density
φe(r) = 1
k1(r − a0)
2 + 1
k2(2r − 2a0)
k12(r − a0)
2 (3.11)
where k12 := k1 + 4k2.
Regarding the misfit part E
i , the above technique cannot be applied since the potential is
not translationally invariant. However, there is a different summation technique to achieve a
computationally efficient formulation which avoids the costly interpolation operator.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 9
To shorten the notation, we let
indicate the sum in which the first term and the last term
are only counted half:
xi :=
i=m+1
xn (3.12)
where m < n and xi ∈ R. It is easy to verify that
2m3 +m
i(m− i) =
m3 −m
(3.13)
for m > 0.
For all pairs (j, j + 1) of continuum repatoms, we now reformulate all terms from (2.7) which
involve the interaction between ℓj and ℓj+1. For j > 0, we get by definition (2.4) of the operator
L, by definition (3.7) of E
i , and by (3.13) that
i=ℓj+1
i−1,iIyqc) + 1
ℓj+1−1
i,i+1Iyqc)
i=ℓj+1
i (Iy
qc) + 1
ℓj+1−1
i (Iy
ℓj+1 − i
i− ℓj
j+1 − ia0
ℓj+1 − i
j − ℓja0) +
i− ℓj
j+1 − ℓj+1a0)
j − ℓja0)
(ℓj+1 − i)
2 + 1
j+1 − ℓj+1a0)
(i− ℓj)
− ℓja0)(y
j+1 − ℓj+1a0)
(ℓj+1 − i)(i − ℓj)
j − ℓja0)
2ν3j + νj
j+1 − ℓj+1a0)
2ν3j + νj
j − ℓja0)(y
j+1 − ℓj+1a0)
ν3j − νj
=: φmj+(y
j , y
j+1).
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 10
For j < 0, we similarly obtain
i=ℓj+1
(Li−1,iIyqc) + 1
ℓj+1−1
(Li,i+1Iyqc)
j − (ℓj − 1)a0)
2ν3j + νj
j+1 − (ℓj+1 − 1)a0)
2ν3j + νj
j − (ℓj − 1)a0)(y
j+1 − (ℓj+1 − 1)a0)
ν3j − νj
=: φmj−(y
j , y
j+1).
Since E
−K (L
−K,−K+1Iyqc) = E
−K (Iy
qc) and E
K+1(L
K,K+1Iyqc) = E
K+1(Iy
qc), the QC approx-
imation of the chain can be given by
Eqc(yqc) =
j=−N+1
j+1 − y
+ φmj−(y
j , y
−K (Iy
qc) +
j=−K+1
Eai (Iy
qc) + 1
K+1(Iy
j+1 − y
+ φmj+(y
j , y
(3.14)
Note that the interpolation Iyqc does not have to be computed here since the relevant terms only
depend on uncoarsened parts of the chain.
Additionally, we will consider the atomistic-continuum approximation
Eac : R2M → R (3.15)
of the atomistic energy without coarsening. It is given exactly like the QC approximation (3.14)
with the only difference being that νj = 1 and ℓj = j everywhere.
4. Primal and Dual Problems
4.1. Problem Setup. We are now ready to set up the problems we will solve. We are interested
in finding the minimum of the energy (3.14) subject to given boundary conditions.
We give the boundary conditions by constraining the deformation of two atoms at each end of the
chain. This guarantees that the potential with next-nearest-neighbor interactions can be directly
applied to all non-boundary atoms without having to neglect interactions. We define the spaces
V a := R2M , V a0 := R
2M−4, V qc := R2N , V
0 := R
2N−4. (4.1)
The spaces V a and V a0 will also be used for the uncoarsened atomistic-continuum potential E
so there is no need to define spaces V ac and V ac0 . We let y
bc ∈ V a denote any vector which has
the desired boundary values ybc−M+1, y
−M+2, y
M−1, and y
M , and we let y
bcq ∈ V qc by any vector
satisfying (recall (3.10))
−N+1 = y
−M+1, y
N−1 = y
−N+2 = y
−M+2, y
= ybcM . (4.2)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 11
For any vector y ∈ V a0 , we denote the extension by zero boundary conditions to be Jy ∈ V
a, so
Jy :=
0 0 yT 0 0
∈ R2M , (4.3)
and similarly we denote the extension by zero boundary conditions of y ∈ V
0 to be J
qcy ∈ V qc.
The spaces of admissible solutions are then given by JV a0 + y
bc ⊂ V a and JqcV
0 + y
bcq ⊂ V qc,
respectively. We note that JT : V a → V a0 is the restriction operator defined by
(JTy)j = yj for j = −M + 3, . . . ,M − 2.
The minima ȳa, ȳac, and ȳqc of the energy functionals Ea, Eac, and Eqc given by (3.6), (3.15),
and (3.14) subject to the above “clamped” boundary conditions are characterized as
a := argmin
y∈JV a
Ea(y) ∈ V a, (4.4)
ac := argmin
y∈JV a
Eac(y) ∈ V a, (4.5)
qc := argmin
y∈JqcV
Eqc(y) ∈ V qc. (4.6)
We note that the minima are uniquely determined because Ea, Eac, and Eqc are strictly convex.
4.2. Matrix Formulation. For the subsequent discussion, it will be convenient to reformulate the
total energies in matrix notation:
Ea(y) = 1
(y − aa)TDaTEaDa(y − aa) + 1
(y − ba)TKa(y − ba), (4.7a)
Eac(y) = 1
(y − aa)TDaTEacDa(y − aa) + 1
(y − ba)TKa(y − ba), (4.7b)
Eqc(y) = 1
(y − aqc)TDqcTEqcDqc(y − aqc) + 1
(y − bqc)TKqc(y − bqc). (4.7c)
The matrices Da ∈ R(2M−1)×2M and Dqc ∈ R(2N−1)×2N compute the distance between two adja-
cent atomistic positions; the matrices Ea ∈ R(2M−1)×(2M−1), Eac ∈ R(2M−1)×(2M−1), and Eqc ∈
(2N−1)×(2N−1) contain the spring constants k1, k2, and k12; and the matricesK
a ∈ R(2M−1)×(2M−1)
and Kqc ∈ R(2N−1)×(2N−1) contain the misfit constant k0. The vectors a
a,ba ∈ R2M and aqc,bqc ∈
2N are constants describing the minimum energy deformations for the elastic energy and mis-
fit energy. The precise and lengthy definitions for all of these matrices and vectors are given in
Appendix A.
If we decompose ȳa = Jya + ybc for ya ∈ V a0 , then the minimization problem is given as
a = argmin
y∈V a
Ea(Jy + ybc)
= argmin
y∈V a
Jy + ybc − aa
DaTEaDa
Jy + ybc − aa
Jy + ybc − ba
Jy + ybc − ba
. (4.8)
We also decompose ȳac = Jyac + ybc and ȳqc = Jqcyqc + ybcq for yac ∈ V a0 and y
qc ∈ V
0 , and
we then formulate similar minimization problems for yac and yqc. Therefore, ya, yac, and yqc are
determined by the linear systems
Maya = fa, (4.9a)
Macyac = fac, (4.9b)
M qcyqc = f qc, (4.9c)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 12
where
Ma := JT (DaTEaDa +Ka)J,
Mac := JT (DaTEacDa +Ka)J,
M qc := JqcT (DqcTEqcDqc +Kqc)Jqc,
a := −JTDaTEaDa(ybc − aa)− JTKa(ybc − ba),
ac := −JTDaTEacDa(ybc − aa)− JTKa(ybc − ba),
qc := −JqcTDqcTEqcDqc(ybcq − aqc)− JqcTKqc(ybcq − bqc).
(4.10)
We note that the matrices Ma, Mac, and M qc are positive definite, so the total energies admit a
single global minimum and no other local minimum.
4.3. Goal-Oriented Error Estimation. To compare the approximate QC model to the original
atomistic model, we have to analyze how much the solution ya of the atomistic model deviates
from the solution yqc of the QC model. This deviation, which can be viewed as an approximation
error, can be measured in different ways, for example as ‖ya−JT IJqcyqc‖ for some norm ‖·‖. Here
we follow a different approach, namely we measure the error of a quantity of interest denoted by
Q(y) for some function Q : R2M−4 → R. Hence, we intend to estimate
Q(ya)−Q(JT IJqcyqc). (4.11)
We will assume for simplicity that Q is linear and thus has a representation Q(y) = qTy for some
vector q ∈ V a0 .
For our application, a natural quantity of interest is the size of the dislocation, that is, the
distance between the two atoms y0 and y1 to the left and right of the dislocation. This gives us
Q(y) = qTy = y1 − y0 with q = [0, . . . , 0,−1, 1, 0, . . . , 0]
T . (4.12)
Two different sources of error arise during the QC approximation, namely the localization of the
potential energy, that is, the passage from the atomistic to the continuum formulation on the one
hand, and the coarsening in the continuum region by the restriction to the repatoms on the other
hand. We denote these two errors by
e := ya − yac and eacqc := yac − JT IJqcyqc. (4.13)
It makes sense to study these sources independently. Employing the linearity of Q, we have that
|Q(ya)−Q(JT IJqcyqc)| = |Q(e) +Q(eacqc)| ≤ |Q(e)|+ |Q(eacqc)|. (4.14)
The error term |Q(e)| will be studied in Section 5, and the error term |Q(eacqc)| will be studied in
the second part of this paper series.
4.4. Dual Problems. To facilitate the goal-oriented error analysis, we introduce the dual problems
Maga = q, (4.15a)
Macgac = q, (4.15b)
M qcgqc = JqcT ITJq, (4.15c)
for ga,gac ∈ R2M−4, and gqc ∈ R2N−4. We note that the dual problems differ from the primal
problems only by the right hand side since the matrices Ma, Mac, and M qc are symmetric.
The solutions ga, gac and gqc can be viewed as influence functions: They describe how the error
at a specific point in the domain influences the error measured in terms of the goal function.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 13
Analogously to the primal errors (4.13), we define the dual errors
ê := ga − gac and êacqc := gac − JT IJqcgqc. (4.16)
In addition, we will need the primal and dual residuals
Ra(y) :=Ma (ya − y) = fa −May,
Rac(y) :=Mac (yac − y) = fac −Macy,
R̂a(g) :=Ma (ga − g) = q−Mag,
R̂ac(g) :=Mac (gac − g) = q−Macg. (4.17)
5. Error Estimation for Atomistic vs. Continuum Modeling
In this section, we estimate the error |Q(e)| arising from the approximation of an atomistic model
by a continuum model. We consider yac and gac to be computable, although in practice we can
only compute the coarsened approximations yqc and gqc.
To this end, we adapt a technique introduced in [16] and [19] to estimate the modeling error for
an elasticity model with rapidly oscillating coefficients and its homogenized version. We generalize
this technique such that it allows for different right hand sides fa and fac of the primal problem
(4.9) instead of a common right hand side as it is used in the above-mentioned works.
We have
Q(ya)−Q(yac) = qTe = gaTMae = (gacT + ê)Mae
= gacTRa(yac) + êTMae. (5.1)
The term gacTRa(yac) can be computed, whereas êTMae cannot because both e and ê are nu-
merically unknown. Instead, we estimate êTMae from above and from below by quantities that
actually can be computed.
We will give two different error estimators η1 and η2. Before, we need to derive some auxiliary
estimates to facilitate their development and analysis.
5.1. Auxiliary Estimates. We reformulate the difference ya − yac of the respective solutions in
terms of a difference of the energy matrices. To this end, we define the perturbation matrix
P := I − (Ea)−1Eac (5.2)
where I denotes the identity matrix. Note that EaP = E
a − Eac.
Lemma 5.1. For any α, β ∈ R, we have that
Ma(αe+ βê) = −JTDaTEaPDa
α(Jyac + ybc − aa) + βJgac
. (5.3)
Proof. We conclude from (4.9) that
Mae =Maya −Macyac + (Mac −Ma)yac
= fa − fac + (Mac −Ma)yac, (5.4)
and similarly since Maga =Macgac = q that
Maê =Ma(ga − gac) = (Mac −Ma)gac. (5.5)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 14
Thus, it follows from (4.10) and (5.2) that
Ma(αe+ βê) = α [(Mac −Ma)yac + fa − fac] + β(Mac −Ma)gac
= JTDaT (Eac − Ea)Da
α(Jyac + ybc − aa) + βJgac
= −JTDaTEaPDa
α(Jyac + ybc − aa) + βJgac
(5.6)
We note that the Ka-related terms cancel here, because they coincide for the atomistic model and
the continuum model. �
Lemma 5.2. We have that
‖αe+ βê‖Ma ≤
α(Jyac + ybc − aa) + βJgac
. (5.7)
We note that the right hand side is numerically computable.
Proof. To shorten the notation, we abbreviate z = α(Jyac +ybc− aa)+ βJgac. By Lemma 5.1, we
‖αe + βê‖Ma = sup
v∈V a
vTMa(αe+ βê)
‖v‖Ma
= sup
v∈V a
−vTJTDaTEaPDaz
‖v‖Ma
≤ sup
v∈V a
‖DaJv‖Ea‖PD
az‖Ea
‖DaJv‖Ea
= ‖PDaz‖Ea .
(5.8)
Here we have used that ‖DaJv‖Ea ≤ ‖v‖Ma because the matrixK
a in (4.10) is positive definite. �
5.2. First Error Estimator. We are now ready to derive the first error estimator, η1. By the
parallelogram identity, we have for all σ 6= 0 that
TMae = (σ−1êT )Ma(σe)
‖σe+ σ−1ê‖2Ma −
‖σe− σ−1ê‖2Ma .
(5.9)
In the following, we will determine computable constants η+low, η
low, η
upp and η
upp such that
η+low ≤ ‖σe+ σ
ê‖Ma ≤ η
η−low ≤ ‖σe− σ
ê‖Ma ≤ η
(5.10)
From Lemma 5.2, we immediately get the upper estimates η+upp and η
η+upp :=
σ(Jyac + ybc − aa) + σ−1Jgac
η−upp :=
σ(Jyac + ybc − aa)− σ−1Jgac
(5.11)
We note that η+low, η
low, η
upp and η
upp will depend on σ, but the estimates (5.10) will hold for any
σ 6= 0. We will now choose σ in such a way that the estimates are as sharp as possible, that is,
such that η+upp and η
upp are smallest.
Lemma 5.3. Both η+upp and η
upp given by (5.11) attain their minima for
σ̄ :=
‖PDaJgac‖Ea
‖PDa(Jyac + ybc − aa)‖Ea
. (5.12)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 15
Proof. We have that
(η±upp)
2 = σ2
∥PDa(Jyac + ybc − aa)
± 2gacTJTDaTP TEaPDa(Jyac + ybc − aa)
+ σ−2
∥PDaJgac
. (5.13)
Setting the first derivative of the mapping σ 7→ (η±upp)
2 to zero, we obtain the condition
∥PDa(Jyac + ybc − aa)
− 2σ̄−3
∥PDaJgac
= 0 (5.14)
for critical points of (η±upp)
2. This equation has the unique positive solution (5.12). Because
lim|σ|→∞ η
upp = limσ→0 η
upp =∞, this point corresponds to a minimum. Hence the quantities η
attain their minima at σ = σ̄. �
Regarding the lower bounds η+low and η
low, we have
‖σ̄e± σ̄−1ê‖Ma = sup
v∈V a
vTMa(σ̄e± σ̄−1ê)
‖v‖Ma
= sup
v∈V a
vT (σ̄Ra(yac)± σ̄−1R̂a(gac))
‖v‖Ma
vT0 (σ̄R
a(yac)± σ̄−1R̂a(gac))
‖v0‖Ma
(5.15)
for any vector v0 ∈ V
0 \ {0}. Numerically, we have the two vectors y
ac and gac at our disposal,
hence it makes sense to take a linear combination v0 = y
ac + θ±gac. Here we follow the strategy
of [19] and choose θ± as the critical points of η±low.
Lemma 5.4. Let
± = σ̄Ra(yac)± σ̄−1R̂a(gac). (5.16)
Then the lower bounds
(yac + θ±gac)T r±
‖yac + θ±gac‖Ma
(5.17)
have a unique critical point for
θ̄± :=
r±Tyac gacTMayac − r±Tgac ‖yac‖2Ma
r±Tgac gacTMayac − r±Tyac ‖gac‖2
. (5.18)
Proof. We have
r±Tgac‖yac + θ±gac‖Ma − r
±T (yac + θ±gac)
(yac+θ±gac)TMagac
‖yac+θ±gac‖Ma
‖yac + θ±gac‖2
. (5.19)
Setting this expression to zero and solving for θ± leads to the above condition. �
However, let us note that this critical point is not necessarily a maximum of η±low, which would
be optimal for bound (5.10). Depending on the actual vectors yac and gac, it can be shown that
this critical point could be a minimum.
Now we have all necessary ingredients to construct the error estimator η1. From (5.1) and (5.9),
we get the computable estimate
gacTRa(yac) + 1
)2 − 1
(η−upp)
2 ≤ Q(ya)−Q(yac) ≤ gacTRa(yac) + 1
(η+upp)
2 − 1
(5.20)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 16
At first sight, this looks like we could get an estimate for |Q(ya)−Q(yac)| from both above and
below. However, this is only true if both the left hand side and the right hand side have the same
sign, which in general does not hold. But we get the following estimate.
Theorem 5.1. We have that
|Q(ya)−Q(yac)| ≤ η1, (5.21)
where the computable error estimator is defined as
η1 := max
∣gacTRa(yac) + 1
)2 − 1
(η−upp)
∣gacTRa(yac) + 1
(η+upp)
2 − 1
. (5.22)
We note that the computation of the η±upp terms involves the solution of a linear system with
matrix Ea as its inverse appears in the operator P . The matrix Ea is not diagonal, but has
condition number O(1). So this is negligible compared to what would be necessary to solve the
original atomistic problem which includes the operator DaTEaDa with condition number O(M2).
5.3. Second Error Estimator. There is no reasonable way to decompose the error estimator η1
into a sum of element-wise or atom-wise contributions due to the η±
terms. Therefore, we derive
another error estimator η2 which allows for such a decomposition, at the price of a less accurate
estimate than η1.
Theorem 5.2. We have that
|Q(ya)−Q(yac)| ≤ η2 ≤
i=−M+3
ηat2,i +
i=−M+1
ηel2,i (5.23)
where the computable global error estimator η2 and the computable local error estimators, η
2,i and
ηel2,i, associated with atoms and elements, respectively, are defined as
η2 :=
∣gacTRa(yac)
∣+ ‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea,
ηat2,i := |g
a(yac)i| , i = −M + 3, . . . ,M − 2,
ηel2,i :=
PDa(Jyac + ybc − aa)
(Ea − Eac)Da(Jyac + ybc − aa)
∣(PDaJgac)i
(Ea − Eac)DaJgac
∣ , i = −M + 1, . . . ,M − 1.
(5.24)
Proof. From (5.1) and Lemma 5.2, we conclude that
|Q(ya)−Q(yac)| ≤
acTRa(yac)
acTRa(yac)
∣+ ‖ê‖Ma‖e‖Ma
acTRa(yac)
∣+ ‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea
= η2,
(5.25)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 17
which gives us the global estimate. For the decomposition into local contributions, we further
estimate
acTRa(yac)
∣+ ‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea
acTRa(yac)
‖PDa(Jyac + ybc − aa)‖2Ea +
‖PDaJgac‖2Ea
i=−M+3
|gaci R
a(yac)i|+
i=−M+1
∣(PDaJgac)i
(Ea − Eac)DaJgac
i=−M+1
PDa(Jyac + ybc − aa)
(Ea − Eac)Da(Jyac + ybc − aa)
i=−M+3
ηat2,i +
i=−M+1
ηel2,i,
(5.26)
which completes the proof. �
Let us remark that instead of the first inequality in (5.26), one can get an apparently better
estimate
‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea
γ‖PDa(Jyac + ybc − aa)‖2Ea +
γ−1‖PDaJgac‖2Ea
(5.27)
by introducing the additional weight factor
‖PDaJgac‖Ea
‖PDa(Jyac + ybc − aa)‖Ea
, (5.28)
and then decomposing the resulting terms similar to the above. However, our numerical results
showed that this modification does not significantly improve the decomposed error estimator for
the application considered here.
6. Numerics
In the preceding sections, we constructed the error estimators η1 and η2. We will now give an
algorithm for adaptive atomistic-continuum modeling based on these error estimators. Then we
will present and discuss some numerical results.
6.1. Algorithm. The error estimator η1 should give a better estimate of the error than η2, because
η2 involves the inequality
∣êTMae
∣ ≤ ‖ê‖Ma‖e‖Ma in (5.25) in contrast to the parallelogram identity
for η1. However, η2 can be decomposed into atom-wise and element-wise contributions η
2,i and
ηel2,i, whereas the η
terms in η1 do not admit a reasonable decomposition that can be used for
atomistic-continuum adaptivity.
We make use of this by employing the sharper estimate η1 to determine whether a given global
error tolerance τgl for the error in an adaptive algorithm has already been achieved or not. If
not, we use the decomposed estimates ηat2,i and η
2,i to determine where the more precise atomistic
modeling is needed. This leads us to the following algorithm:
(1) Choose τgl. Model all atoms as a continuum. Set τat ← τgl.
(2) Solve primal problem (4.9b) for yac and dual problem (4.15b) for gac.
(3) Compute error estimator η1 from (5.22).
(4) If η1 ≤ τgl, then stop.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 18
M iteration K τat η1
100 1 0 1.000000e-10 3.899207e-02
2 28 1.000000e-11 5.915080e-10
3 32 1.000000e-12 4.878532e-11
1000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.915100e-10
3 32 1.000000e-12 4.878548e-11
10000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.915100e-10
3 32 1.000000e-12 4.878548e-11
100000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.915099e-10
3 32 1.000000e-12 4.878540e-11
1000000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.914422e-10
3 32 1.000000e-12 4.871775e-11
Table 1: Convergence of the algorithm for τgl = 10
−10 and different values of M .
(5) Compute local error estimators ηat2,i and η
2,i from (5.24).
(6) Set τat ←
(7) Make all atoms i atomistic for which
ηtot2,i := η
2,i +
ηel2,i−1 + η
≥ τat. (6.1)
(8) Go to (2).
Here τdiv > 1 is a constant factor which describes how fast the atom-wise tolerance τat should
decrease during adaption. Our experience has been that τdiv = 10 is a reasonable choice.
The crucial adaption step is (7). The adaption criterion (6.1) deems all atoms to be modeled
atomistically if the associated error from the decomposition of η2 exceeds the atomistic error toler-
ance τat. Here the element-wise errors η
2,i are distributed equally to the two adjacent atoms i and
i+ 1.
For the dislocation at the center of the chain and the chosen goal function, we expect that the
atomistic repatoms always form a symmetric interval around the center. We have used the above
adaptive atomistic-continuum algorithm to approximate our Frenkel-Kontorova model and have
always found that the atomistic region is the set of atoms −K + 1, . . . ,K for some K depending
on M and τgl. Thus, the modeling approach given in Section 3 of restricting to an atomistic region
consisting of atoms −K + 1, . . . ,K for some K rather than considering a more general atomistic
region is justified a posteriori.
6.2. Numerical Results. The algorithm has been implemented as described above. The boundary
conditions were chosen as
ybc−M+1 = −M, y
−M+2 = −M + 1, y
M−1 =M − 1, y
M =M. (6.2)
The elastic constants are k0 = 1 and k1 = k2 = 2.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 19
−500 −400 −300 −200 −100 0 100 200 300 400 500
−500 −400 −300 −200 −100 0 100 200 300 400 500
Figure 4: The error estimators ηtot
2,i (left) and η
2,i, η
2,i (right) for M = 500, K = 20.
Table 1 shows how the algorithm given above performs. After 3 iterations, the desired accuracy
τgl = 10
−10 is achieved. Moreover, we can see from the table that the number of iterations are
independent of M , that means the algorithm behaves robustly with respect to the problem size M .
Figure 4 (left) shows the decomposition of the error estimator η2 for a typical setting M = 500,
K = 20. One can clearly see that the error in the atomistic region is small, whereas the error is
large in the continuum regions that border the atomistic region. It then decreases exponentially
towards the endpoints.
The error in both the atomistic region around the center and the continuum regions far away
from the center are in the range of the (relative) machine precision εmach, which accounts for the
fluctuations in these regions. The error can be considered to be numerically zero in these regions.
In the continuum regions, we observe an error of magnitude O(ε2
), whereas in the continuum
region we have O(εmach), which leads to the different magnitudes of the fluctuations.
Figure 4 (right) shows the element-wise contributions ηel2,i and the atom-wise contributions η
of the decomposed error estimator ηtot2,i = η
2,i +
ηel2,i−1 + η
. The atomistic part ηat2,i, which
corresponds to the gacTRa(yac) term, is dominant in the sense that it is about ten times larger
than ηel2,i, which comes from the estimate for the perturbation term ê
TMae. The fluctuations
due to the limited machine precion in the atomistic region come from ηat2,i, whereas those in the
continuum region away from the defect stem from ηel2,i. Let us note that in other applications of
duality-based error estimation, the first term might not always be the dominant term. For example,
in mesh refinement for classical linear finite elements, the first term even vanishes due to Galerkin
orthogonality.
Table 2 and Figure 5, which display the same data, show the efficiency of the error estimators,
η1/|Q(y
a − yac)| and η2/|Q(y
a − yac)|, for M = 1000. For comparison, the actual error is given as
well. For the relatively small 1D problem, the actual error can be easily computed, whereas in real
world applications it is of course not available. One can clearly see that η1 gives a better estimate
than η2, which numerically confirms our conjecture that η1 is a better estimator than η2. We see
that η1 overestimates the actual error by a factor of 1.4, while η2 is in a still acceptable range of
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 20
K |Q(ya − yac)| η1 η1/|Q(y
a − yac)| η2 η2/|Q(y
a − yac)|
0 3.627633e-02 3.899208e-02 1.074863 3.999783e-02 1.102588
2 3.375762e-02 3.872272e-02 1.147081 5.101700e-02 1.511274
4 3.468605e-03 4.343595e-03 1.252260 5.422007e-03 1.563166
6 5.418585e-04 7.156249e-04 1.320686 9.187940e-04 1.695635
8 1.227067e-04 1.675383e-04 1.365356 2.193196e-04 1.787348
10 3.287188e-05 4.540984e-05 1.381419 5.984186e-05 1.820457
15 1.416914e-06 1.966114e-06 1.387603 2.597488e-06 1.833201
20 6.267636e-08 8.695824e-08 1.387417 1.148736e-07 1.832805
25 2.770161e-09 3.843388e-09 1.387424 5.077204e-09 1.832819
30 1.224369e-10 1.698739e-10 1.387440 2.244073e-10 1.832840
35 5.410783e-12 7.508365e-12 1.387667 9.918687e-12 1.833133
40 2.379208e-13 3.318024e-13 1.394592 4.383361e-13 1.842362
45 8.992806e-15 1.430733e-14 1.590975 1.901601e-14 2.114580
50 7.771561e-16 4.120094e-16 0.530150 6.201285e-16 0.797946
Table 2: Efficiency of the error estimators, η1/|Q(y
a − yac)| and η2/|Q(y
a − yac)|, for M = 1000. For
K = 45 and K = 50 the results become inaccurate due to limited machine precision.
0 5 10 15 20 25 30 35 40 45 50
|Q(ya−yac)|
Figure 5: Efficiency of the error estimators for M = 1000.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 21
τgl optimal K K by η1 K by η2
1e-02 3 3 3
1e-03 5 5 5
1e-04 9 9 10
1e-05 12 13 13
1e-06 16 17 17
1e-07 20 20 21
1e-08 23 24 24
1e-09 27 28 28
1e-10 31 31 32
1e-11 35 35 35
1e-12 38 39 39
1e-13 42 42 43
1e-14 45 46 47
Table 3: Efficiency of the error estimators for M = 1000.
2 times the actual error. Moreover, we can see from Table 2 and Figure 5 that the error decreases
exponentially with K.
Finally, we compare the optimal (smallest) value of K which is needed to achieve a given accuracy
τgl with the values for K determined by the error estimators η1 and η2, again taking into account
the precise error which is available for the model problem. We see from Table 3 that even η2 only
overestimates K by at most 2 atoms. Thus, we get an efficient estimate of the required atomistic
region for both error estimators.
Appendix A. Matrix Definitions
We describe the matrices from Section 4.2. The matrix
. . .
. . .
∈ R(2M−1)×2M (A.1)
transforms atomistic positions to distances between adjacent atoms. Similarly,
Dqc =
−ν−1−N+1 ν
−ν−1−N+2 ν
. . .
. . .
N−1 ν
∈ R(2N−1)×2N (A.2)
transforms repatom positions from a coarsened chain to normalized distances between adjacent
repatoms. The matrices
(Ea)ij =
k1 + k2 i = j ∈ {−M + 1,M − 1}
k1 + 2k2 i = j ∈ {−M + 2, . . . ,M − 2}
k2 |j − i| = 1
0 otherwise,
(A.3)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 22
(Eac)ij =
k12(δ
i + δ
i+1) +
i + δ
i+1) +
i−1 + δ
i + δ
i+1 + δ
i+2) i = j
i + δ
i+2) j = i+ 1
i−1 + δ
i+1) j = i− 1
0 otherwise,
(A.4)
(Eqc)ij =
ωik12 +
i + δ
i+1) +
i−1 + δ
i + δ
i+1 + δ
i+2) i = j
i + δ
i+2) j = i+ 1
i−1 + δ
i+1) j = i− 1
0 otherwise,
(A.5)
for i, j = −M + 1, . . . ,M − 1 and i, j = −N + 1, . . . , N − 1, respectively, describe the spring
interactions in terms of the distances between atoms or repatoms. Accordingly, the matrices
. . .
∈ R2M×2M (A.6)
(Kqc)ij =
(2νi−1 + ν
i−1) + (2νi + ν
i = j ∈ {−N + 2, . . . , N − 1}
k0(2ν−N+1 + ν
−N+1) i = j = −N + 1
k0(2νN−1 + ν
N−1) i = j = N
k0(νi − ν
i ) j = i+ 1
k0(νj − ν
j ) j = i− 1
0 otherwise,
(A.7)
for i, j = −N + 1, . . . , N describe the misfit interactions for the original atomistic system and the
QC approximation. Finally, the constant vectors
(−M + 1)a0 (−M + 2)a0 · · · (M − 1)a0 Ma0
∈ R2M , (A.8a)
ℓ−N+1a0 ℓ−N+2a0 · · · ℓN−1a0 ℓNa0
∈ R2N , (A.8b)
−Ma0 (−M + 1)a0 · · · −a0 a0 · · · (M − 1)a0 Ma0
∈ R2M , (A.8c)
(ℓ−N+1 − 1)a0 (ℓ−N+2 − 1)a0 · · · (ℓ0 − 1)a0 ℓ1a0 · · · ℓN−1a0 ℓNa0
∈ R2M ,
(A.8d)
fix the equilibrium positions for the spring interactions and the misfit energy, respectively.
References
[1] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York,
2000.
[2] M. Arndt and M. Griebel, Derivation of higher order gradient continuum models from atomistic models for
crystalline solids, Multiscale Model. Simul., 4 (2005), pp. 531–562.
[3] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Lectures in
Mathematics, ETH Zürich, Birkhäuser, Basel, 2003.
[4] X. Blanc, C. Le Bris, and F. Legoll, Analysis of a prototypical multiscale method coupling atomistic and
continuum mechanics, Math. Model. Numer. Anal., 39 (2005), pp. 797–826.
[5] P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press,
2000.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 23
[6] M. S. Daw and M. I. Baskes, Semiempirical, quantum mechanical calculation of hydrogen embrittlement in
metals, Phys. Rev. Lett., 50 (1983), pp. 1285–1288.
[7] , Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals,
Phys. Rev. B, 29 (1984), pp. 6443–6453.
[8] M. Dobson and M. Luskin, Analysis of a force-based quasicontinuum approximation, 2006,
arXiv:math.NA/0611543.
[9] W. E, J. Lu, and J. Z. Yang, Uniform accuracy of the quasicontinuum method, Phys. Rev. B, 74 (2006),
p. 214115.
[10] W. E and P. Ming, Analysis of multiscale methods, J. Comput. Math., 22 (2004), pp. 210–219.
[11] E. Kaxiras, Atomic and Electronic Structure of Solids, Cambridge University Press, 2003.
[12] J. Knap and M. Ortiz, An analysis of the quasicontinuum method, J. Mech. Phys. Solids, 49 (2001), pp. 1899–
1923.
[13] P. Lin, Theoretical and numerical analysis for the quasi-continuum approximation of a material particle model,
Math. Comput., 72 (2003), pp. 657–675.
[14] , Convergence analysis of a quasi-continuum approximation for a two-dimensional material without defects,
SIAM J. Numer. Anal., 45 (2007), pp. 313–332.
[15] M. Marder, Condensed Matter Physics, John Wiley & Sons, 2000.
[16] J. T. Oden and S. Prudhomme, Estimation of modeling error in computational mechanics, J. Comput. Phys.,
182 (2002), pp. 496–515.
[17] J. T. Oden, S. Prudhomme, and P. Bauman, Error control for molecular statics problems, Int. J. Multiscale
Comput. Eng., 4 (2006), pp. 647–662.
[18] J. T. Oden, S. Prudhomme, A. Romkes, and P. Bauman, Multiscale modeling of physical phenomena:
Adaptive control of models, SIAM J. Sci. Comput., 28 (2006), pp. 2359–2389.
[19] J. T. Oden and K. S. Vemaganti, Estimation of local modeling error and goal-oriented adaptive modeling
of heterogeneous materials: Part I: Error estimates and adaptive algorithms, J. Comput. Phys., 164 (2000),
pp. 22–47.
[20] C. Ortner and E. Süli, A-posteriori analysis and adaptive algorithms for the quasicontinuum method in one
dimension, Research Report NA-06/13, Oxford University Computing Laboratory, 2006.
[21] , A-priori analysis of the quasicontinuum method in one dimension, Research Report NA-06/12, Oxford
University Computing Laboratory, 2006.
[22] E. B. Tadmor, R. Miller, R. Phillips, and M. Ortiz, Nanoindentation and incipient plasticity, J. Mater.
Res., 14 (1999), pp. 2233–2250.
[23] E. B. Tadmor, M. Ortiz, and R. Phillips, Quasicontinuum analysis of defects in solids, Philos. Mag. A, 73
(1996), pp. 1529–1563.
[24] E. B. Tadmor, R. Phillips, and M. Ortiz, Mixed atomistic and continuum models of deformation in solids,
Langmuir, 12 (1996), pp. 4529–4534.
	1. Introduction
	2. Quasicontinuum Approximation
	2.1. The Atomistic System
	2.2. The Atomistic-Continuum Energy
	2.3. Repatoms: Reduction of Degrees of Freedom
	2.4. The Quasicontinuum Energy
	3. Frenkel-Kontorova Model
	3.1. Atomistic Frenkel-Kontorova Model
	3.2. Quasicontinuum Approximation of the Frenkel-Kontorova Model
	4. Primal and Dual Problems
	4.1. Problem Setup
	4.2. Matrix Formulation
	4.3. Goal-Oriented Error Estimation
	4.4. Dual Problems
	5. Error Estimation for Atomistic vs. Continuum Modeling
	5.1. Auxiliary Estimates
	5.2. First Error Estimator
	5.3. Second Error Estimator
	6. Numerics
	6.1. Algorithm
	6.2. Numerical Results
	Appendix A. Matrix Definitions
	References

  We propose and analyze a goal-oriented a posteriori error estimator for the
atomistic-continuum modeling error in the quasicontinuum method. Based on this
error estimator, we develop an algorithm which adaptively determines the
atomistic and continuum regions to compute a quantity of interest to within a
given tolerance. We apply the algorithm to the computation of the structure of
a crystallographic defect described by a Frenkel-Kontorova model and present
the results of numerical experiments. The numerical results show that our
method gives an efficient estimate of the error and a nearly optimal
atomistic-continuum modeling strategy.

Introduction
The quasicontinuum (QC) method [22, 23, 24] has been successfully used to efficiently cou-
ple atomistic and continuum models for crystalline solids and offers the possibility of computing
mesoscale or macroscale properties by a nearly minimal number of degrees of freedom. Accurate
modeling requires that an atomistic model be used in regions with highly non-uniform deformations
such as around dislocations, whereas a continuum model can be used in regions with nearly uniform
deformations to reduce the number of degrees of freedom.
It is usually not known a priori which regions of some specimen undergo uniform deformations
and which do not, so a posteriori error estimation is important for the design of efficient numerical
approximations by the quasicontinuum method. Since the purpose of a computation is often to
obtain the value of a (usually local) quantity of interest to a desired error tolerance rather than to
obtain a solution to a desired error tolerance for a global norm, there has been great interest in
the development of goal-oriented error estimators for many problems. They are based on duality
techniques and have been developed and used to adaptively refine finite element approximations of
continuum problems [1, 3] and to study and control modeling error [19].
In this paper, we extend this approach to develop an a posteriori error estimator for the qua-
sicontinuum method which quantifies the atomistic-continuum modeling error for a goal function
and allows for an adaptive decision about which regions can be accurately modeled as a continuum
and which regions need to be modeled atomistically. Methods to determine the optimal mesh size
within the continuum region will be studied in a forthcoming paper.
Crystallographic defects [5] provide a challenge to validate atomistic-continuum error estimators
and adaptivity. No such error estimators and adaptive methods currently exist for fully three-
dimensional crystals. As a step in this direction, we develop a rigorous theory for a simple one-
dimensional atomistic model for a defect that is a modification of the Frenkel-Kontorova model [15].
We add next-nearest-neighbor harmonic interactions between the atoms to the nearest-neighbor
harmonic interactions between the atoms in the classical Frenkel-Kontorova model.
A priori analyses for various quasicontinuum approximations have been given in [4, 8, 9, 10,
12, 13, 14, 21]. An a posteriori analysis for a slightly different one-dimensional quasicontinuum
http://arxiv.org/abs/0704.1924v1
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 2
approximation is given in [20]. The development and application of a goal-oriented error estimator
for mesh coarsening in a two-dimensional quasicontinuum method is reported in [17, 18].
Let us mention that the continuum model used in the QC method, which coincides with the
model obtained by the classical thermodynamic limit, is by far not the only reasonable continuum
model to use. A method to derive continuum models which approximate atomistic models up to
an arbitrarily high order has been proposed in [2].
The paper is organized as follows. In Section 2, we give a general formulation of the one-
dimensional quasicontinuum approximation [23] that includes not only two-body and three-body
potentials, but also many body potentials such as the embedded atom potential [6, 7]. In Section 3,
we describe our extension of the Frenkel-Kontorova model and its quasicontinuum approximation.
In Section 4, we introduce the primal and dual problems for our model and formulate our approach
to goal-oriented error estimation.
Next, in Section 5 we extend the approach in [16] to develop an error estimator for atomistic-
continuum modeling. This first error estimator does not allow a decomposition among the atoms
that can be used for atomistic-continuum adaptivity, so we propose and analyze a less accurate
second error estimator that does allow such a decomposition.
Finally, in Section 6 we propose an adaptive atomistic-continuum modeling algorithm and show
that it gives an efficient estimate of the modeling error and a nearly optimal atomistic-continuum
modeling strategy for the computation of defect structure.
2. Quasicontinuum Approximation
The departure point for the QC approximation is the potential energy of the atomistic system.
The potential energy that is utilized fully models the properties of the system. The local minima
of the potential energy model the metastable states of the system, and the potential energy can be
used in Newton’s equations of motion to model the dynamical behavior.
The QC method approximates the potential energy of the atomistic system in two steps. First,
we develop a continuum potential energy that will be used in the adaptively determined continuum
region, and we then show how to reduce the degrees of freedom in the continuum region.
2.1. The Atomistic System. We assume that the atomistic system has 2M atoms with deforma-
tion given by ya = (ya−M+1, . . . y
M) ∈ R
2M . Without loss of generality, we assume that the atoms
are ordered so that their positions satisfy yai < y
i+1. Furthermore, we assume that the atomistic
total potential energy, Ea(ya), can be written as a sum over potential energies associated with each
atom, Eai (y
a), so that
Ea(ya) =
i=−M+1
Eai (y
a). (2.1)
This decomposition can be found for most empirical potentials, including embedded atom potential
energies [6, 7]. For example, if the atomistic total potential energy Ea(ya) is given by
Ea(ya) =
ψ(yaj − y
i ), (2.2)
where ψ(r) is an empirical two-body potential energy, then we can obtain the decomposition (2.1)
by taking
Eai (y
j 6=i
ψ(yaj − y
i ). (2.3)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 3
We note that Eai (y
a) can also contain contributions from external forces, such as for the Frenkel-
Kontorova model described in Section 3, and can thus depend on i.
2.2. The Atomistic-Continuum Energy. For any deformation ya ∈ R2M , we let Li,i+1ya ∈ RZ
denote the linear extrapolation of the atomistic positions yai and y
i+1 given by
(Li,i+1ya)k = (k − i)y
i+1 + (i+ 1− k)y
i for k = −∞, . . . ,∞. (2.4)
The continuum potential energy Eci (y
a) of atom i is obtained from the average of the atomistic
potential energy Eai evaluated at the extrapolations L
i−1,iya and Li,i+1ya by
Eci (y
a) := 1
Eai (L
i−1,i
a) + 1
Eai (L
i,i+1
a), (2.5)
where we note that the domain of Eai has been expanded to the infinite periodic atomistic systems
in the range of Li−1,i and Li,i+1. We assume that Eai is finite for infinite periodic atomistic systems,
which is true for (2.3) when the two-body potential ψ(r) decays fast enough so that
k=1 ψ(kr) is
finite for r 6= 0. At the endpoints of the chain, the extrapolation can be done only to one side, so
we neglect the undefined part and define
Ec−M+1(y
a) := 1
Ea−M+1(L
−M+1,−M+2
a) and EcM (y
a) := 1
EaM (L
M−1,M
a). (2.6)
We then decide for each atom i whether to model its energy atomistically by Eai (y
a) or as a
continuum by Eci (y
a). We thus obtain for the whole chain the atomistic-continuum energy
Eac(ya) :=
i=−M+1
δai E
i=−M+1
δci E
i=−M+1
δai E
a) + 1
i=−M+2
i−1,i
a) + 1
i=−M+1
δci E
i,i+1
(2.7)
where
δai =
1 if atom i is modeled atomistically,
0 if atom i is modeled as continuum,
and δci = 1− δ
i . (2.8)
This approximation allows for a slightly faster evaluation of the energy and its derivatives, especially
if Eai is long-ranged. However, it reveals its full strength only after the quasicontinuum coarsening to
be described next. We note that sometimes atomistic degrees of freedom and energies are referred
to as nonlocal and continuum degrees of freedom and energies are referred to as local [23].
2.3. Repatoms: Reduction of Degrees of Freedom. The quasicontinuum method allows a
reduction of the number of degrees of freedom in the continuum region. To this end, we choose
so-called representative atoms, or more briefly called repatoms. The repatoms are a subset of the
original atoms. The quasicontinuum approximation of the energy is defined completely in terms of
the repatoms.
We choose the repatoms by defining indices ℓj for j = −N + 1, . . . , N where
−M + 1 = ℓ−N+1 < · · · < ℓj < ℓj+1 < · · · < ℓN =M.
The atoms at yai for i = ℓ−N+1, ℓ−N+2, . . . , ℓN are repatoms, and all of the remaining atoms are
non-repatoms. We have that
νj = ℓj+1 − ℓj (2.9)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 4
gives the number of atomistic intervals between the repatoms ℓj and ℓj+1.We require that the chain
not be coarsened in the atomistic regions, which precisely means that δc
= . . . = δc
whenever νj > 1.
Finally, the interactions of the atomistic energy only partially reach into the continuum part if
the atomistic potential has a finite cutoff radius. To allow for an exact calculation of this energy
without atomistic interpolation, we require that these regions are not coarsened as well. As we
will see in the next subsection, the atomistic next-nearest-neighbor interactions from the Frenkel-
Kontorova model studied in this paper reach two atoms into the continuum part. Hence, we require
that νj−2 = νj−1 = νj = νj+1 = 1 whenever δ
= 1. Other potential energies in general require
similar conditions that depend on their cut-off radius.
We denote the position of the j-th repatom by y
j = y
and the vector of all repatoms by
yqc ∈ R2N .
2.4. The Quasicontinuum Energy. Now we define the quasicontinuum energy. To this end, the
missing non-repatoms are implicitly reconstructed. We will see later that this helps to set up the
QC model, but needs not be done for the actual computation.
The reconstruction is done by a linear interpolation between the nearest repatom to the right
and to the left. That is, the vector of all atomistic positions is computed from the vector yqc of
repatom positions by
I : R2N → R2M , (Iyqc)ℓj+m :=
νj −m
j+1, m = 0, . . . , νj . (2.10)
We note that
j = (Iy
qc)ℓj . (2.11)
The underlying idea is that in regions where the lattice spacing of the atoms is nearly constant,
this interpolation is very close to the actual atomistic positions and therefore leads to a good
approximation of the total energy. Only a few repatoms are needed in these regions. This exactly
corresponds to mesh coarsening in classical finite element approximations of continuum models. On
the other hand, in regions where the lattice spacing is non-uniform, such as around a dislocation,
all atoms must be chosen to be repatoms to obtain sufficient accuracy. This guarantees that the
full resolution of the atomistic model in the critical regions is retained and corresponds to a high
refinement in classical finite element continuum models.
We define the QC approximation of the total energy to be
Eqc(yqc) := Eac(Iyqc). (2.12)
Now (2.12) has to be reformulated such that it can be computed efficiently, without the overhead
of the interpolation. Most atomistic potentials are invariant to translations, a property that al-
lows us to simplify (2.12) considerably. For any translationally invariant energy Eai , we have that
Eai (L
i,i+1ya) = φi(y
i+1 − y
i ) and E
i−1,iya) = φi(y
i − y
i−1) for some function φi. If these
functions φi coincide, that is, φi = φj for all i and j, we can write
Eai (L
i−1,i
a) = φ
yai − y
and Eai (L
i,i+1
a) = φ
yai+1 − y
(2.13)
for some function φ : R→ R. Here φ plays the role of a continuum energy density and is given for
the two-body potential (2.2) by
φ(r) =
ψ(kr).
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 5
Equations (2.7), (2.12), and (2.13) lead to
Eqc(yqc) =
i=−M+1
δai E
i (Iy
qc) + 1
i=−M+2
i−1,iIyqc) + 1
i=−M+1
δci E
i,i+1Iyqc)
i=−M+1
δai E
i (Iy
qc) + 1
i=−M+2
δciφ((Iy
qc)i − (Iy
qc)i−1) (2.14)
i=−M+1
δciφ((Iy
qc)i+1 − (Iy
qc)i).
Because Iyqc is the linear interpolation between two repatoms y
j and y
j+1, we have
(Iyqc)i+1 − (Iy
qc)i =
j+1 − y
, i = ℓj, . . . , ℓj+1 − 1. (2.15)
Hence,
Eqc(yqc) =
i=−M+1
δai E
i (Iy
qc) +
j=−N+1
j+1 − y
(2.16)
with weight factors
δcℓj + δ
0 if both y
j and y
j+1 are atomistic,
if exactly one of y
j and y
j+1 is continuum,
νj if both y
and y
j+1 are continuum.
(2.17)
The first sum corresponds to the atomistic region which will be a small region and is thus com-
putationally inexpensive. The second sum only involves at most 2N terms which is a considerable
reduction when N ≪M.
Note that the second term in formula (2.16) coincides with an integral over the energy density φ
as it occurs in finite element discretizations of classical continuum mechanical models. Hence the
apparently unmotivated definitions (2.7) and (2.5) of the continuum energy here result in what is
commonly understood as a continuum energy. The linear interpolation operator I resembles the
Cauchy-Born hypothesis.
3. Frenkel-Kontorova Model
Dislocations are lines in crystals which represent a defect in the lattice structure [15], see Fig-
ure 1. Typically, there is a core of small radius surrounding the dislocation line where the lattice
structure is highly deformed, but the lattice structure is nearly uniform outside the core. A sim-
ple one-dimensional model for a defect such as a dislocation is given by the Frenkel-Kontorova
model [15]. Here, the elastic energy is modeled by harmonic interactions between the atoms in the
one-dimensional chain and the misfit energy of the slip plane is modeled by a periodic potential. A
more accurate model of the same form is given by the Peierls-Nabarro model [11].
3.1. Atomistic Frenkel-Kontorova Model. We study a single defect in the middle of the chain
of 2M atoms. To achieve a symmetric description in terms of bonds, we number the atoms from
−M + 1 to M . The defect is situated between the atoms numbered 0 and 1 (Figure 2).
Recall that the atomistic positions are denoted by ya = (ya−M+1, . . . , y
M ) ∈ R
2M . The total
potential energy for this atomistic system is then a function Ea : R2M → R of the atomistic
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 6
Figure 1: Cross-section through a dislocation in a three-dimensional cubic lattice. The displayed plane
repeats periodically in the three-dimensional crystal. Vertical bonds are shown by lines to emphasize the
topological defect.
(−K−1)a
(K+1)a
= atomistic = continuum
Figure 2: Numbering of the atoms. The dislocation is situated in the middle of the chain between atoms ya
and ya
positions. For the Frenkel-Kontorova model, the energy, Ea = Ea,e + Ea,m, consists of two parts,
namely the part which models the elastic energy of the defect, Ea,e, and the part which models the
misfit energy on the slip plane, Ea,m.
The elastic energy is modeled by Hookean (harmonic) springs between nearest-neighbors (NN)
and next-nearest neighbors (NNN), and the total elastic energy is given by
Ea,e(ya) =
i=−M+1
i+1 − y
i − a0)
i=−M+2
i+1 − y
i−1 − 2a0)
2, (3.1)
where the moduli k1 > 0 and k2 > 0 describe the strength of the elastic interactions, and where
a0 ∈ R denotes the equilibrium distance.
We note that the asymptotic expansion to second order of any nonlinear NN/NNN potential
energy
E(ya) =
i=−M+1
ψ(yai+1 − y
i ) +
i=−M+2
ψ(yai+1 − y
i−1) (3.2)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 7
substrate
modeled
layer
misfit
energy
Figure 3: Frenkel-Kontorova model. The wells depict the misfit energy (3.4).
about aa = [(−M + 1)a0, (−M + 2)a0, · · · , (M − 1)a0,Ma0]
∈ R2M has the form
E(ya) ≈ E(aa) +
ψ′(a0) + 2ψ
′(2a0)
i=−M+1
(yai+1 − y
i − a0)
− ψ′(2a0)(y
M − y
M−1 − a0)− ψ
′(2a0)(y
−M+2 − y
−M+1 − a0)
ψ′′(a0)
i=−M+1
(yai+1 − y
i − a0)
2 + 1
ψ′′(2a0)
i=−M+2
(yai+1 − y
i−1 − 2a0)
(3.3)
We thus see that the elastic energy (3.1) with k1 = ψ
′′(a0) and k2 = ψ
′′(2a0) approximates the
energy (3.2) to second order if ψ′(a0) + 2ψ
′(2a0) = 0 and if we ignore the boundary terms in the
second line of (3.3).
The misfit energy of the slip plane is modeled by a periodic potential (Figure 3). We model this
misfit energy by
Ea,m(ya) =
i=−M+1
yai − a0
, (3.4)
where ⌊x⌋ denotes the largest integer smaller than or equal to x, and where the constant k0 > 0
determines the strength of the misfit energy.
Altogether, the total potential energy of the atomistic system is given by
Ea(ya) = Ea,e(ya) + Ea,m(ya)
i=−M+1
(yai+1 − y
i − a0)
2 + 1
i=−M+2
(yai+1 − y
i−1 − 2a0)
i=−M+1
yai − a0
(3.5)
We restrict ourselves to configurations in which the M leftmost atoms yai for −M + 1 ≤ i ≤ 0 are
situated in the interval
, whereas the M rightmost atoms yai for 1 ≤ i ≤ M
are situated in the interval
. The defect is situated between atoms ya0 and
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 8
ya1 . In this case, the total energy simplifies to
Ea(ya) = 1
i=−M+1
(yai+1 − y
i − a0)
2 + 1
i=−M+2
(yai+1 − y
i−1 − 2a0)
i=−M+1
(yai − (i− 1)a0)
(yai − ia0)
(3.6)
3.2. Quasicontinuum Approximation of the Frenkel-Kontorova Model. We now apply
the quasicontinuum method to the dislocation model described in Section 3.1. The total energy
(3.6) is split up into atom-wise contributions, separately for the elastic interactions and the misfit
interactions:
a) = 1
i − y
i−1 − a0)
2 + 1
i+1 − y
i − a0)
i − y
i−2 − 2a0)
2 + 1
i+2 − y
i − 2a0)
k0 (y
i − (i− 1)a0)
, i = −M + 1, . . . , 0,
k0 (y
i − ia0)
, i = 1, . . . ,M.
(3.7)
To simplify notation, we use the convention that the undefined terms at the endpoints of the chain
are neglected. We thus have that
Ea(ya) = Ea,m(ya) + Ea,e(ya) =
i=−M+1
a) + E
a)] . (3.8)
Since the largest displacement of the atoms is to be expected near the defect, we deem the atoms
−K +1, . . . ,K atomistic and the remaining atoms −M +1, . . . ,−K and K +1, . . . ,M continuum.
Here K < M is some constant whose optimal value will be determined by the algorithm given in
Section 5.
The optimal choice of the repatoms for coarsening is investigated in the second paper of this
series, so we work with a general formulation which holds for any values of ℓj for now. However,
there are two restrictions on the coarsening. Since the atomistic region must not be coarsened and
since we need full refinement in the vicinity of two atoms around the atomistic region due to the
NNN interactions, we have that
ℓj = j, j = −K − 1, . . . ,K + 2. (3.9)
Second, we require that
ℓ−N+1 = −M + 1, ℓ−N+2 = −M + 2, ℓN−1 =M − 1, and ℓN =M (3.10)
to incorporate the boundary conditions later.
The elastic part E
i is translationally invariant, so we perform its QC approximation as described
in the previous section. This leads to the continuum energy density
φe(r) = 1
k1(r − a0)
2 + 1
k2(2r − 2a0)
k12(r − a0)
2 (3.11)
where k12 := k1 + 4k2.
Regarding the misfit part E
i , the above technique cannot be applied since the potential is
not translationally invariant. However, there is a different summation technique to achieve a
computationally efficient formulation which avoids the costly interpolation operator.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 9
To shorten the notation, we let
indicate the sum in which the first term and the last term
are only counted half:
xi :=
i=m+1
xn (3.12)
where m < n and xi ∈ R. It is easy to verify that
2m3 +m
i(m− i) =
m3 −m
(3.13)
for m > 0.
For all pairs (j, j + 1) of continuum repatoms, we now reformulate all terms from (2.7) which
involve the interaction between ℓj and ℓj+1. For j > 0, we get by definition (2.4) of the operator
L, by definition (3.7) of E
i , and by (3.13) that
i=ℓj+1
i−1,iIyqc) + 1
ℓj+1−1
i,i+1Iyqc)
i=ℓj+1
i (Iy
qc) + 1
ℓj+1−1
i (Iy
ℓj+1 − i
i− ℓj
j+1 − ia0
ℓj+1 − i
j − ℓja0) +
i− ℓj
j+1 − ℓj+1a0)
j − ℓja0)
(ℓj+1 − i)
2 + 1
j+1 − ℓj+1a0)
(i− ℓj)
− ℓja0)(y
j+1 − ℓj+1a0)
(ℓj+1 − i)(i − ℓj)
j − ℓja0)
2ν3j + νj
j+1 − ℓj+1a0)
2ν3j + νj
j − ℓja0)(y
j+1 − ℓj+1a0)
ν3j − νj
=: φmj+(y
j , y
j+1).
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 10
For j < 0, we similarly obtain
i=ℓj+1
(Li−1,iIyqc) + 1
ℓj+1−1
(Li,i+1Iyqc)
j − (ℓj − 1)a0)
2ν3j + νj
j+1 − (ℓj+1 − 1)a0)
2ν3j + νj
j − (ℓj − 1)a0)(y
j+1 − (ℓj+1 − 1)a0)
ν3j − νj
=: φmj−(y
j , y
j+1).
Since E
−K (L
−K,−K+1Iyqc) = E
−K (Iy
qc) and E
K+1(L
K,K+1Iyqc) = E
K+1(Iy
qc), the QC approx-
imation of the chain can be given by
Eqc(yqc) =
j=−N+1
j+1 − y
+ φmj−(y
j , y
−K (Iy
qc) +
j=−K+1
Eai (Iy
qc) + 1
K+1(Iy
j+1 − y
+ φmj+(y
j , y
(3.14)
Note that the interpolation Iyqc does not have to be computed here since the relevant terms only
depend on uncoarsened parts of the chain.
Additionally, we will consider the atomistic-continuum approximation
Eac : R2M → R (3.15)
of the atomistic energy without coarsening. It is given exactly like the QC approximation (3.14)
with the only difference being that νj = 1 and ℓj = j everywhere.
4. Primal and Dual Problems
4.1. Problem Setup. We are now ready to set up the problems we will solve. We are interested
in finding the minimum of the energy (3.14) subject to given boundary conditions.
We give the boundary conditions by constraining the deformation of two atoms at each end of the
chain. This guarantees that the potential with next-nearest-neighbor interactions can be directly
applied to all non-boundary atoms without having to neglect interactions. We define the spaces
V a := R2M , V a0 := R
2M−4, V qc := R2N , V
0 := R
2N−4. (4.1)
The spaces V a and V a0 will also be used for the uncoarsened atomistic-continuum potential E
so there is no need to define spaces V ac and V ac0 . We let y
bc ∈ V a denote any vector which has
the desired boundary values ybc−M+1, y
−M+2, y
M−1, and y
M , and we let y
bcq ∈ V qc by any vector
satisfying (recall (3.10))
−N+1 = y
−M+1, y
N−1 = y
−N+2 = y
−M+2, y
= ybcM . (4.2)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 11
For any vector y ∈ V a0 , we denote the extension by zero boundary conditions to be Jy ∈ V
a, so
Jy :=
0 0 yT 0 0
∈ R2M , (4.3)
and similarly we denote the extension by zero boundary conditions of y ∈ V
0 to be J
qcy ∈ V qc.
The spaces of admissible solutions are then given by JV a0 + y
bc ⊂ V a and JqcV
0 + y
bcq ⊂ V qc,
respectively. We note that JT : V a → V a0 is the restriction operator defined by
(JTy)j = yj for j = −M + 3, . . . ,M − 2.
The minima ȳa, ȳac, and ȳqc of the energy functionals Ea, Eac, and Eqc given by (3.6), (3.15),
and (3.14) subject to the above “clamped” boundary conditions are characterized as
a := argmin
y∈JV a
Ea(y) ∈ V a, (4.4)
ac := argmin
y∈JV a
Eac(y) ∈ V a, (4.5)
qc := argmin
y∈JqcV
Eqc(y) ∈ V qc. (4.6)
We note that the minima are uniquely determined because Ea, Eac, and Eqc are strictly convex.
4.2. Matrix Formulation. For the subsequent discussion, it will be convenient to reformulate the
total energies in matrix notation:
Ea(y) = 1
(y − aa)TDaTEaDa(y − aa) + 1
(y − ba)TKa(y − ba), (4.7a)
Eac(y) = 1
(y − aa)TDaTEacDa(y − aa) + 1
(y − ba)TKa(y − ba), (4.7b)
Eqc(y) = 1
(y − aqc)TDqcTEqcDqc(y − aqc) + 1
(y − bqc)TKqc(y − bqc). (4.7c)
The matrices Da ∈ R(2M−1)×2M and Dqc ∈ R(2N−1)×2N compute the distance between two adja-
cent atomistic positions; the matrices Ea ∈ R(2M−1)×(2M−1), Eac ∈ R(2M−1)×(2M−1), and Eqc ∈
(2N−1)×(2N−1) contain the spring constants k1, k2, and k12; and the matricesK
a ∈ R(2M−1)×(2M−1)
and Kqc ∈ R(2N−1)×(2N−1) contain the misfit constant k0. The vectors a
a,ba ∈ R2M and aqc,bqc ∈
2N are constants describing the minimum energy deformations for the elastic energy and mis-
fit energy. The precise and lengthy definitions for all of these matrices and vectors are given in
Appendix A.
If we decompose ȳa = Jya + ybc for ya ∈ V a0 , then the minimization problem is given as
a = argmin
y∈V a
Ea(Jy + ybc)
= argmin
y∈V a
Jy + ybc − aa
DaTEaDa
Jy + ybc − aa
Jy + ybc − ba
Jy + ybc − ba
. (4.8)
We also decompose ȳac = Jyac + ybc and ȳqc = Jqcyqc + ybcq for yac ∈ V a0 and y
qc ∈ V
0 , and
we then formulate similar minimization problems for yac and yqc. Therefore, ya, yac, and yqc are
determined by the linear systems
Maya = fa, (4.9a)
Macyac = fac, (4.9b)
M qcyqc = f qc, (4.9c)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 12
where
Ma := JT (DaTEaDa +Ka)J,
Mac := JT (DaTEacDa +Ka)J,
M qc := JqcT (DqcTEqcDqc +Kqc)Jqc,
a := −JTDaTEaDa(ybc − aa)− JTKa(ybc − ba),
ac := −JTDaTEacDa(ybc − aa)− JTKa(ybc − ba),
qc := −JqcTDqcTEqcDqc(ybcq − aqc)− JqcTKqc(ybcq − bqc).
(4.10)
We note that the matrices Ma, Mac, and M qc are positive definite, so the total energies admit a
single global minimum and no other local minimum.
4.3. Goal-Oriented Error Estimation. To compare the approximate QC model to the original
atomistic model, we have to analyze how much the solution ya of the atomistic model deviates
from the solution yqc of the QC model. This deviation, which can be viewed as an approximation
error, can be measured in different ways, for example as ‖ya−JT IJqcyqc‖ for some norm ‖·‖. Here
we follow a different approach, namely we measure the error of a quantity of interest denoted by
Q(y) for some function Q : R2M−4 → R. Hence, we intend to estimate
Q(ya)−Q(JT IJqcyqc). (4.11)
We will assume for simplicity that Q is linear and thus has a representation Q(y) = qTy for some
vector q ∈ V a0 .
For our application, a natural quantity of interest is the size of the dislocation, that is, the
distance between the two atoms y0 and y1 to the left and right of the dislocation. This gives us
Q(y) = qTy = y1 − y0 with q = [0, . . . , 0,−1, 1, 0, . . . , 0]
T . (4.12)
Two different sources of error arise during the QC approximation, namely the localization of the
potential energy, that is, the passage from the atomistic to the continuum formulation on the one
hand, and the coarsening in the continuum region by the restriction to the repatoms on the other
hand. We denote these two errors by
e := ya − yac and eacqc := yac − JT IJqcyqc. (4.13)
It makes sense to study these sources independently. Employing the linearity of Q, we have that
|Q(ya)−Q(JT IJqcyqc)| = |Q(e) +Q(eacqc)| ≤ |Q(e)|+ |Q(eacqc)|. (4.14)
The error term |Q(e)| will be studied in Section 5, and the error term |Q(eacqc)| will be studied in
the second part of this paper series.
4.4. Dual Problems. To facilitate the goal-oriented error analysis, we introduce the dual problems
Maga = q, (4.15a)
Macgac = q, (4.15b)
M qcgqc = JqcT ITJq, (4.15c)
for ga,gac ∈ R2M−4, and gqc ∈ R2N−4. We note that the dual problems differ from the primal
problems only by the right hand side since the matrices Ma, Mac, and M qc are symmetric.
The solutions ga, gac and gqc can be viewed as influence functions: They describe how the error
at a specific point in the domain influences the error measured in terms of the goal function.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 13
Analogously to the primal errors (4.13), we define the dual errors
ê := ga − gac and êacqc := gac − JT IJqcgqc. (4.16)
In addition, we will need the primal and dual residuals
Ra(y) :=Ma (ya − y) = fa −May,
Rac(y) :=Mac (yac − y) = fac −Macy,
R̂a(g) :=Ma (ga − g) = q−Mag,
R̂ac(g) :=Mac (gac − g) = q−Macg. (4.17)
5. Error Estimation for Atomistic vs. Continuum Modeling
In this section, we estimate the error |Q(e)| arising from the approximation of an atomistic model
by a continuum model. We consider yac and gac to be computable, although in practice we can
only compute the coarsened approximations yqc and gqc.
To this end, we adapt a technique introduced in [16] and [19] to estimate the modeling error for
an elasticity model with rapidly oscillating coefficients and its homogenized version. We generalize
this technique such that it allows for different right hand sides fa and fac of the primal problem
(4.9) instead of a common right hand side as it is used in the above-mentioned works.
We have
Q(ya)−Q(yac) = qTe = gaTMae = (gacT + ê)Mae
= gacTRa(yac) + êTMae. (5.1)
The term gacTRa(yac) can be computed, whereas êTMae cannot because both e and ê are nu-
merically unknown. Instead, we estimate êTMae from above and from below by quantities that
actually can be computed.
We will give two different error estimators η1 and η2. Before, we need to derive some auxiliary
estimates to facilitate their development and analysis.
5.1. Auxiliary Estimates. We reformulate the difference ya − yac of the respective solutions in
terms of a difference of the energy matrices. To this end, we define the perturbation matrix
P := I − (Ea)−1Eac (5.2)
where I denotes the identity matrix. Note that EaP = E
a − Eac.
Lemma 5.1. For any α, β ∈ R, we have that
Ma(αe+ βê) = −JTDaTEaPDa
α(Jyac + ybc − aa) + βJgac
. (5.3)
Proof. We conclude from (4.9) that
Mae =Maya −Macyac + (Mac −Ma)yac
= fa − fac + (Mac −Ma)yac, (5.4)
and similarly since Maga =Macgac = q that
Maê =Ma(ga − gac) = (Mac −Ma)gac. (5.5)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 14
Thus, it follows from (4.10) and (5.2) that
Ma(αe+ βê) = α [(Mac −Ma)yac + fa − fac] + β(Mac −Ma)gac
= JTDaT (Eac − Ea)Da
α(Jyac + ybc − aa) + βJgac
= −JTDaTEaPDa
α(Jyac + ybc − aa) + βJgac
(5.6)
We note that the Ka-related terms cancel here, because they coincide for the atomistic model and
the continuum model. �
Lemma 5.2. We have that
‖αe+ βê‖Ma ≤
α(Jyac + ybc − aa) + βJgac
. (5.7)
We note that the right hand side is numerically computable.
Proof. To shorten the notation, we abbreviate z = α(Jyac +ybc− aa)+ βJgac. By Lemma 5.1, we
‖αe + βê‖Ma = sup
v∈V a
vTMa(αe+ βê)
‖v‖Ma
= sup
v∈V a
−vTJTDaTEaPDaz
‖v‖Ma
≤ sup
v∈V a
‖DaJv‖Ea‖PD
az‖Ea
‖DaJv‖Ea
= ‖PDaz‖Ea .
(5.8)
Here we have used that ‖DaJv‖Ea ≤ ‖v‖Ma because the matrixK
a in (4.10) is positive definite. �
5.2. First Error Estimator. We are now ready to derive the first error estimator, η1. By the
parallelogram identity, we have for all σ 6= 0 that
TMae = (σ−1êT )Ma(σe)
‖σe+ σ−1ê‖2Ma −
‖σe− σ−1ê‖2Ma .
(5.9)
In the following, we will determine computable constants η+low, η
low, η
upp and η
upp such that
η+low ≤ ‖σe+ σ
ê‖Ma ≤ η
η−low ≤ ‖σe− σ
ê‖Ma ≤ η
(5.10)
From Lemma 5.2, we immediately get the upper estimates η+upp and η
η+upp :=
σ(Jyac + ybc − aa) + σ−1Jgac
η−upp :=
σ(Jyac + ybc − aa)− σ−1Jgac
(5.11)
We note that η+low, η
low, η
upp and η
upp will depend on σ, but the estimates (5.10) will hold for any
σ 6= 0. We will now choose σ in such a way that the estimates are as sharp as possible, that is,
such that η+upp and η
upp are smallest.
Lemma 5.3. Both η+upp and η
upp given by (5.11) attain their minima for
σ̄ :=
‖PDaJgac‖Ea
‖PDa(Jyac + ybc − aa)‖Ea
. (5.12)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 15
Proof. We have that
(η±upp)
2 = σ2
∥PDa(Jyac + ybc − aa)
± 2gacTJTDaTP TEaPDa(Jyac + ybc − aa)
+ σ−2
∥PDaJgac
. (5.13)
Setting the first derivative of the mapping σ 7→ (η±upp)
2 to zero, we obtain the condition
∥PDa(Jyac + ybc − aa)
− 2σ̄−3
∥PDaJgac
= 0 (5.14)
for critical points of (η±upp)
2. This equation has the unique positive solution (5.12). Because
lim|σ|→∞ η
upp = limσ→0 η
upp =∞, this point corresponds to a minimum. Hence the quantities η
attain their minima at σ = σ̄. �
Regarding the lower bounds η+low and η
low, we have
‖σ̄e± σ̄−1ê‖Ma = sup
v∈V a
vTMa(σ̄e± σ̄−1ê)
‖v‖Ma
= sup
v∈V a
vT (σ̄Ra(yac)± σ̄−1R̂a(gac))
‖v‖Ma
vT0 (σ̄R
a(yac)± σ̄−1R̂a(gac))
‖v0‖Ma
(5.15)
for any vector v0 ∈ V
0 \ {0}. Numerically, we have the two vectors y
ac and gac at our disposal,
hence it makes sense to take a linear combination v0 = y
ac + θ±gac. Here we follow the strategy
of [19] and choose θ± as the critical points of η±low.
Lemma 5.4. Let
± = σ̄Ra(yac)± σ̄−1R̂a(gac). (5.16)
Then the lower bounds
(yac + θ±gac)T r±
‖yac + θ±gac‖Ma
(5.17)
have a unique critical point for
θ̄± :=
r±Tyac gacTMayac − r±Tgac ‖yac‖2Ma
r±Tgac gacTMayac − r±Tyac ‖gac‖2
. (5.18)
Proof. We have
r±Tgac‖yac + θ±gac‖Ma − r
±T (yac + θ±gac)
(yac+θ±gac)TMagac
‖yac+θ±gac‖Ma
‖yac + θ±gac‖2
. (5.19)
Setting this expression to zero and solving for θ± leads to the above condition. �
However, let us note that this critical point is not necessarily a maximum of η±low, which would
be optimal for bound (5.10). Depending on the actual vectors yac and gac, it can be shown that
this critical point could be a minimum.
Now we have all necessary ingredients to construct the error estimator η1. From (5.1) and (5.9),
we get the computable estimate
gacTRa(yac) + 1
)2 − 1
(η−upp)
2 ≤ Q(ya)−Q(yac) ≤ gacTRa(yac) + 1
(η+upp)
2 − 1
(5.20)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 16
At first sight, this looks like we could get an estimate for |Q(ya)−Q(yac)| from both above and
below. However, this is only true if both the left hand side and the right hand side have the same
sign, which in general does not hold. But we get the following estimate.
Theorem 5.1. We have that
|Q(ya)−Q(yac)| ≤ η1, (5.21)
where the computable error estimator is defined as
η1 := max
∣gacTRa(yac) + 1
)2 − 1
(η−upp)
∣gacTRa(yac) + 1
(η+upp)
2 − 1
. (5.22)
We note that the computation of the η±upp terms involves the solution of a linear system with
matrix Ea as its inverse appears in the operator P . The matrix Ea is not diagonal, but has
condition number O(1). So this is negligible compared to what would be necessary to solve the
original atomistic problem which includes the operator DaTEaDa with condition number O(M2).
5.3. Second Error Estimator. There is no reasonable way to decompose the error estimator η1
into a sum of element-wise or atom-wise contributions due to the η±
terms. Therefore, we derive
another error estimator η2 which allows for such a decomposition, at the price of a less accurate
estimate than η1.
Theorem 5.2. We have that
|Q(ya)−Q(yac)| ≤ η2 ≤
i=−M+3
ηat2,i +
i=−M+1
ηel2,i (5.23)
where the computable global error estimator η2 and the computable local error estimators, η
2,i and
ηel2,i, associated with atoms and elements, respectively, are defined as
η2 :=
∣gacTRa(yac)
∣+ ‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea,
ηat2,i := |g
a(yac)i| , i = −M + 3, . . . ,M − 2,
ηel2,i :=
PDa(Jyac + ybc − aa)
(Ea − Eac)Da(Jyac + ybc − aa)
∣(PDaJgac)i
(Ea − Eac)DaJgac
∣ , i = −M + 1, . . . ,M − 1.
(5.24)
Proof. From (5.1) and Lemma 5.2, we conclude that
|Q(ya)−Q(yac)| ≤
acTRa(yac)
acTRa(yac)
∣+ ‖ê‖Ma‖e‖Ma
acTRa(yac)
∣+ ‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea
= η2,
(5.25)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 17
which gives us the global estimate. For the decomposition into local contributions, we further
estimate
acTRa(yac)
∣+ ‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea
acTRa(yac)
‖PDa(Jyac + ybc − aa)‖2Ea +
‖PDaJgac‖2Ea
i=−M+3
|gaci R
a(yac)i|+
i=−M+1
∣(PDaJgac)i
(Ea − Eac)DaJgac
i=−M+1
PDa(Jyac + ybc − aa)
(Ea − Eac)Da(Jyac + ybc − aa)
i=−M+3
ηat2,i +
i=−M+1
ηel2,i,
(5.26)
which completes the proof. �
Let us remark that instead of the first inequality in (5.26), one can get an apparently better
estimate
‖PDa(Jyac + ybc − aa)‖Ea‖PD
aJgac‖Ea
γ‖PDa(Jyac + ybc − aa)‖2Ea +
γ−1‖PDaJgac‖2Ea
(5.27)
by introducing the additional weight factor
‖PDaJgac‖Ea
‖PDa(Jyac + ybc − aa)‖Ea
, (5.28)
and then decomposing the resulting terms similar to the above. However, our numerical results
showed that this modification does not significantly improve the decomposed error estimator for
the application considered here.
6. Numerics
In the preceding sections, we constructed the error estimators η1 and η2. We will now give an
algorithm for adaptive atomistic-continuum modeling based on these error estimators. Then we
will present and discuss some numerical results.
6.1. Algorithm. The error estimator η1 should give a better estimate of the error than η2, because
η2 involves the inequality
∣êTMae
∣ ≤ ‖ê‖Ma‖e‖Ma in (5.25) in contrast to the parallelogram identity
for η1. However, η2 can be decomposed into atom-wise and element-wise contributions η
2,i and
ηel2,i, whereas the η
terms in η1 do not admit a reasonable decomposition that can be used for
atomistic-continuum adaptivity.
We make use of this by employing the sharper estimate η1 to determine whether a given global
error tolerance τgl for the error in an adaptive algorithm has already been achieved or not. If
not, we use the decomposed estimates ηat2,i and η
2,i to determine where the more precise atomistic
modeling is needed. This leads us to the following algorithm:
(1) Choose τgl. Model all atoms as a continuum. Set τat ← τgl.
(2) Solve primal problem (4.9b) for yac and dual problem (4.15b) for gac.
(3) Compute error estimator η1 from (5.22).
(4) If η1 ≤ τgl, then stop.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 18
M iteration K τat η1
100 1 0 1.000000e-10 3.899207e-02
2 28 1.000000e-11 5.915080e-10
3 32 1.000000e-12 4.878532e-11
1000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.915100e-10
3 32 1.000000e-12 4.878548e-11
10000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.915100e-10
3 32 1.000000e-12 4.878548e-11
100000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.915099e-10
3 32 1.000000e-12 4.878540e-11
1000000 1 0 1.000000e-10 3.899208e-02
2 28 1.000000e-11 5.914422e-10
3 32 1.000000e-12 4.871775e-11
Table 1: Convergence of the algorithm for τgl = 10
−10 and different values of M .
(5) Compute local error estimators ηat2,i and η
2,i from (5.24).
(6) Set τat ←
(7) Make all atoms i atomistic for which
ηtot2,i := η
2,i +
ηel2,i−1 + η
≥ τat. (6.1)
(8) Go to (2).
Here τdiv > 1 is a constant factor which describes how fast the atom-wise tolerance τat should
decrease during adaption. Our experience has been that τdiv = 10 is a reasonable choice.
The crucial adaption step is (7). The adaption criterion (6.1) deems all atoms to be modeled
atomistically if the associated error from the decomposition of η2 exceeds the atomistic error toler-
ance τat. Here the element-wise errors η
2,i are distributed equally to the two adjacent atoms i and
i+ 1.
For the dislocation at the center of the chain and the chosen goal function, we expect that the
atomistic repatoms always form a symmetric interval around the center. We have used the above
adaptive atomistic-continuum algorithm to approximate our Frenkel-Kontorova model and have
always found that the atomistic region is the set of atoms −K + 1, . . . ,K for some K depending
on M and τgl. Thus, the modeling approach given in Section 3 of restricting to an atomistic region
consisting of atoms −K + 1, . . . ,K for some K rather than considering a more general atomistic
region is justified a posteriori.
6.2. Numerical Results. The algorithm has been implemented as described above. The boundary
conditions were chosen as
ybc−M+1 = −M, y
−M+2 = −M + 1, y
M−1 =M − 1, y
M =M. (6.2)
The elastic constants are k0 = 1 and k1 = k2 = 2.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 19
−500 −400 −300 −200 −100 0 100 200 300 400 500
−500 −400 −300 −200 −100 0 100 200 300 400 500
Figure 4: The error estimators ηtot
2,i (left) and η
2,i, η
2,i (right) for M = 500, K = 20.
Table 1 shows how the algorithm given above performs. After 3 iterations, the desired accuracy
τgl = 10
−10 is achieved. Moreover, we can see from the table that the number of iterations are
independent of M , that means the algorithm behaves robustly with respect to the problem size M .
Figure 4 (left) shows the decomposition of the error estimator η2 for a typical setting M = 500,
K = 20. One can clearly see that the error in the atomistic region is small, whereas the error is
large in the continuum regions that border the atomistic region. It then decreases exponentially
towards the endpoints.
The error in both the atomistic region around the center and the continuum regions far away
from the center are in the range of the (relative) machine precision εmach, which accounts for the
fluctuations in these regions. The error can be considered to be numerically zero in these regions.
In the continuum regions, we observe an error of magnitude O(ε2
), whereas in the continuum
region we have O(εmach), which leads to the different magnitudes of the fluctuations.
Figure 4 (right) shows the element-wise contributions ηel2,i and the atom-wise contributions η
of the decomposed error estimator ηtot2,i = η
2,i +
ηel2,i−1 + η
. The atomistic part ηat2,i, which
corresponds to the gacTRa(yac) term, is dominant in the sense that it is about ten times larger
than ηel2,i, which comes from the estimate for the perturbation term ê
TMae. The fluctuations
due to the limited machine precion in the atomistic region come from ηat2,i, whereas those in the
continuum region away from the defect stem from ηel2,i. Let us note that in other applications of
duality-based error estimation, the first term might not always be the dominant term. For example,
in mesh refinement for classical linear finite elements, the first term even vanishes due to Galerkin
orthogonality.
Table 2 and Figure 5, which display the same data, show the efficiency of the error estimators,
η1/|Q(y
a − yac)| and η2/|Q(y
a − yac)|, for M = 1000. For comparison, the actual error is given as
well. For the relatively small 1D problem, the actual error can be easily computed, whereas in real
world applications it is of course not available. One can clearly see that η1 gives a better estimate
than η2, which numerically confirms our conjecture that η1 is a better estimator than η2. We see
that η1 overestimates the actual error by a factor of 1.4, while η2 is in a still acceptable range of
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 20
K |Q(ya − yac)| η1 η1/|Q(y
a − yac)| η2 η2/|Q(y
a − yac)|
0 3.627633e-02 3.899208e-02 1.074863 3.999783e-02 1.102588
2 3.375762e-02 3.872272e-02 1.147081 5.101700e-02 1.511274
4 3.468605e-03 4.343595e-03 1.252260 5.422007e-03 1.563166
6 5.418585e-04 7.156249e-04 1.320686 9.187940e-04 1.695635
8 1.227067e-04 1.675383e-04 1.365356 2.193196e-04 1.787348
10 3.287188e-05 4.540984e-05 1.381419 5.984186e-05 1.820457
15 1.416914e-06 1.966114e-06 1.387603 2.597488e-06 1.833201
20 6.267636e-08 8.695824e-08 1.387417 1.148736e-07 1.832805
25 2.770161e-09 3.843388e-09 1.387424 5.077204e-09 1.832819
30 1.224369e-10 1.698739e-10 1.387440 2.244073e-10 1.832840
35 5.410783e-12 7.508365e-12 1.387667 9.918687e-12 1.833133
40 2.379208e-13 3.318024e-13 1.394592 4.383361e-13 1.842362
45 8.992806e-15 1.430733e-14 1.590975 1.901601e-14 2.114580
50 7.771561e-16 4.120094e-16 0.530150 6.201285e-16 0.797946
Table 2: Efficiency of the error estimators, η1/|Q(y
a − yac)| and η2/|Q(y
a − yac)|, for M = 1000. For
K = 45 and K = 50 the results become inaccurate due to limited machine precision.
0 5 10 15 20 25 30 35 40 45 50
|Q(ya−yac)|
Figure 5: Efficiency of the error estimators for M = 1000.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 21
τgl optimal K K by η1 K by η2
1e-02 3 3 3
1e-03 5 5 5
1e-04 9 9 10
1e-05 12 13 13
1e-06 16 17 17
1e-07 20 20 21
1e-08 23 24 24
1e-09 27 28 28
1e-10 31 31 32
1e-11 35 35 35
1e-12 38 39 39
1e-13 42 42 43
1e-14 45 46 47
Table 3: Efficiency of the error estimators for M = 1000.
2 times the actual error. Moreover, we can see from Table 2 and Figure 5 that the error decreases
exponentially with K.
Finally, we compare the optimal (smallest) value of K which is needed to achieve a given accuracy
τgl with the values for K determined by the error estimators η1 and η2, again taking into account
the precise error which is available for the model problem. We see from Table 3 that even η2 only
overestimates K by at most 2 atoms. Thus, we get an efficient estimate of the required atomistic
region for both error estimators.
Appendix A. Matrix Definitions
We describe the matrices from Section 4.2. The matrix
. . .
. . .
∈ R(2M−1)×2M (A.1)
transforms atomistic positions to distances between adjacent atoms. Similarly,
Dqc =
−ν−1−N+1 ν
−ν−1−N+2 ν
. . .
. . .
N−1 ν
∈ R(2N−1)×2N (A.2)
transforms repatom positions from a coarsened chain to normalized distances between adjacent
repatoms. The matrices
(Ea)ij =
k1 + k2 i = j ∈ {−M + 1,M − 1}
k1 + 2k2 i = j ∈ {−M + 2, . . . ,M − 2}
k2 |j − i| = 1
0 otherwise,
(A.3)
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 22
(Eac)ij =
k12(δ
i + δ
i+1) +
i + δ
i+1) +
i−1 + δ
i + δ
i+1 + δ
i+2) i = j
i + δ
i+2) j = i+ 1
i−1 + δ
i+1) j = i− 1
0 otherwise,
(A.4)
(Eqc)ij =
ωik12 +
i + δ
i+1) +
i−1 + δ
i + δ
i+1 + δ
i+2) i = j
i + δ
i+2) j = i+ 1
i−1 + δ
i+1) j = i− 1
0 otherwise,
(A.5)
for i, j = −M + 1, . . . ,M − 1 and i, j = −N + 1, . . . , N − 1, respectively, describe the spring
interactions in terms of the distances between atoms or repatoms. Accordingly, the matrices
. . .
∈ R2M×2M (A.6)
(Kqc)ij =
(2νi−1 + ν
i−1) + (2νi + ν
i = j ∈ {−N + 2, . . . , N − 1}
k0(2ν−N+1 + ν
−N+1) i = j = −N + 1
k0(2νN−1 + ν
N−1) i = j = N
k0(νi − ν
i ) j = i+ 1
k0(νj − ν
j ) j = i− 1
0 otherwise,
(A.7)
for i, j = −N + 1, . . . , N describe the misfit interactions for the original atomistic system and the
QC approximation. Finally, the constant vectors
(−M + 1)a0 (−M + 2)a0 · · · (M − 1)a0 Ma0
∈ R2M , (A.8a)
ℓ−N+1a0 ℓ−N+2a0 · · · ℓN−1a0 ℓNa0
∈ R2N , (A.8b)
−Ma0 (−M + 1)a0 · · · −a0 a0 · · · (M − 1)a0 Ma0
∈ R2M , (A.8c)
(ℓ−N+1 − 1)a0 (ℓ−N+2 − 1)a0 · · · (ℓ0 − 1)a0 ℓ1a0 · · · ℓN−1a0 ℓNa0
∈ R2M ,
(A.8d)
fix the equilibrium positions for the spring interactions and the misfit energy, respectively.
References
[1] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York,
2000.
[2] M. Arndt and M. Griebel, Derivation of higher order gradient continuum models from atomistic models for
crystalline solids, Multiscale Model. Simul., 4 (2005), pp. 531–562.
[3] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Lectures in
Mathematics, ETH Zürich, Birkhäuser, Basel, 2003.
[4] X. Blanc, C. Le Bris, and F. Legoll, Analysis of a prototypical multiscale method coupling atomistic and
continuum mechanics, Math. Model. Numer. Anal., 39 (2005), pp. 797–826.
[5] P. M. Chaikin and T. C. Lubensky, Principles of Condensed Matter Physics, Cambridge University Press,
2000.
ADAPTIVITY FOR THE QUASICONTINUUM APPROXIMATION OF A FRENKEL-KONTOROVA MODEL 23
[6] M. S. Daw and M. I. Baskes, Semiempirical, quantum mechanical calculation of hydrogen embrittlement in
metals, Phys. Rev. Lett., 50 (1983), pp. 1285–1288.
[7] , Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals,
Phys. Rev. B, 29 (1984), pp. 6443–6453.
[8] M. Dobson and M. Luskin, Analysis of a force-based quasicontinuum approximation, 2006,
arXiv:math.NA/0611543.
[9] W. E, J. Lu, and J. Z. Yang, Uniform accuracy of the quasicontinuum method, Phys. Rev. B, 74 (2006),
p. 214115.
[10] W. E and P. Ming, Analysis of multiscale methods, J. Comput. Math., 22 (2004), pp. 210–219.
[11] E. Kaxiras, Atomic and Electronic Structure of Solids, Cambridge University Press, 2003.
[12] J. Knap and M. Ortiz, An analysis of the quasicontinuum method, J. Mech. Phys. Solids, 49 (2001), pp. 1899–
1923.
[13] P. Lin, Theoretical and numerical analysis for the quasi-continuum approximation of a material particle model,
Math. Comput., 72 (2003), pp. 657–675.
[14] , Convergence analysis of a quasi-continuum approximation for a two-dimensional material without defects,
SIAM J. Numer. Anal., 45 (2007), pp. 313–332.
[15] M. Marder, Condensed Matter Physics, John Wiley & Sons, 2000.
[16] J. T. Oden and S. Prudhomme, Estimation of modeling error in computational mechanics, J. Comput. Phys.,
182 (2002), pp. 496–515.
[17] J. T. Oden, S. Prudhomme, and P. Bauman, Error control for molecular statics problems, Int. J. Multiscale
Comput. Eng., 4 (2006), pp. 647–662.
[18] J. T. Oden, S. Prudhomme, A. Romkes, and P. Bauman, Multiscale modeling of physical phenomena:
Adaptive control of models, SIAM J. Sci. Comput., 28 (2006), pp. 2359–2389.
[19] J. T. Oden and K. S. Vemaganti, Estimation of local modeling error and goal-oriented adaptive modeling
of heterogeneous materials: Part I: Error estimates and adaptive algorithms, J. Comput. Phys., 164 (2000),
pp. 22–47.
[20] C. Ortner and E. Süli, A-posteriori analysis and adaptive algorithms for the quasicontinuum method in one
dimension, Research Report NA-06/13, Oxford University Computing Laboratory, 2006.
[21] , A-priori analysis of the quasicontinuum method in one dimension, Research Report NA-06/12, Oxford
University Computing Laboratory, 2006.
[22] E. B. Tadmor, R. Miller, R. Phillips, and M. Ortiz, Nanoindentation and incipient plasticity, J. Mater.
Res., 14 (1999), pp. 2233–2250.
[23] E. B. Tadmor, M. Ortiz, and R. Phillips, Quasicontinuum analysis of defects in solids, Philos. Mag. A, 73
(1996), pp. 1529–1563.
[24] E. B. Tadmor, R. Phillips, and M. Ortiz, Mixed atomistic and continuum models of deformation in solids,
Langmuir, 12 (1996), pp. 4529–4534.
	1. Introduction
	2. Quasicontinuum Approximation
	2.1. The Atomistic System
	2.2. The Atomistic-Continuum Energy
	2.3. Repatoms: Reduction of Degrees of Freedom
	2.4. The Quasicontinuum Energy
	3. Frenkel-Kontorova Model
	3.1. Atomistic Frenkel-Kontorova Model
	3.2. Quasicontinuum Approximation of the Frenkel-Kontorova Model
	4. Primal and Dual Problems
	4.1. Problem Setup
	4.2. Matrix Formulation
	4.3. Goal-Oriented Error Estimation
	4.4. Dual Problems
	5. Error Estimation for Atomistic vs. Continuum Modeling
	5.1. Auxiliary Estimates
	5.2. First Error Estimator
	5.3. Second Error Estimator
	6. Numerics
	6.1. Algorithm
	6.2. Numerical Results
	Appendix A. Matrix Definitions
	References
