# Efficient a Posteriori Error Control of a Consistent Atomistic/Continuum
Coupling Method for Two Dimensional Crystalline Defects

Yangshuai Wang  Hao Wang<EMAIL_ADDRESS>Department of Mathematics,
University of British Columbia, 1984 Mathematics Road, Vancouver, Canada
School of Mathematics, Sichuan University, No. 24 Yihuan Road, Chengdu, China

###### Abstract

Adaptive atomistic/continuum (a/c) coupling method is an important method for
the simulation of material and atomistic systems with defects to achieve the
balance of accuracy and efficiency. Residual based a posteriori error
estimator is often employed in the adaptive algorithm to provide an estimate
of the error of the strain committed by applying the continuum approximation
for the atomistic system and the finite element discretization in the
continuum region. In this work, we propose a theory based approximation for
the residual based a posteriori error estimator which greatly improves the
efficiency of the adaptivity. In particular, the numerically expensive
modeling residual is only computed exactly in a small region around the
coupling interface but replaced by a theoretically justified approximation by
the coarsening residual outside that region. We present a range of adaptive
computations based on our modified a posteriori error estimator and its
variants for different types of crystalline defects some of which are not
considered in previous related literature of the adaptive a/c methods. The
numerical results show that, compared with the original residual based error
estimator, the adaptive algorithm using the modified error estimator with
properly chosen parameters leads to the same optimal convergence rate of the
error but reduces the computational cost by one order with respect to the
number of degrees of freedom.

###### keywords:

Atomistic-to-continuum coupling, crystalline defects, a posteriori error
estimate, adaptivity.

###### MSC:

[2020] 65N12, 65N15, 70C20, 82D25

##  1 Introduction

Atomistic/continuum (a/c) coupling methods are a class of concurrent multi-
scale schemes coupling molecular mechanics models of atomistic processes with
continuum models of long-ranged elastic fields [1, 2, 3, 4, 5, 6, 7]. The goal
of the a/c methods is to combine the accuracy of the atomistic models and the
efficiency of the continuum models for the simulation of atomistic systems,
especially those with defects. In particular, the atomistic model is often
applied in a small neighborhood of the localized defects such as vacancies and
dislocations, while the continuum model is employed away from the defect
cores. We refer to [3] for an extensive overview and benchmarking of different
a/c methods.

In the past two decades, a/c coupling methods have attracted considerable
attention from both the engineering community and the mathematical community
[3, 8, 9, 10, 11, 12, 13]. The quantitative estimates of the approximation
error for the a/c coupling methods motivate and establish analytical
frameworks for similar multiscale computational methods [14, 15, 16, 17]. The
rigorous mathematical and numerical analysis for the a/c coupling methods
provides the theoretical foundation for achieving (quasi-)optimal balance
between accuracy and computational cost for the simulation of various
crystalline defects [6, 7, 9, 10, 18, 19, 20, 21].

While the modeling and the a priori analysis of a/c coupling methods have
received extensive and comprehensive discussion by numerous previous works [3,
6, 9, 19, 22], only a few researches concern the a posteriori analysis and
adaptivity which is the cornerstone for the applications of the methods in
real world material system.

The goal-oriented a posteriori error estimate is the first approach applied to
the adaptive a/c problems. Oden and Prudhomme propose a general framework for
the estimation of modeling error in computational mechanics in [23] which is
later applied to the adaptive computation of a practical three dimensional
nanoindentation problem in [24] together with a fine discussion and an
extension to the problem of molecular dynamics in [14]. Arndt and Luskin
derive the goal-oriented error estimator for the one-dimensional Frenkel-
Kontorova model in a series of works [25, 26, 27] where a detailed study of
the model adaptivity and the mesh adaptivity using the goal-oriented approach
is given. However, the a/c methods employed in the works just mentioned are
the original energy-based quasicontinuum method (QCE) which, despite its
comparatively easy implementation, is inconsistent and suffers from the ghost-
force at the atomistic to continuum coupling interface [3, 19, 28]. More
recent advances in this direction are given in [11] and [29] where an adaptive
modeling error control for an a/c method based on the Arlequin framework [30]
are analyzed. The advantage of the goal-oriented a posteriori error estimate
is that it concentrates on the minimization of the error of some prescribed
quantities of interests which often involve certain geometric parameters of
the defect core, though the exact value of such quantities may not be easy to
estimate so that a proper stopping criterion for the adaptive algorithm is yet
to be proposed [11]. In contrast, the residual based a posteriori error
estimate, which essentially follows the classic adaptive finite element
methods, provides a quantitative estimate of the error in certain global norms
that are related more directly to the strain and the energy of the system [31,
32, 33]. The residual based a posteriori error control is analyzed in [31, 32]
for consistent a/c coupling methods in one-dimension. The extension to two
dimensional system with nearest neighbor interaction is first carried out in
[33] where rigorous a posteriori estimates of the residual and the stability
constant are derived. The analysis in [33] also clarifies the so-called
“stress tensor correction” which essentially distinguishes the high
dimensional results with previous one dimensional ones. With the a posteriori
error estimates and the corresponding adaptive algorithm, the numerical
experiments in [33] achieve the adaptive movement of the a/c interface and the
mesh refinement in the continuum region. The two dimensional analysis is then
extended to the case of finite range interactions in [34] by introducing a
computable approximated stress tensor formulation. However, both of the two
works consider only simple point defects (di-vacancy in [33] while di-vacancy
and micro-crack in [34]). More importantly, the main drawback of the original
residual based a posteriori error estimates proposed in these two works is the
inefficiency of the resulting adaptive algorithms which comes from the high
computational cost of the evaluation of the modeling residual. The numerical
experiments in [34] find that the CPU time of computing the modeling residual
dramatically exceeds that of solving the a/c coupling problem.

The purpose of this work is to develop an efficient adaptive algorithm based
on a theory based approximation for the residual based a posteriori error
estimator for a consistent a/c coupling scheme and apply the algorithm to the
adaptive simulations of different defects in two dimensions. In particular, we
prove that the modeling residual, which is expensive to compute, and the
coarsening residual, whose computational cost is marginal compared with that
of solving the a/c problem on a given mesh, essentially follow a ratio
relation under certain assumptions. Based on this theoretical result, we then
propose various modified residual based a posteriori error estimator and mesh
structures so that the modeling residual is only exactly computed in a
relatively small neighborhood of the coupling interface but properly
approximated by the coarsening error away from the defect core. We finally
conduct a range of adaptive computations, using the modified estimators as
well as the original residual based estimator and a priori graded meshes if
applicable, for the three typical types of defects, namely micro-crack, anti-
plane screw dislocation and multiple vacancies, where the later two are
considered for the first time in the research of the adaptive a/c methods. The
numerical results show that our modified error estimator with properly chosen
parameters leads to the same convergence behavior as the original error
estimator while gain the efficiency by reducing the computational cost by one
order with respect to the number of degrees of freedom. In addition, the
advantage of the adaptive a/c methods is clearly demonstrated in the case of
multiple vacancies where no a priori graded mesh may be derived at first. Such
defect can also be considered as a simplification of the more complicated but
realistic models which consider the interactions of multiple defects where the
adaptive a/c may better reveal its power.

We constrain ourselves to the GRAC23 derived in [21] with the nearest neighbor
interactions to clearly present the ideas and steps, though the analysis and
constructions can be in principle extended to finite range interactions and to
other a/c methods. We give a further remark for this in § 6.

###  1.1 Outline

In § 2 we introduce the crystalline defects and set up the atomistic, the
continuum and the geometric reconstruction based atomistic/continuum (GRAC)
coupling scheme considered in this work. In § 3 we first review the a priori
and a posteriori error estimates for the a/c coupling scheme previously given
in [21, 33, 35]. We then explain where the high computational cost previously
noticed in [34] comes and present the theoretical analysis as well as the
numerical justification of the relation between the modeling and the
coarsening residual estimators. In § 4 we propose the modified error estimator
and the corresponding adaptive algorithm based on the theoretical and
numerical results in the previous section. In § 5 we present the numerical
results for the adaptive computations for the crystalline defects we consider
and give a thorough discussion and explanation. We draw conclusions and
discuss possible future works in § 6.

###  1.2 Notations

We use the symbol $\langle\cdot,\cdot\rangle$ to denote an abstract duality
pairing between a Banach space and its dual space. The symbol $|\cdot|$
normally denotes the Euclidean or Frobenius norm, while $\|\cdot\|$ denotes an
operator norm. For $E\in C^{2}(X)$, the first and second variations are
denoted by $\langle\delta E(u),v\rangle$ and $\langle\delta^{2}E(u)v,w\rangle$
for $u,v,w\in X$. The closed ball with radius $r>0$ and center $x$ is denoted
by $B_{r}(x)$, or $B_{r}$ if the center is the origin.

##  2 Formulation of the Atomistic Model and its A/C Approximation

We present the atomistic fomulation and its a/c approximation of crystalline
defects in an two-dimensional infinite lattice in this section. We keep the
presentation as concise as possible since much of the detail can be found in
various earlier works [10, 21, 33, 34, 35].

###  2.1 Atomistic model

We denote a perfect single lattice which possesses no defects by $\Lambda^{\rm
hom}:=\mathsf{A}\mathbb{Z}^{2}$ where $\mathsf{A}\in\mathbb{R}^{2\times 2}$
can be considered as a macroscopic deformation gradient that often represents
a uniform stretch or compression of a crystalline. To model local defects of a
crystal lattice, we employ $\Lambda\subset\mathbb{R}^{2}$ which satisfies the
following assumption:

######  Assumption 1.

$\exists~{}\Omega^{\rm def}:=\cup_{i=1}^{M}B_{R^{\rm def}_{i}}(x_{i})$, such
that $\Lambda\backslash\Omega^{\rm def}=\Lambda^{\rm hom}\backslash\Omega^{\rm
def}$ and $\Lambda\cap\Omega^{\rm def}$ is finite.

The region $B_{R^{\rm def}_{i}}(x_{i})$ often contains localized defects with
a center positioned at $x_{i}$ and is called a defect core. Such formulation
can represent different types of defects. In particular, for point defects
such as vacancies and micro-cracks we have $\Lambda\subset\Lambda^{\rm hom}$,
and for an anti-plane screw dislocation we have $\Lambda\equiv\Lambda^{\rm
hom}$.

We denote a deformed configuration by a lattice function $y\in\mathscr{U}$
where $\mathscr{U}:=\\{v:\Lambda\to\mathbb{R}^{2}\\}$ is the set of vector-
valued lattice functions. The displacement $u\in\mathscr{U}$ is then defined
accordingly by
$u(\ell)=y(\ell)-y_{0}(\ell)=y(\ell)-x(\ell)-u_{0}(\ell)=y(\ell)-\ell-
u_{0}(\ell)$ for any $\ell\in\Lambda$ if we define $x\in\mathscr{U}$ to be the
identity map. We consider $u_{0}\in\mathscr{U}$ as the predictor, which is
usually a far-field crystalline environment, and $u\in\mathscr{U}$ the
corrector [8, 36]. We simply take $u_{0}=0$ for point defects and the
derivation of $u_{0}$ for anti-plane screw dislocation is reviewed in Appendix
A. We refer to Figure 1 for the illustration of the geometry of the three
different types of defects we consider in this work, namely micro-crack, anti-
plane screw dislocation and multiple vacancies.

(a) Micro-crack

(b) Anti-plane screw dislocation (colors by predictor $u_{0}$)

(c) Multiple separate vacancies

Figure 1: Illustration of three typical defects considered in this work.

######  Remark 1.

Assumption 1 is presented in a slightly different form as that in various
earlier works such as [8, 33, 36, 37], where only a single round shaped defect
core $B_{R^{\rm def}}$ is assumed. We decide to use the current formulation to
include the case of multiple vacancies which has not been considered before.
In a/c adaptive computations, the defect core is often contained in the
atomistic region and the allowance of separating the defect core into a number
of round shaped regions $B_{R^{\rm def}_{i}}(x_{i})$ facilitates the
description of the initial steps of the adaptivity which are clearly
visualized in Figure 11 in § 5.3. However, we note that whether the defects
are separated or “centralized” have no fundamental influence on the analytical
results we derive in the current work. Only the condition that
$\Lambda\cap\Omega^{{\rm def}}=\Lambda\cap\cup_{i=1}^{M}B_{R^{\rm
def}_{i}}(x_{i})$ (or $\Lambda\cap B_{R^{\rm def}}$ in the case of a single
round shaped defect core) plays the role which essentially means that the
defect core can not appear at infinity.

We define the interaction neighborhood for each $\ell\in\Lambda$ by
$\mathcal{N}_{\ell}:=\\{\ell^{\prime}\in\Lambda\,|\,0<|\ell^{\prime}-\ell|\leq
r_{\rm cut}\\}$, where $r_{\rm cut}$ is called the cut-off radius and is
always equal to $1$ in the current work since we only consider the nearest
neighbor interaction. We also denote the interaction range
$\mathcal{R}_{\ell}:=\\{\ell^{\prime}-\ell\,|\,\ell^{\prime}\in\mathcal{N}_{\ell}\\}$
as the union of lattice vectors defined by the finite differences between the
atom $\ell$ and the lattice points in its neighborhood $\mathcal{N}_{\ell}$.

The finite difference stencil for each $\ell\in\Lambda$ is then defined by
$Dv(\ell):=\\{D_{\rho}v(\ell)\\}_{\rho\in\mathcal{R}_{\ell}}:=\\{v(\ell+\rho)-v(\ell)\\}_{\rho\in\mathcal{R}_{\ell}}$.
The lattice energy-norm (essentially a discrete $H^{1}$-semi-norm), which we
use to measure the difference or error for lattice functions, is then defined
by

$\|Dv\|_{\ell^{2}}:=\bigg{(}\sum_{\ell\in\Lambda}\sum_{\rho\in\mathcal{R}_{\ell}}|D_{\rho}v(\ell)|^{2}\bigg{)}^{1/2},\
\ \forall v\in\mathscr{U}.$ (1)

To measure the local “regularity” of a displacement function
$u\in\mathscr{U}$, it is convenient to assume the existence of a background
mesh, that is, a canonical triangulation $\mathcal{T}_{{\rm a}}$ of
$\mathbb{R}^{2}$ into triangles whose nodes are the reference lattice sites
$\Lambda$. We refer to [8, Figure 1] for the illustration of such
triangulation in 2D. We then introduce a standard piecewise affine
interpolation of $u$ with respect to $\mathcal{T}_{{\rm a}}$ by $I_{{\rm
a}}u$. Notice that $\Lambda$ is a lattice with some vacancies when we consider
micro-crack and multiple vacancies. We can thus construct the piecewise
interpolant with respect to $\Lambda^{\rm hom}$ by extending $u$ to the
vacancy sites (see [33, Appendix 1] for detail).

We identify $u=I_{\rm a}u$ if there is no confusion and denote the piecewise
constant gradient $\nabla u:=\nabla I_{\rm a}u$. We then introduce the
discrete homogeneous Sobolev space

$\mathscr{U}^{1,2}:=\\{u\in\mathscr{U}~{}|~{}\nabla u\in L^{2}\\},$

with the semi-norm $|u|_{\mathscr{U}^{1,2}}:=\|\nabla u\|_{L^{2}}$.

The atomistic system we consider is modeled by the site energy $V_{\ell}$
associated with the lattice points where
$V_{\ell}:(\mathbb{R}^{2})^{\mathcal{R}_{\ell}}\rightarrow\mathbb{R}$
represents the energy contribution by atom $\ell$. We assume that $V_{\ell}\in
C^{k}((\mathbb{R}^{2})^{\mathcal{R}_{\ell}}),k\geq 2$ and $V_{\ell}$ is
homogeneous outside the defect core $\Omega^{{\rm def}}$, namely, $V_{\ell}=V$
and $\mathcal{R}_{\ell}=\mathcal{R}$ for $\ell\in\Lambda\setminus\Omega^{{\rm
def}}$. A further condition of point symmetry $\mathcal{R}=-\mathcal{R}$, and
$V(\\{-g_{-\rho}\\}_{\rho\in\mathcal{R}})=V(g)$ is also assumed to keep the
homogeneity of the system outside the defect core [8].

We note that $V_{\ell}$ may vary according to $\ell$ inside each ${B_{R^{\rm
def}_{i}}}(x_{i})$ so that certain defects can be modeled. We also note that a
great number of practical interatomic potentials satisfy the assumptions of
smoothness and symmetry listed in the last paragraph including the widely used
embedded atom method (EAM) [38] and Finnis-Sinclair model [39] so that it is
not a serious restriction for realistic atomistic or material systems. We will
use the empirical EAM potential throughout our numerical experiments, see § 5
for detail.

Since the energy of an infinite lattice is typically ill-defined, with a
slight abuse of notation, we redefine the potential $V_{\ell}(y)$ as the
energy difference $V_{\ell}(y)-V_{\ell}(x)$ for point defects, which is
equivalent to assuming $V_{\ell}(x)=0$. As for anti-plane screw dislocation,
we denote $V_{\ell}(y):=V(y)-V(y_{0})$ (cf. [8, Section 3.1]). The energy
functional of the atomistic system

$\mathscr{E}^{\rm a}(y)=\sum_{\ell\in\Lambda}V_{\ell}(y)$ (2)

is then a meaningful object. By the assumptions of smoothness and symmetry on
$V_{\ell}$, $\mathscr{E}^{\rm a}(y)$ is well-defined for
$y-y^{B}\in\mathscr{U}^{1,2}$ with certain Fréchet differentiability, where
$y^{B}(x)=Bx+u_{0}(x)$. We refer to [8] for a detailed discussion.

The atomistic problem we would like to solve is to find a strongly stable
equilibrium $y$, such that, given a macroscopic applied strain
$B\in\mathbb{R}^{2\times 2}$,

$y\in\arg\min\big{\\{}\mathscr{E}^{\rm
a}(y)~{}\,\big{|}\,~{}y-y^{B}\in\mathscr{U}^{1,2}\big{\\}}.$ (3)

The strong stability on $y$ means that there exists a constant $c_{0}>0$ such
that

$\langle\delta^{2}\mathscr{E}^{\rm a}(y)v,v\rangle\geq c_{0}\|\nabla
v\|_{L^{2}}^{2},\quad\forall v\in\mathscr{U}^{1,2}.$

###  2.2 GRAC coupling method

The atomistic problem (3) is not computable in practice since it is defined on
an infinite domain and it considers every atom as a degree of freedom. The
former problem is easily circumvented by restricting the computational domain
to be $\Omega\subset\mathbb{R}^{2}$ which is a simply connected polygonal and
we assume that $\exists R\in\mathbb{R}$ such that $B_{R}\subset\Omega\subset
B_{c_{0}R}$ for some $c_{0}>1$ where $R$ is the radius of $\Omega$.

To deal with the later issue, namely the overwhelming number of degrees of
freedom, we first define a continuum approximation by the Cauchy-Born rule
which is both compatible with (2) and related to the strain. The Cauchy-Born
energy density functional $W:\mathbb{R}^{2\times 2}\to\mathbb{R}$ is simply
given by

$W({\sf F}):=\det(\mathsf{A}^{-1})\cdot V({\sf F}\mathcal{R}),$

where $V$ is the homogeneous site energy potential on $\Lambda^{\rm hom}$ [40,
41]. An atomistic-to-continuum (a/c) coupling method then essentially couples
the atomistic model with its Cauchy-Born approximation together. To formulate
our coupling scheme, we first decompose the computational domain into two
parts: $\Omega=\Omega^{\rm a}\bigcup\Omega^{\rm c}$ where the defect core
$\Omega^{{\rm def}}$ is contained in $\Omega^{\rm a}$. Next, we decompose the
set of atoms inside the atomistic region into a core atomistic set
$\Lambda^{\rm a}$ and an interface set $\Lambda^{\rm i}$ (usually a few
“layers” of atoms surrounding $\Lambda^{\rm a}$ for each $B_{R^{\rm
def}_{i}}(x_{i})$) such that $\Lambda^{{\rm a},{\rm
i}}:=\Lambda\bigcap\Omega^{\rm a}=\Lambda^{\rm a}\bigcup\Lambda^{\rm i}$ and
$\Lambda\cap\Omega^{{\rm def}}$ is contained in $\Lambda^{\rm a}$.

To reduce the number of degrees of freedom, we introduce the finite element
partition of the computational domain $\mathcal{T}_{h}$ so that
$\mathcal{T}_{h}=\mathcal{T}^{{\rm c}}_{h}\bigcup\mathcal{T}^{{\rm a}}_{h}$
where $\mathcal{T}^{\rm a}_{h}$ is the canonical triangulation induced by
$\Lambda^{{\rm a},{\rm i}}$ and $\mathcal{T}^{{\rm c}}_{h}$ be a shape-regular
simplicial partition (triangles for 2D) of the continuum region $\Omega^{\rm
c}$. We note that $\mathcal{T}^{\rm a}_{h}$ may contain ”holes” due to the
existence of defects (see Figure 6 for an illustration). We define
$\mathcal{N}_{h}$ be the set of nodes in $\mathcal{T}_{h}$ and
$\mathcal{F}_{h}^{{\rm c}}$ be the set of edges in $\mathcal{T}_{h}^{{\rm
c}}$.

Assuming $P_{1}$ finite element are used on $\mathcal{T}_{h}$, the general
formulation of an energy-based atomistic-to-continuum coupling energy
functional can be written as

$\displaystyle\mathscr{E}^{\rm ac}(y_{h}):=$
$\displaystyle\sum_{\ell\in\Lambda^{\rm
a}}V_{\ell}(y_{h})+\sum_{\ell\in\Lambda^{\rm i}}\omega_{\ell}V^{\rm
i}_{\ell}(y_{h})+\sum_{T\in\mathcal{T}_{h}}\omega_{T}W(\nabla y_{h}|_{T}),$
(4)

where $V_{\ell}^{\rm i}$ is a modified interface site potential and
$\omega_{\ell}$ and $\omega_{T}$ are certain coefficients which are usually
called effective volumes of lattice sites or elements (see [21, Section 2.2]
for a detailed discussion). We note that $\omega_{T}=0$ for
$T\in\mathcal{T}_{h}^{{\rm a}}\setminus\mathcal{T}_{h}^{\rm i}$ with
$\mathcal{T}_{h}^{\rm i}:=\\{T\in\mathcal{T}_{h}:\Lambda^{\rm i}\cap
T\neq\emptyset\\}$.

The interface potential in (4) is the key for constructing a consistent a/c
coupling method. In the current work, we restrict ourselves to the geometric
reconstruction based consistent a/c (GRAC) coupling method which is first
proposed in [21] and then further developed in [35]. The essential idea for
GRAC method is to introduce the geometric reconstruction parameters
$C_{\ell;\rho,\varsigma}$ so that for each $\ell\in\Lambda^{\rm
i},\rho,\varsigma\in\mathcal{R}_{\ell}$ the interface potential is redefined
by

$V_{\ell}^{\rm
i}(y):=V\Big{(}\big{(}{\textstyle\sum_{\varsigma\in\mathcal{R}_{\ell}}C_{\ell;\rho,\varsigma}D_{\varsigma}y(\ell)}\big{)}_{\rho\in\mathcal{R}_{\ell}}\Big{)}.$
(5)

The geometric reconstruction parameters are then determined by solving the
following equations:

$V_{\ell}^{\rm i}(y^{{\sf F}})=V(y^{{\sf
F}})\qquad\quad\forall\ell\in\Lambda^{\rm i},\quad{\sf
F}\in\mathbb{R}^{d\times d},$ (6)

and

$\langle\delta\mathscr{E}^{\rm ac}(y^{{\sf F}}),v\rangle=0\qquad\forall
v\in\mathscr{U}_{h},\quad{\sf F}\in\mathbb{R}^{d\times d},$ (7)

which are termed as the energy patch test consistency and the force patch test
consistency [42]. According to [21, Proposition 3.7], under the setting of the
current work, for $\ell\in\Lambda^{{\rm i}}$, the parameters
$C_{\ell;\rho,\zeta}$ in (5) are determined by: (a)
$C_{\ell,\rho,\varsigma}=0$ for $|(\rho-\varsigma)~{}\textrm{mod}~{}6|>1$; (b)
$C_{\ell,\rho,\rho-1}=C_{\ell,\rho,\rho+1}=1-C_{\ell,\rho,\rho}$; (c)
$C_{\ell,\rho,\rho}=C_{\ell+\rho,-\rho,-\rho}=1$ for
$\ell+\rho\in\Lambda^{{\rm i}}$; (d) $C_{\ell,\rho,\rho}=1$ for
$\ell+\rho\in\Lambda^{{\rm a}}$; (e) $C_{\ell,\rho,\rho}=2/3$ for
$\ell+\rho\in\Lambda^{{\rm c}}$. The coefficient $2/3$ given in (e) introduces
the name of the specific GRAC method we consider in the rest of this paper,
which is the GRAC23 method. However, we need to note that these coefficients
are not unique and how they are optimally determined are discussed in depth in
[35].

Let $\Omega_{h}=\bigcup_{T\in\mathcal{T}_{h}}T$. We define the space of
coarse-grained displacements to be

$\displaystyle\mathscr{U}_{h}:=\big{\\{}u_{h}:\Omega_{h}\to\mathbb{R}^{2}\,\big{|}\,$
$u_{h}$ is continuous and p.w. affine w.r.t. $\mathcal{T}_{h}$,
$\displaystyle\text{ $u_{h}=0$ on $\partial\Omega_{h}$ }\big{\\}}.$

The a/c coupling problem we would like to solve is to find

$y_{h}\in\arg\min\big{\\{}\mathscr{E}^{\rm
ac}(y_{h})~{}\,\big{|}\,~{}y_{h}-y^{B}\in\mathscr{U}_{h}\big{\\}}.$ (8)

####  2.2.1 The formulation of the stress tensors

The stress tensor, which is in principle defined by the first variation of the
energy, plays the key role in both a priori and a posteriori error analysis of
the a/c coupling methods. For any $y-y^{B}\in\mathscr{U}^{1,2}$, and
$y_{h}-y^{B}\in\mathscr{U}_{h}$, we define the atomistic stress tensor
$\sigma^{\rm a}(y;\cdot)\in{\rm P}_{0}(\mathcal{T}_{\rm a})^{2\times 2}$, the
continuum stress tensor $\sigma^{\rm c}(y_{h};\cdot)\in{\rm
P}_{0}(\mathcal{T}_{h})^{2\times 2}$ and the a/c stress tensor $\sigma^{{\rm
ac}}(y_{h};\cdot)\in{\rm P}_{0}(\mathcal{T}_{h})^{2\times 2}$ respectively by
the following identities

$\displaystyle\langle\delta\mathscr{E}^{\rm a}(y),v\rangle$
$\displaystyle=\sum_{T\in\mathcal{T}_{\rm a}}|T|\sigma^{\rm
a}(y;T)\nabla_{T}v,\qquad\forall v\in\mathscr{U}^{1,2},$ (9)
$\displaystyle\langle\delta\mathscr{E}^{\rm c}(y_{h}),v_{h}\rangle$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}|T|\sigma^{\rm
c}(y_{h};T)\nabla_{T}v_{h},\qquad\forall v_{h}\in\mathscr{U}_{h},$ (10)
$\displaystyle\langle\delta\mathscr{E}^{\rm ac}(y_{h}),v_{h}\rangle$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}|T|\sigma^{{\rm
ac}}(y_{h};T)\nabla_{T}v_{h},\qquad\forall v_{h}\in\mathscr{U}_{h},$ (11)

where $\mathcal{T}_{\rm a}$ is the canonical triangulation induced by the
reference lattice $\Lambda$ and $\mathscr{E}^{\rm
c}(y_{h}):=\sum\limits_{T\in\mathcal{T}_{h}}|T|W(\nabla_{T}y_{h})$. Since we
only consider nearest neighbor interactions, we can choose the following
stress tensors from the first variations (9)-(11),

$\displaystyle\sigma^{\rm a}(y;T):=~{}$
$\displaystyle\frac{1}{\det\mathsf{A}}\sum_{b=(\ell,\ell+\rho)\in\partial
T}\partial_{\rho}V_{\ell}\otimes a_{\rho},$ (12) $\displaystyle\sigma^{\rm
c}(y_{h};T):=~{}$ $\displaystyle\partial
W(\nabla_{T}y_{h})=\frac{1}{\det\mathsf{A}}\sum_{j=1}^{6}\partial_{j}V(\nabla_{T}y_{h})\otimes
a_{j},$ (13) $\displaystyle\sigma^{{\rm
ac}}(y_{h};T):=~{}\frac{1}{\det\mathsf{A}}$
$\displaystyle\sum_{b=(\ell,\ell+\rho)\in\partial
T}\partial_{\rho}V_{\ell}(I_{\rm a}y_{h})\otimes
a_{\rho}+\omega_{T}\sigma^{\rm c}(y_{h};T).$ (14)

We note that the stress tensors defined in (9)-(11) are not unique. Equation
(9)-(11) may hold with any of the stress tensors differ from (12)-(14) with a
piecewise constant tensor field $\sigma\in{\rm P}_{0}(\mathcal{T})^{2\times
2}$ that satisfies the divergence free condition

$\sum_{T\in\mathcal{T}}|T|\sigma(T):\nabla_{T}v\equiv 0,\qquad\forall
v\in({\rm P}_{1}(\mathcal{T}))^{2}.$ (15)

We postpone the detail of the construction and optimization of such divergence
free field to § 4.1 where a local stress tensor correction is applied to the
coupling stress tensor $\sigma^{\rm ac}$.

##  3 The Analysis of the Residual Based A Posteriori Error Estimates

In this section, we give the a posteriori error estimates for the GRAC method
introduced in § 2.2. We first present some useful a priori results in § 3.1
that motivates the separation of the a posteriori error estimator. We then
give a brief review of the original residual based a posteriori error
estimates in § 3.1 which is derived and analyzed in [33]. The error estimates
provide a solid foundation of the current work but the problem of high
computational cost arise in computing the modeling part of the a posteriori
error estimator. We finally concentrate on the analysis and numerical
verification of the relation between the modeling and the coarsening parts of
the estimator whose quotient in fact possesses a certain ratio under suitable
assumptions. Such ratio relation serves as the theoretical basis for deriving
modified a posteriori error estimators with much better computational
efficiency that are discussed in § 4 and § 5. We note that the content in §
3.1 and § 3.2 can be found in detail in earlier works [33, 34] and we include
them only for completeness and refer to related literature where necessary.

###  3.1 A Priori Error Estimates

If $y$ is a strongly stable solution of (3) and is sufficiently smooth in the
continuum region, there exists a strongly stable solution $y_{h}$ to (8) and a
constant $C^{a-priori}$ for GRAC method such that [8, 3],

$\displaystyle\|\nabla u_{h}-\nabla u\|_{L^{2}}\leq
C^{a-priori}\big{(}\|D^{2}u\|_{\ell^{2}(\mathcal{I}_{\rm
ext})}+\|D^{3}u\|_{\ell^{2}(\Lambda\cap\Omega^{\rm
c})}+\|D^{2}u\|^{2}_{\ell^{4}(\Lambda\cap\Omega^{\rm c})}$
$\displaystyle+\|hD^{2}u\|_{\ell^{2}(\Lambda\cap\Omega^{\rm
c})}+\|Du\|_{\ell^{2}(\Lambda\setminus B_{R/2})}\big{)},$ (16)

where $u=y-y^{B},u_{h}=y_{h}-y^{B}$ and $\mathcal{I}_{\rm ext}$ is the
extended interface region whose definition is given in [21, § 5]. The right
hand side of (3.1) can be essentially split into three parts and we give a
discussion here as it helps to understand the correspondence between the a
priori error estimates and the a posteriori error estimates.

The first part is the modeling error defined by $\zeta_{\rm
mo}:=\|D^{2}u\|_{\ell^{2}(\mathcal{I}_{\rm
ext})}+\|D^{3}u\|_{\ell^{2}(\Lambda\cap\Omega^{\rm
c})}+\|D^{2}u\|^{2}_{\ell^{4}(\Lambda\cap\Omega^{\rm c})}$. The error is due
to the difference between the atomistic model and the continuum model which
lies in both the discrepancy of the nonlocality and the locality of the
atomistic and the coupling energy and the lattice and mesh on which we solve
the problems. This is the fundamentally different from the discretization of a
PDE using certain numerical schemes, see [21, 33] for further discussion.

The second part is the coarsening error denoted by $\zeta_{\rm
cg}:=\|hD^{2}u\|_{\ell^{2}(\Lambda\cap\Omega^{\rm c})}$. This is caused by the
application of the finite element method in the continuum region and is a
reminiscence of the discretization error of the numerical solution to an
elliptic differential equation.

The third part is the truncation error denoted by $\zeta_{\rm
tr}:=\|Du\|_{\ell^{2}(\Lambda\setminus B_{R/2})}$, which arises from the
finite size of our computational domain. We note that the region $B_{R/2}$
comes from the construction of the truncation operator, see [8, Lemma 7.3] for
a better detail.

###  3.2 Residual based a posterior error estimates

We derive the a posteriori error estimates based on the residual for the a/c
method in this section. To avoid the unnecessary complexity of the
presentation while keeping the main ideas and steps, we restrict our
discussion to the case of the GRAC23 method with nearest neighbor multibody
interactions. We refer to [34] for the extension to the finite range
interaction. We also note the derivation is concise as the detail can be found
in [33].

The first variation of the atomistic variational problem (3) and that of the
a/c coupling problem (8) are to find $y-y^{B}\in\mathscr{U}^{1,2}$ and
$y_{h}-y^{B}\in\mathscr{U}_{h}$ such that

$\langle\delta\mathscr{E}^{\rm a}(y),v\rangle=0,\quad\forall
v\in\mathscr{U}^{1,2},\text{ and }\langle\delta\mathscr{E}^{\rm
ac}(y_{h}),v_{h}\rangle=0,\quad\forall v_{h}\in\mathscr{U}_{h},$ (17)

respectively. Since we limit ourselves in a finite computational domain
$\Omega$, we choose the same truncation operator
$T_{R}:\mathscr{U}^{1,2}\to\mathscr{U}_{R}$ where
$\mathscr{U}_{R}:=\\{u\in\mathscr{U}^{1,2}~{}|~{}u(x)=0,~{}\forall
x\in\Lambda\backslash\Omega\\}$.

The residual operator ${\sf R}$ is then defined as an operator on
$\mathscr{U}^{1,2}$ and ${\sf R}[v]$ is given by

${\sf R}[v]=\langle\delta\mathscr{E}^{\rm a}(I_{\rm
a}y_{h}),v\rangle,\quad\forall v\in\mathscr{U}^{1,2}.$ (18)

Let $v_{R}=T_{R}v$, and
$v_{h}=\mathcal{C}_{h}T_{R}v:\mathscr{U}^{1,2}\to\mathscr{U}_{h}$, where
$\mathcal{C}_{h}:\mathscr{U}_{R}\to\mathscr{U}_{h}$ is the modified Clément
operator [33, 43, 44]. By the second equation in (17), we separate the
residual ${\sf R}[v]$ into three groups

$\displaystyle{\sf R}[v]=\langle\delta\mathscr{E}^{\rm a}(I_{\rm
a}y_{h}),v\rangle=$ $\displaystyle\langle\delta\mathscr{E}^{\rm a}(I_{\rm
a}y_{h}),v\rangle-\langle\delta\mathscr{E}^{\rm ac}(y_{h}),v_{h}\rangle$
$\displaystyle=$ $\displaystyle\big{\\{}\langle\delta\mathscr{E}^{\rm
a}(I_{\rm a}y_{h}),v_{R}\rangle-[\delta\mathscr{E}^{\rm
ac}(y_{h}),v_{R}]\big{\\}}+\big{\\{}[\delta\mathscr{E}^{\rm
ac}(y_{h}),v_{R}]-\langle\delta\mathscr{E}^{\rm
ac}(y_{h}),v_{h}\rangle\big{\\}}$
$\displaystyle+\big{\\{}\langle\delta\mathscr{E}^{\rm a}(I_{\rm
a}y_{h}),v\rangle-\langle\delta\mathscr{E}^{\rm a}(I_{\rm
a}y_{h}),v_{R}\rangle\big{\\}}=:{\sf R}_{\rm mo}+{\sf R}_{\rm cg}+{\sf R}_{\rm
tr},$ (19)

where the operator $[\cdot,\cdot]$ is defined by

$\displaystyle[\delta\mathscr{E}^{\rm
ac}(y_{h}),v_{R}]:=\sum_{T\in\mathcal{T}_{h}}\int_{T}\sigma^{{\rm
ac}}(y_{h},T)\nabla v_{R}\,{\rm d}x=$
$\displaystyle\sum_{T\in\mathcal{T}_{h}}\sigma^{{\rm
ac}}(y_{h},T)(\sum_{T^{\prime}\in\mathcal{T}_{\rm a},T^{\prime}\bigcap
T\neq\emptyset}|T\bigcap T^{\prime}|\nabla v_{R})$ $\displaystyle=$
$\displaystyle\sum_{T\in\mathcal{T}_{\rm
a}}|T|\big{(}\sum_{T^{\prime}\in\mathcal{T}_{h},T^{\prime}\bigcap
T\neq\emptyset}\frac{|T^{\prime}\bigcap T|}{|T|}\sigma^{{\rm
ac}}(y_{h},T^{\prime})\big{)}\nabla v_{R},$ (20)

which is different from $\langle\cdot,\cdot\rangle$ since
$v_{R}\notin\mathscr{U}_{h}$. We essentially separate the residual into the
modeling part ${\sf R}_{\rm mo}$, the the coarsening part ${\sf R}_{\rm cg}$
and the truncation part ${\sf R}_{\sf tr}$.

The estimate of the truncation residual ${\sf R}_{\rm tr}$ is the easiest
among the three and is simply given by

$|{\sf R}_{\rm tr}|\leq C^{{\rm tr}}\|\sigma^{\rm a}(I_{\rm
a}u_{h})-\sigma^{0}\|_{L^{2}(\Omega\setminus B_{R/2})}\|\nabla
v\|_{L^{2}}=:\eta_{\rm tr}(y_{h})\|\nabla v\|_{L^{2}},$ (21)

where $C^{{\rm tr}}$ is the constant appeared when estimating the difference
of $Dv_{R}$ and $Dv$ in terms of $Dv$, and we can specify the divergence-free
stress tensor $\sigma^{0}=\partial W(\nabla y^{B})$ (cf. [33, §3.1.1]).

The coarsening residual ${\sf R}_{\rm cg}$ is estimated by

$|{\sf R}_{\rm cg}|\leq\sqrt{3}C^{\rm
tr}C_{\mathcal{T}_{h}}(\sum_{f\in\mathcal{F}_{h}^{\rm c}\bigcap{\rm
int}(\Omega)}(h_{f}\llbracket\sigma^{{\rm
ac}}\rrbracket_{f})^{2})^{\frac{1}{2}}\|\nabla v\|_{L^{2}(\Omega)}=:\eta_{\rm
cg}(y_{h})\|\nabla v\|_{L^{2}(\Omega)},$ (22)

where $\mathcal{F}_{h}^{\rm c}\bigcap{\rm int}(\Omega)$ is the set of interior
edges, $\llbracket\sigma^{{\rm ac}}\rrbracket_{f}:=\sigma^{{\rm
ac}}(y_{h},T^{+}_{f})\nu^{+}+\sigma^{{\rm ac}}(y_{h},T^{-}_{f})\nu^{-}$
denotes the jump of $\sigma^{{\rm ac}}$ across each interior edge $f$ and
$C_{\mathcal{T}_{h}}$ is the constant for the modified Clément operator that
only depends on the shape regularity of the mesh. We refer to [33, §3.1.3] for
a thorough discussion.

For the analysis of the modeling error ${\sf R}_{\rm mo}$, the formulation of
$\langle\delta\mathscr{E}^{\rm a}(I_{\rm a}y_{h}),v_{R}\rangle$ and
$[\delta\mathscr{E}^{\rm ac}(y_{h}),v_{R}]$ given in (11) and (20) lead to
[33, §3.1.2],

$|{\sf R}_{\rm mo}|\leq C^{{\rm
tr}}\big{\\{}\sum_{T^{\prime}\in\mathcal{T}_{\rm a}}|T|\big{[}\sigma^{\rm
a}(I_{\rm a}y_{h},T^{\prime})-\sum_{T\in\mathcal{T}_{h},T\bigcap
T^{\prime}\neq\emptyset}\frac{|T\bigcap T^{\prime}|}{|T^{\prime}|}\sigma^{{\rm
ac}}(y_{h},T)\big{]}^{2}\big{\\}}^{\frac{1}{2}}\|\nabla v\|_{L^{2}}=:\eta_{\rm
mo}(y_{h})\|\nabla v\|_{L^{2}}.$ (23)

######  Remark 2.

We have to note that though the estimate of the modeling residual $\eta_{\rm
mo}$ is probably the most straight forward one among the three, it consumes
the most computational cost in the adaptive computation. It is easily observe
from (23) that to compute $\eta_{\rm mo}$ we need to find $|T^{\prime}\bigcap
T|$ for every $T\in\mathcal{T}_{\rm a}$ and $T^{\prime}\in\mathcal{T}_{h}$.
Thus the complex geometry, which is illustrated soon in Figure 2, enters into
the computation when we assign the global error estimator into local
contributions in each adaptive step. We give an estimate of the computational
cost of the modeling residual in § 3.3 and from then on, we will discuss how
to deal with this inefficiency that eventually reduces the computational cost
by one order with respect to the number of degrees of freedom.

Combining the residual estimates in (21), (22) and (23), the stability results
in [33, Section 3.2] and a quantitative version of the inverse function
theorem [3, Lemma B.1], we can obtain the following theorem for the rigorous
global a posteriori existence and error estimate [33, Theorem 3.7].

######  Theorem 1.

Let $y_{h}$ be the a/c solution defined in (17), $\eta(y_{h})=\eta_{{\rm
tr}}(y_{h})+\eta_{\rm mo}(y_{h})+\eta_{\rm cg}(y_{h})$ be the total residual
estimator, $\gamma(y_{h})$ be the stability constant defined in [33, § 3.2.4],
and $L_{1}$ be the Lipschitz constant of $\delta^{2}\mathscr{E}^{\rm a}$ whose
definition can be found in [3, Lemma 3.9]. Under the assumption that
$\gamma(y_{h})>0$ and $2L_{1}\eta(y_{h})<\gamma(y_{h})^{2}$, there exists a
unique $y$ satisfying $y-y^{B}\in\mathscr{U}^{1,2}$ which solves the atomistic
variational problem (17), and satisfies the following a posteriori error
bound,

$\|\nabla I_{\rm a}y_{h}-\nabla y\|_{L^{2}}\equiv\|\nabla I_{\rm
a}u_{h}-\nabla u\|_{L^{2}}\leq 2\frac{\eta(y_{h})}{\gamma(y_{h})}.$ (24)

###  3.3 Local contribution and analysis of the modeling and the coarsening
residual estimators

The a posteriori error estimates we just proposed is global. In order to
design and implement the adaptive algorithm, we need to assign the global
residual estimator $\eta(y_{h})$ into elementwise local contributions except
for the truncation part. The truncation residual estimator is often computed
directly in practice since it is relatively small if the computational domain
$\Omega$ is sufficiently large. We then check whether it exceeds the threshold
during the main adaptive process. See Algorithm 2 in § 4.2 for detail.

Following standard adaptive finite element procedure, the local coarsening
residual estimator for an element $T\in\mathcal{T}_{h}$ can be defined as

$\eta_{{\rm cg}}(y_{h};T):=\Bigg{(}\sum_{f\in\mathcal{F}_{h}\bigcap
T\in\mathcal{T}_{h}}\frac{1}{2}(h_{f}\llbracket\sigma^{{\rm
ac}}\rrbracket_{f})^{2}\Bigg{)}^{\frac{1}{2}}.$ (25)

The local modeling residual estimator for an element $T\in\mathcal{T}_{h}$ is
defined as the sum of the differences between stress tenors of $T$ and those
of $T^{\prime}\in\mathcal{T}_{\rm a}$ satisfying $T\bigcap
T^{\prime}\neq\emptyset$ so that

$\eta_{{\rm mo}}(y_{h};T)=\Bigg{(}\sum_{T^{\prime}\in\mathcal{T}_{\rm
a},T^{\prime}\bigcap T\neq\emptyset}\eta_{{\rm
mo}}(y_{h};T^{\prime},T)\Bigg{)}^{\frac{1}{2}},\quad\forall
T\in\mathcal{T}_{h},$ (26)

where

$\eta_{{\rm mo}}(y_{h};T^{\prime},T):=|T^{\prime}\bigcap T|\Bigg{(}\sigma^{\rm
a}(I_{\rm a}y_{h},T^{\prime})-\sum_{T\in\mathcal{T}_{h},T^{\prime}\bigcap
T\neq\emptyset}\frac{|T\bigcap T^{\prime}|}{|T^{\prime}|}(\sigma^{{\rm
ac}}(y_{h},T))\Bigg{)}^{2}.$

We note that the global modeling residual estimator $\eta_{{\rm mo}}(y_{h})$
defined in (23) is summed over $\mathcal{T}_{\rm a}$ for the simplicity of
presentation while the elementwise modeling residual estimator is defined on
each $T\in\mathcal{T}_{h}$ for the purpose of adaptive mesh refinements, and
the summation of $\eta_{{\rm mo}}(T)$ over $T\in\mathcal{T}_{h}$ and
$\eta_{{\rm mo}}(y_{h})$ are essentially equal up to a constant $C^{\rm tr}$.
We omit the constants $C^{\rm tr}$ and $C_{\mathcal{T}_{h}}$ here only for the
sake of presentation and will determine the values of these constants
empirically in our adaptive computations (cf. § 4.2).

As we have already pointed out in Remark 2, evaluating $\eta_{{\rm
mo}}(T^{\prime})$ requires computing $|T\cap T^{\prime}|$ for every
$T\in\mathcal{T}_{\rm a}$ and $T^{\prime}\in\mathcal{T}_{h}$ where finding the
geometric relationship between $T$ and $T^{\prime}$ is essential. See Figure 2
for an illustration of the mismatch of $T\in\mathcal{T}_{\rm a}$ and
$T^{\prime}\in\mathcal{T}_{h}$. We note that the computational cost of
evaluating $\eta_{{\rm mo}}$ is proportional to the number of atoms in the
continuum region since in our implementation we loop over all the elements in
$\mathcal{T}_{{\rm a}}$ lying in the continuum region to record the geometric
relationship. The number of atoms in the continuum region is about the square
of the number of elements in $\mathcal{T}_{h}$. Moreover, the ratio of the
number of atoms in the atomistic region and the number of elements in
$\mathcal{T}_{h}$ remains approximately a constant under certain assumptions
(see [3, Section 8.1] for detail). Hence, if $N$ is the total number of
degrees of freedom of the a/c coupling problem (8), then the computational
cost of evaluating $\eta_{{\rm mo}}(T^{\prime})$ can be roughly estimated by
$O(N^{2})$. In contrast, the cost of evaluating the coarsening residual
$\eta_{{\rm cg}}(T)$ is proportional to the number of edges
$f\in\mathscr{F}^{\rm c}_{h}$ of the elements in the continuum region, which
makes it only of $O(N)$. We give a numerical verification of the computational
cost in § 5 (cf. Figure 8a for micro-crack, Figure 10a for anti-plane screw
dislocation and Figure 13a for multiple vacancies).

Figure 2: Illustration of $|T^{\prime}\cap T|$. The edges of
$T\in\mathcal{T}_{h}$ are marked with long blue solid lines whereas the
elements $T^{\prime}\in\mathcal{T}_{\rm a}$ at the boundaries of the elements
$T\in\mathcal{T}_{h}$ are marked in grey.

To avoid the exceedingly high cost of computation, we apply an approximation
of the elementwise residual estimator $\eta_{{\rm mo}}(y_{h};T)$ in our
adaptive computations. The idea of the approximation comes from a simple
observation that the elementwise modeling residual estimator $\eta_{{\rm
mo}}(y_{h};T)$ should somehow be “negligible” compared with the elementwise
coarsening residual estimator $\eta_{{\rm cg}}(y_{h};T)$. Such observation may
be easily verified if we take the spring potential (essentially a linearized
model) as an example. With the assumption that $u$ is relative smooth in the
regions that are far-away from the defect core, we can directly estimate
$\eta_{{\rm cg}}(y_{h};T)\approx O(h_{T}^{2})$ and $\eta_{{\rm
mo}}(y_{h};T)\approx O(h_{T})$ by some algebraic manipulations according to
the definitions (25) and (26). Since $h_{T}\gg 1$ for any
$T\in\mathcal{T}_{h}$ that is sufficiently far away from the defect core, we
have a good reason to believe that $\eta_{{\rm mo}}(y_{h};T)$ should be a
quantity of higher order compared with $\eta_{{\rm cg}}(y_{h};T)$. For the
rest of this section, we make this argument rigorous and numerically verify
its validity.

We start with the a priori error estimates (3.1). We can further split the
modeling error by

$\displaystyle\zeta_{{\rm mo}}=\|D^{2}u\|_{\ell^{2}(\mathcal{I}^{\rm
ext})}+\|D^{3}u\|_{\ell^{2}(\Lambda\cap\Omega^{{\rm
c}})}+\|D^{2}u\|^{2}_{\ell^{4}(\Lambda\cap\Omega^{{\rm c}})}=:\zeta_{\rm
mo}^{\rm i}(y)+\zeta_{\rm mo}^{\rm c}(y),$ (27)

where $\zeta_{\rm mo}^{\rm i}(y)$ and $\zeta_{\rm mo}^{\rm c}(y)$ are
essentially the modeling error at the interface and in the continuum region.
If we assume that the computational domain contains only a single atomistic
region $\Omega^{{\rm a}}$ whose radius is denoted by $R_{{\rm a}}$, then the
following theorem states that for $R_{\rm a}$ sufficiently large, $\zeta_{\rm
mo}^{\rm c}(y)$ is negligible compared with $\zeta_{\rm cg}(y)$ and their
elementwise counterparts obey a certain ratio relation.

######  Theorem 2.

Suppose Assumption 1 holds and the computational domain only contains a single
atomistic region. Let $y$ be a strongly stable solution of (3), $u=y-y^{B}$,
and $\zeta_{\rm mo}^{\rm c}(y)$ and $\zeta_{\rm cg}(y)$ be defined in (27). If
the triangulation $\mathcal{T}_{h}$ is a graded mesh such that the mesh size
function $h(x)$ for $x\in T\in\mathcal{T}_{h}$ satisfies

$|h(x)|\leq C^{\rm{mesh}}\big{(}\frac{|x|}{R_{\rm
a}}\big{)}^{\beta},\quad\text{ with }1<\beta<\frac{d+2}{2},$ (28)

we then have

$\displaystyle\lim_{R^{{\rm a}}\rightarrow\infty}\frac{\zeta^{\rm c}_{{\rm
mo}}(y)}{\zeta_{\rm cg}(y)}=0.$ (29)

Moreover, for $T\in\mathcal{T}_{h}$, if we define the local error
contributions by $\zeta_{\rm cg}(y;T):=\|hD^{2}u\|_{\ell^{2}(\Lambda_{T})}$
and $\zeta^{\rm c}_{{\rm
mo}}(y;T):=\|D^{3}u\|_{\ell^{2}(\Lambda_{T})}+\|D^{2}u\|^{2}_{\ell^{4}(\Lambda_{T})}$
respectively with $\Lambda_{T}:=\Lambda\cap T$. Then for $R_{{\rm a}}$
sufficiently large, there exits a constant $C^{a-priori}_{{\rm mo}-{\rm
cg}}>0$ such that

$\frac{\zeta^{\rm c}_{{\rm mo}}(y;T)}{\zeta_{\rm
cg}(y;T)}\leq\frac{C^{a-priori}_{{\rm mo}-{\rm cg}}}{h_{T}},$ (30)

where $h_{T}:=\mathrm{diam}(T)$ for $T\in\mathcal{T}_{h}$.

###### Proof.

Let $y$ be a strongly stable solution of (3) and $u=y-y^{B}$. We have the
following sharp decay estimates of $u$ for both point defects and dislocations
[8, Theorem 2.3, 3.6], for $|\ell|$ sufficiently large and $j=1,2,3$,

$\displaystyle\big{|}D^{j}u(\ell)\big{|}$ $\displaystyle\leq
C\big{(}1+|\ell|\big{)}^{-d-j+1},\qquad\quad\textrm{point defects},$
$\displaystyle\big{|}D^{j}u(\ell)\big{|}$ $\displaystyle\leq
C\big{(}1+|\ell|\big{)}^{-j-1}\log(|\ell|),\qquad\textrm{dislocations},$ (31)

where the stencil norm is defined by
$|Du(\ell)|:=\big{(}\sum_{\rho\in\mathcal{R}_{\ell}}D_{\rho}u(\ell)\big{)}^{1/2}$.
For ${\bm{\rho}}\in(\mathcal{R}_{\ell})^{j}$,
$D_{\bm{\rho}}u:=D_{\rho_{1}}\cdots D_{\rho_{j}}u$ denotes a $j$-th order
derivative, and $D^{j}u$ is defined recursively by $D^{j}u:=DD^{j-1}u$ which
denotes the $j$-th order collection of derivatives.

By the definition of $\zeta_{\rm mo}^{{\rm c}}(y)$, $\zeta_{\rm cg}(y)$,
together with the sharp decay estimates (3.3) and the restriction of the mesh
size function $h(x)$ by (28), for $R_{\rm a}$ sufficiently large, which is
equivalent to $|\ell|$ sufficiently large, we obtain for point defects,

$\displaystyle\zeta_{\rm mo}^{{\rm c}}(y)$ $\displaystyle\leq C_{\zeta_{\rm
mo}}\Big{(}\int_{R_{{\rm a}}}^{\infty}r^{d-1}(r^{-d-2})^{2}\,{\rm
d}r\Big{)}^{1/2}\leq C_{\zeta_{\rm mo}}(R_{{\rm a}})^{-d/2-2},$
$\displaystyle\zeta_{\rm cg}(y)$ $\displaystyle\leq C_{\zeta_{\rm
cg}}\Big{(}\int_{R_{{\rm a}}}^{\infty}r^{d-1}\Big{(}\frac{r}{R_{{\rm
a}}}\Big{)}^{2\beta}(r^{-d-1})^{2}\,{\rm d}r\Big{)}^{1/2}\leq C_{\zeta_{\rm
cg}}(R_{{\rm a}})^{-d/2-1}.$ (32)

For dislocations, similarly we have the estimates [8, Proposition 2]

$\displaystyle\zeta_{\rm mo}^{{\rm c}}(y)\leq C^{\prime}_{\zeta_{\rm
mo}}(R_{{\rm a}})^{-2}\quad\text{and}\quad\zeta_{\rm cg}(y)\leq
C^{\prime}_{\zeta_{\rm cg}}(R_{{\rm a}})^{-1}.$ (33)

Hence, we have

$\lim_{R^{{\rm a}}\rightarrow\infty}\frac{\zeta^{{\rm c}}_{{\rm
mo}}(y)}{\zeta_{\rm cg}(y)}=\lim_{R^{{\rm
a}}\rightarrow\infty}C^{a-priori}_{{\rm mo}-{\rm cg}}R_{{\rm a}}^{-1}=0,$

with $C^{a-priori}_{{\rm mo}-{\rm cg}}:=\max\\{C_{\zeta_{\rm
mo}}/C_{\zeta_{\rm cg}},C^{\prime}_{\zeta_{\rm mo}}/C^{\prime}_{\zeta_{\rm
cg}}\\}$, which yields the stated result (29) for both point defects and
dislocations.

Let $\mathcal{T}_{{\rm a}}$ be the canonical triangulation induced by
$\Lambda$. Given $T\in\mathcal{T}^{\rm c}_{h}$, we define the set of “inner
finest elements” by
$\mathcal{T}_{\varepsilon}^{T}:=\\{T_{\varepsilon}\in\mathcal{T}_{{\rm
a}}~{}\big{|}~{}|T_{\varepsilon}\cap T|=|T_{\varepsilon}|\\}$.

To obtain the estimate (30), we first have for any $T\in\mathcal{T}^{\rm
c}_{h}$ and $T_{\varepsilon}\in\mathcal{T}^{T}_{\varepsilon}$

$\frac{\zeta^{{\rm c}}_{{\rm mo}}(y;T)}{\zeta_{\rm cg}(y;T)}\leq
C\frac{\|D^{3}u\|_{\ell^{2}(\Lambda_{T})}}{\|hD^{2}u\|_{\ell^{2}(\Lambda_{T})}}\leq
C^{a-priori}_{{\rm mo}-{\rm
cg}}h^{-1}_{T}\frac{\|D^{3}u\|_{\ell^{2}(\Lambda_{T_{\varepsilon}})}}{\|D^{2}u\|_{\ell^{2}(\Lambda_{T_{\varepsilon}})}}.$
(34)

Denoting the barycenter of $T_{\varepsilon}$ by $x_{T_{\varepsilon}}$ and
applying the decay estimates (3.3) again, we arrive at

$\displaystyle\frac{\zeta^{{\rm c}}_{{\rm mo}}(y;T)}{\zeta_{\rm cg}(y;T)}\leq
C^{a-priori}_{{\rm mo}-{\rm
cg}}h^{-1}_{T}\frac{(x_{T_{\varepsilon}}+1)^{-3}-x_{T_{\varepsilon}}^{-3}}{(x_{T_{\varepsilon}}+1)^{-2}-x_{T_{\varepsilon}}^{-2}}\leq\frac{C^{a-priori}_{{\rm
mo}-{\rm cg}}}{h_{T}},$ (35)

where in the last inequality we use the fact $x_{T_{\varepsilon}}>1$.
Combining (34) and (35) we achieve the estimate (30) which finishes the proof
of Theorem 2. ∎

We then give a similar a posteriori estimate under suitable assumptions which
is computable and can be employed in the adaptive computations. We note that
similar to $\zeta^{\rm c}_{{\rm mo}}(y;T)$, $\eta^{\rm c}_{\rm mo}(y_{h};T)$
represents the elementwise residual in the continuum region.

######  Corollary 3.

Under the assumptions of Theorem 2, if $y_{h}$ is a strongly stable solution
of (8) satisfying (3.1) and we additionally assume that $u_{h}$ satisfies the
decay estimates (3.3), then for $T\in\mathcal{T}_{h}$, there exists constants
$C^{\rm a-post}$ and $C^{\rm a-post}_{\rm mo-cg}$ such that

$\frac{\eta^{\rm c}_{\rm mo}(y_{h};T)}{\eta_{\rm cg}(y_{h};T)}\leq C^{\rm
a-post}\frac{\zeta^{\rm c}_{\rm mo}(y;T)}{\zeta_{\rm cg}(y;T)}\leq\frac{C^{\rm
a-post}_{\rm mo-cg}}{h_{T}}.$ (36)

###### Proof.

We sketch out the proof here for the sake of simplicity. If $u$ and $u_{h}$
are the strongly stable solutions of (3) and (8) respectively satisfying
(3.1), for $R_{\rm a}$ sufficiently large and $T$ lying in the region that far
away from the defect core, we can roughly (abuse of notations) estimate that
$\sigma^{\rm a}(u;T)/\sigma^{\rm ac}(u;T)\lesssim\sigma^{\rm
a}(u_{h};T)/\sigma^{\rm ac}(u_{h};T)$. Then, following the a priori
consistency estimates presented in [21, Section 5] and the a posterior
residual estimates given in Theorem 1, we have an asymptotic result that
$\eta^{{\rm c}}_{{\rm mo}}(y_{h};T)/\eta_{\rm cg}(y_{h};T)\leq C^{\rm
a-post}\zeta^{{\rm c}}_{{\rm mo}}(y;T)/\zeta_{\rm cg}(y;T)$ under an
additional assumption of $u_{h}$ satisfying the same decay estimates as $u$.
According to the estimate of $\zeta^{{\rm c}}_{{\rm mo}}(y;T)/\zeta_{\rm
cg}(y;T)$ in Theorem 2, we obtain the stated result by choosing $C^{\rm
a-post}_{\rm mo-cg}=C^{a-priori}_{{\rm mo}-{\rm cg}}C^{\rm a-post}$, where
$C^{a-priori}_{{\rm mo}-{\rm cg}}$ is introduced in (30). ∎

We note that there are three key assumptions for the analytical result (36) in
Corollary 3 to hold. The first assumption is that $u_{h}$ shares the same
decay estimates as $u$ which results in the first inequality of (36). Such
assumption is not easy to prove rigorously in principle and we refer the
interested reader to the analysis in [45, § 5.1] in the context of QM/MM
coupling methods. The second and the third assumptions are that only a single
atomistic region $\Omega^{{\rm a}}$ exists and $R_{{\rm a}}$ is sufficiently
large which result in the a priori decay estimate (3.3) and consequently the
second inequality of (36) (or essentially (30)). The second assumption is
certainly violated by the defect of multiple vacancies in the first steps of
the adaptive computation where multiple disjoint atomistic regions exist (see
the first four subplots in Figure 11 in § 5.3). The third assumption may not
be satisfied by all three cases of defects in the first steps of the adaptive
computations and the quantitative characterization of how large $R_{\rm a}$
should be in order to guarantee the decay estimates may involve substantial
additional technicalities.

Therefore, instead of pursuing the theoretical consideration further, we
verify the first assumption and the estimate (36) numerically in Figure 3 and
Figure 4 respectively whose further detail can be found later in § 5.

(a) Micro-crack

(b) Anti-plane screw dislocation

(c) Multiple vacancies

Figure 3: Numerical verification of the decay rate of $u$ and $u_{h}$.

The blue dots in each subplot in Figure 3 are the modules of the atomistic
solution $u$ sampled in six different directions in a spokewise form from the
defect core $\Omega^{\rm def}$ while the red asterisks are taken from the a/c
solution $u_{h}$ in one of the adaptive steps when $R_{{\rm a}}\geq 10$ for
the first time for micro-crack and anti-plane screw dislocation and when
$R_{{\rm a}}\geq 8$ (or more precisely, the radius of each atomistic region is
greater than or equal to $8$) for the first time for multiple vacancies while
the atomistic regions are still disjoint. The $x$-axis denotes the distance
from the “center” of the defects. We note that we define the “center” of the
defect core of the separated vacancies as their geometric center and an
increase of $|u|$ and $|u_{h}|$ appears in the first half of the $x-$axis in
the figure since it gets closer to the true defects as the distance $|x|$
grows initially. We observe in all three cases that $u_{h}$ essentially
follows the same decay (or increase) rates as $u$ just as we assumed in
Corollary 3. In addition, we observe a close match between $u$ and $u_{h}$ now
but the detailed convergence of the error together with a thorough discussion
is presented in § 5.

(a) Micro-crack

(b) Anti-plane screw dislocation

(c) Multiple vacancies

Figure 4: Elementwise modeling error over elementwise coarsening error versus
the diameter of corresponding element.

Figure 4 directly verifies the estimate in (36) in Corollary 3. The three
subplots are the snapshots at the same mesh as those in Figure 3. We clearly
see that with the increasing diameters of the elements in the continuum
region, the ratio of the elementwise modeling error over the elementwise
coarsening error, namely $\eta^{{\rm c}}_{{\rm mo}}(y_{h};T)/\eta_{\rm
cg}(y_{h};T)$, decays asymptotically with the rate $h^{-1}_{T}$. Though the
exact value of $C^{\rm a-post}_{\rm mo-cg}$ is difficult to obtain, we use the
largest mo-cg ratio for the elements $T$ near $\Omega^{\rm a}\cup\Omega^{\rm
buf}$ as the estimate of $C^{\rm a-post}_{\rm mo-cg}$ in practice, where
$\Omega^{\rm buf}$ is the buffer region whose definition is given in (39). The
upper straight line validates our choice of this constant as all the points
are under this line.

##  4 Modified Local Error Estimators and the Adaptive Algorithms

In this section, we propose several modified residual based local error
estimators based on the theoretical results derived in the last section and
design the algorithm for the a posteriori error control and adaptive
computations for the system we consider. We note that our adaptive procedure
is the classic one, namely

$\mbox{Solve}~{}\rightarrow~{}\mbox{Estimate}~{}\rightarrow~{}\mbox{Mark}~{}\rightarrow~{}\mbox{Refine},$

and a common practice is to include the adaptive algorithm and the numerical
results in one section. However, since the analysis in § 3.3 gives us plenty
of possibilities to design various modified error estimators and mesh
structures, we decide to separate out the elaboration of these estimators and
mesh structures in this independent section so that the ideas behind are well
understood.

###  4.1 Stress tensor correction

Before we go into the detail of the modified estimators and mesh structures,
we first briefly review the stress tensor correction which is a preliminary
step for specifying a correct a posteriori error estimator.

As we mentioned in § 2.2.1, the stress tensors $\sigma^{{\rm ac}}$ and
$\sigma^{\rm a}$, on which $\eta_{\rm tr}$, $\eta_{\rm mo}$, and $\eta_{\rm
cg}$ all depend, are unique only up to a divergence free tensor field whose
definition is given in (15). Therefore, we need to minimize
$\eta(y_{h})=\eta_{\rm tr}(y_{h})+\eta_{\rm mo}(y_{h})+\eta_{\rm cg}(y_{h})$
with respect to all the admissible stress tensors such that

$\langle\delta\mathscr{E}^{\rm a}(I_{\rm a}y_{h}),v\rangle\leq\min_{c_{\rm
a}\in N_{1}(\mathcal{T}_{\rm a})^{2},c_{h}\in
N_{1}(\mathcal{T}_{h})^{2}}\tilde{\eta}\big{(}\sigma^{\rm a}(I_{\rm
a}y_{h})+\nabla c_{\rm a}{\sf J},\sigma^{{\rm ac}}(y_{h})+\nabla c_{h}{\sf
J}\big{)}\|\nabla v\|_{L^{2}},$ (37)

where ${\sf J}$ is the rotation matrix defined in [33, § 2.4.3]. Instead of
minimizing the total error estimator $\eta$ with respect to $c_{\rm a}$ and
$c_{h}$ introduced in (37), [33] introduced an approximation of the stress
tensor correction which only requires minimizing the modeling error $\eta_{\rm
mo}$ with respect to the degrees of freedom of $\sigma^{{\rm ac}}$ adjacent to
the interface. This method dramatically reduces the computational cost of such
correction but greatly improves the accuracy of the estimator in the adaptive
computation.

For the sake of completeness, we give the stress tensor correction algorithm
which essentially follows from [33, Algorithm 1].

Algorithm 1 Stress tensor correction

  1. 1.

Take $\sigma^{\rm a}(I_{\rm a}y_{h})$ and $\sigma^{{\rm ac}}(y_{h})$ as the
canonical forms in (12) and (14) respectively.

  2. 2.

Denote $q_{f}$ as the midpoint of $f\in\mathcal{F}_{h}$. $c_{h}$ minimizes the
following sum

$\sum_{T\in\mathcal{T}^{\rm i}}|T|\left[\sigma^{\rm a}(I_{\rm
a}y_{h},T)-\big{(}\sigma^{{\rm ac}}(I_{\rm a}y_{h},T)+\nabla c_{h}{\sf
J}\big{)}\right]^{2}$ (38)

subject to the constrain such that $c_{h}(q_{f})=0$, for $f\bigcap\Lambda_{\rm
i}=\emptyset$.

  3. 3.

Let $\sigma^{{\rm ac}}(y_{h})=\sigma^{{\rm ac}}(y_{h})+\nabla c_{h}{\sf J}$,
compute $\eta_{\rm mo}$, $\eta_{\rm tr}$ and $\eta_{\rm cg}$ with $\sigma^{\rm
a}(I_{\rm a}y_{h})$ and $\sigma^{{\rm ac}}(y_{h})$.

###  4.2 Modified local error estimators and the adaptive algorithm

We are finally at the stage of proposing various modified local a posteriori
error estimators and design our main adaptive algorithm. Before that, we
should first recall the difficulty we face and the theoretical results in hand
that may be able to deal with the difficulty so that the ideas behind the
modified error estimators are easily understood.

Our main difficulty comes from the overwhelming cost of computing the
elementwise modeling residual defined in (26) which is due to the mismatch of
$\mathcal{T}_{{\rm a}}$ and $\mathcal{T}_{h}$ in $\Omega^{\rm c}$ as
illustrated in Figure 2. We then prove in Theorem 2 and Corollary 3 that,
under suitable assumptions, the elementwise modeling residual is proportional
to the elementwise coarsening residual by a factor of $O({h^{-1}_{T}})$ in the
continuum region, where $h_{T}$ is the diameter of the element. Moreover, we
numerically verify that this ratio relationship holds for the cases of defects
we consider even without the stringent assumptions in the theorem and
corollary. We may conclude that both the theoretical analysis and the
numerical evidence show that, in principle, we are able to approximate the
elementwise modeling residual by its coarsening counterpart. However, we
should be conservative to such approximation around the atomistic-to-continuum
interface, since a good number of researches in the analysis of the a/c
coupling method have shown that the interface error may be large and subtly
important for the construction of an efficient a/c mesh [3, 21, 33].

To justify our concern, we show in Figure 5 the decay of the modeling residual
and the coarsening residual given in (23) and (22) respectively and those in
different regions, with respect to the increasing number of degrees of
freedom. The computations use the original elementwise modeling residual and
its coarsening counterpart defined in (26) and (25) respectively. Figure 5
shows that, somehow as we predicted, the overall modeling residual
(represented by the red line with asterisk markers) is almost negligible
(about two orders smaller) compared with the overall coarsening residual
(represented by the dark blue line with round markers). However, the modeling
residual (represented by the purple line with square markers) and its
coarsening counterpart (represented by the light blue line with triangle
markers) are comparable around the atomistic-to-continuum interface which
affirms that the modeling residual may not be neglected or even be
approximated at the interface. In addition, we clearly observe that the
modeling residual concentrates at the interface which perfectly coincides with
various earlier analysis of the interface influence of different coupling
schemes [3, 9, 10]. The numerical observations listed above essentially
motivate us to propose the following constructions of the modified a
posteriori error estimators and mesh structures.

(a) Micro-crack

(b) Anti-plane screw dislocation

(c) Multiple vacancies

Figure 5: The modeling residual and the coarsening residual in different
regions.

To keep the accuracy of the interface modeling residual while reducing the
computational cost, we introduce an “interface buffer” region $\Omega^{\rm
buf}$ such that

$\displaystyle\Omega^{\rm buf}:=\bigcup_{T\in\mathcal{T}^{\rm
buf}_{h}}T\qquad\textrm{with}\quad\mathcal{T}^{\rm
buf}_{h}:=\\{T\in\mathcal{T}^{{\rm c}}_{h}~{}|~{}\exists
T^{\prime}\in\mathcal{T}_{{\rm a}}~{}~{}\textrm{s.t.}~{}~{}T\cap
T^{\prime}=T^{\prime}\\}.$ (39)

We denote $\Lambda^{\rm buf}:=\Lambda\cap\Omega^{\rm buf}$. See Figure 6 for
an illustration of the structure of the mesh $\mathcal{T}_{h}$ and the
“interface buffer” region $\Omega^{\rm buf}$ (the grey region) for a single
vacancy.

Figure 6: Illustration of “interface buffer” region $\Omega^{\rm buf}$ which
is the domain marked grey.

We note that the modeling residual can be computed exactly by (26) in
$\Omega^{\rm buf}$ without the need of finding the geometric relation of
$\mathcal{T}_{{\rm a}}$ and $\mathcal{T}_{h}$ as they coincide in this region
and thus the save of the computational cost. However, the length of
$\Omega^{\rm buf}$, $R_{\rm buf}$ should be kept small since it will introduce
unnecessary degrees of freedom otherwise.

Outside $\Omega^{\rm buf}$, namely, the region $\Omega\setminus(\Omega^{\rm
a}\cup\Omega^{\rm buf})$, we apply certain modification for the elementwise a
posteriori error estimator. A smooth way to introduce the modification is to
employ a blending function $\beta\in C^{0}(\mathbb{R};\mathbb{R})$ so that the
blended approximated modeling residual $\tilde{\eta}_{\rm mo}(u_{h};T)$ for
$T\in\mathcal{T}_{h}^{\rm c}\setminus\mathcal{T}_{h}^{\rm buf}$ is defined by

$\displaystyle\tilde{\eta}_{\rm mo}(u_{h};T):=\beta(r^{\rm
buf}_{T})\cdot\eta_{\rm mo}(u_{h};T)+\big{(}1-\beta(r^{\rm
buf}_{T})\big{)}\cdot\frac{C^{\rm a-post}_{\rm mo-cg}}{h_{T}}\eta_{\rm
cg}(u_{h};T),$ (40)

where $r^{\rm buf}_{T}:={\rm dist}(T,\Lambda^{\rm buf}):=\inf\\{|\ell-
x_{T}|,\forall\ell\in\Lambda^{\rm buf}\\}$ is the distance between $T$ and
$\Lambda^{\rm buf}$ and $x_{T}$ is the barycenter of $T$. $C^{\rm a-post}_{\rm
mo-cg}$ is the constant defined in (36) which is in practice determined
empirically (cf. Figure 4). Given the blending width $R_{\rm bld}>0$, a
typical linear blending function $\beta\in C^{0}(\mathbb{R};\mathbb{R})$ is
defined by

$\displaystyle\beta(r):=\left\\{\begin{array}[]{ll}1,&r\geq R_{\rm bld}\\\
\frac{r-1}{R_{\rm bld}-1},&1<r<R_{\rm bld}\\\ 0,&r\leq 1\end{array}\right..$
(44)

A more straightforward way for approximating the elementwise modeling error
residual for $T\in\mathcal{T}^{\rm c}_{h}\setminus\mathcal{T}_{h}^{\rm buf}$
is to directly define

$\displaystyle\tilde{\eta}_{\rm mo}(u_{h};T):=\frac{C^{\rm a-post}_{\rm mo-
cg}}{h_{T}}\eta_{\rm cg}(u_{h};T).$ (45)

The advantage of such approximation is apparent that no geometric
consideration is involved even for a small layer of elements near the
“interface buffer” region.

We denote $\tilde{\eta}^{2}_{\rm mo}(u_{h}):=(C^{\rm
tr})^{2}\sum_{T\in\mathcal{T}_{h}}\tilde{\eta}^{2}_{\rm mo}(u_{h};T)$. Recall
that the local coarsening residual is already defined by (25). We thus are
ready to define the modified elementwise local estimator $\rho_{T}(u_{h};T)$
for $T\in\mathcal{T}_{h}$ as

$\rho_{T}(u_{h};T):=(C^{\rm tr})^{2}\frac{\tilde{\eta}_{\rm
mo}^{2}(u_{h};T)}{\tilde{\eta}_{{\rm mo}}(u_{h})}+(\sqrt{3}C^{\rm
tr}C_{\mathcal{T}_{h}})^{2}\frac{\eta^{2}_{{\rm cg}}(u_{h};T)}{\eta_{\rm
cg}(u_{h})}.$ (46)

It is straightforward to see that
$\sum_{T\in\mathcal{T}_{h}}\rho_{T}(u_{h};T)=\tilde{\eta}_{\rm
mo}(u_{h})+\eta_{\rm cg}(u_{h})$. We note that the constants $C^{\rm tr}$ and
$C_{\mathcal{T}_{h}}$ are not known a priori and we use empirical estimates in
the computation instead.

As mentioned earlier, we do not assign the truncation residual to local
contributions but compute it by (21) directly and check its value during the
adaptive process since it is much smaller than the total local estimator
$\rho_{T}$ when $\Omega$ is sufficiently large (cf. § 3.3).

We give the main adaptive algorithm below. The algorithm essentially follows
from [33, Algorithm 3] but we adapt it to the current setting.

Algorithm 2 A posteriori mesh refinement with control of the computational
domain.

  1. Step 0

Prescribe: Set $\Omega$, $\mathcal{T}_{h}$, $N_{\max}$, $\rho_{\rm tol}$,
$\tau_{1}$, $\tau_{2}$ and $R_{\max}$.

  2. Step 1

Solve: Solve the GRAC solution $u_{h}$ of (8) on the current mesh
$\mathcal{T}_{h}$.

  3. Step 2

Estimate: Carry out the stress tensor correction step in Algorithm 1, compute
the error estimator $\rho_{T}$ by (46) for each $T\in\mathcal{T}_{h}$, and the
truncation error estimator $\eta_{\rm tr}$ by (21). Compute the degrees of
freedom $N$ and $\rho=\sum_{T}\rho_{T}$. If $\eta_{\rm tr}>\tau_{2}\rho$,
enlarge the computational domain $\Omega$. Stop if $N>N_{\max}$ or
$\rho<\rho_{\rm tol}$ or $R>R_{\max}$.

  4. Step 3

Mark:

    1. Step 3.1

: Choose a minimal subset $\mathcal{M}\subset\mathcal{T}_{h}$ such that

$\sum_{T\in\mathcal{M}}\rho_{T}\geq\textrm{mean}(\rho).$

    2. Step 3.2

: We can find the interface elements which are within $k$ layers of atomistic
distance, $\mathcal{M}^{k}_{\rm
i}:=\\{T\in\mathcal{M}\bigcap\mathcal{T}_{h}^{\rm
c}:{\rm{dist}}(T,\Lambda^{{\rm i}})\leq k\\}$. Choose $K\geq 1$, find the
first $k\leq K$ such that

$\sum_{T\in\mathcal{M}^{k}_{\rm
i}}\rho_{T}\geq\tau_{1}\sum_{T\in\mathcal{M}}\rho_{T},$ (47)

with tolerance $0<\tau_{1}<1$. If such a $k$ can be found, let
$\mathcal{M}=\mathcal{M}\setminus\mathcal{M}^{k}_{\rm i}$. Then in step 3,
expand the interface $\Lambda^{\rm i}$ and the “interface buffer”
$\Lambda^{\rm buf}$ outward by $k$ layers.

  5. Step 4

Refine: Bisect all elements $T\in\mathcal{M}$. Go to Step 1.

##  5 Numerical Experiments

In this section, we conduct the adaptive computations for the three types of
defects, namely micro crack, anti-plane screw dislocation and multiple
vacancies, based on our adaptive algorithm proposed in Algorithm 2 with
different a posteriori error estimators and mesh structures. The geometries of
these defects are visualized in Figure 1.

We consider a two dimensional triangular lattice $\Lambda^{\rm
hom}:=\mathsf{A}\mathbb{Z}^{2}$ with

$\displaystyle\mathsf{A}=\Bigg{(}\begin{array}[]{cc}1&1/2\\\
0&\sqrt{3}/2\end{array}\Bigg{)}.$ (50)

Note that upon rotating and possibly dilating, the projection of the BCC
crystal is in fact a triangular lattice [8] and hence it is of practical
interest to be considered.

We adopt the well-known EAM potential as the site energy potential $V_{\ell}$
through our numerical experiments, where

$\displaystyle V_{\ell}(y):=$
$\displaystyle\sum_{\ell^{\prime}\in\mathcal{N}_{\ell}}\phi(|y(\ell)-y(\ell^{\prime})|)+F\Big{(}{\textstyle\sum_{\ell^{\prime}\in\mathcal{N}_{\ell}}\psi(|y(\ell)-y(\ell^{\prime})|)}\Big{)}$
$\displaystyle=$
$\displaystyle\sum_{\rho\in\mathcal{R}_{\ell}}\phi\big{(}|D_{\rho}y(\ell)|\big{)}+F\Big{(}{\textstyle\sum_{\rho\in\mathcal{R}_{\ell}}}\psi\big{(}|D_{\rho}y(\ell)|\big{)}\Big{)},$
(51)

for a pair potential $\phi$, an electron density function $\psi$ and an
embedding function $F$. In particular, we choose

$\phi(r)=\exp(-2a(r-1))-2\exp(-a(r-1)),\quad\psi(r)=\exp(-br),$ $\text{ and
}\quad
F(\tilde{\rho})=C\left[(\tilde{\rho}-\tilde{\rho_{0}})^{2}+(\tilde{\rho}-\tilde{\rho_{0}})^{4}\right],$

with parameters $a=4,b=3,c=10$ and $\tilde{\rho_{0}}=6\exp(0.9b)$, which is
exactly the same as the numerical experiments in [35] and [33].

We test the performance of the modified a posteriori error estimators and the
mesh structures with different parameters, i.e. the width of the buffer region
and the width of the blending region, which we believe may influence the
adaptive process. In particular, we carry out the adaptive computations for
the modified elementwise local estimator $\rho_{T}(u_{h};T)$ defined in (46)
first by using the directly approximated modeling residual $\tilde{\eta}_{\rm
mo}(u_{h};T)$ defined in (45) with the widths of the buffer region being
$R_{\rm buf}=1,3$ and $10$. We then fix $R_{\rm buf}=3$ and apply the blended
approximated modeling residual defined in (40) with the widths of the blending
region $R_{\rm bld}=2$ and $6$. To better reveal the effectiveness and the
efficiency of our error estimators and mesh structures, we compare the results
of our adaptive computations with those using an a priori graded mesh [8, 21],
the original residual based a posteriori error estimator given in (26) and
(25) and also the purely coarsening error estimator by totally neglecting the
modeling residual.

In all numerical experiments, the radius of the computational domain $\Omega$
is set to be $300$ initially. The adaptive processes start with the initial
configuration with $R_{\rm a}=6$, where for the case of multiple vacancies, it
means that the initial radius of each disjoint atomistic region is 6. The
adaptive parameters in Algorithm 2 are fixed to be $\tau_{1}=0.7$ and
$\tau_{2}=1.0$.

###  5.1 Micro-crack

The first defect we consider is the micro-crack, which is a prototypical
example of point defects. We note that, technically speaking, this is not a
“crack”. It only serves as an example of a localized defect with an
anisotropic shape. To generate this defect, we remove $k$ atoms from
$\Lambda^{\rm hom}$,

$\displaystyle\Lambda_{k}^{\rm def}:=\\{-(k/2)e_{1},\ldots,(k/2-1)e_{1})\\},$
$\displaystyle\qquad{\rm if}\quad k\quad\textrm{is even},$
$\displaystyle\Lambda_{k}^{\rm
def}:=\\{-(k-1)/2e_{1},\ldots,(k-1)/2e_{1})\\},$ $\displaystyle\qquad{\rm
if}\quad k\quad\textrm{is odd},$

and $\Lambda=\Lambda^{\rm hom}\setminus\Lambda_{k}^{\rm def}$. In our example,
we set $k=11$. We apply an isotropic stretch $\mathrm{S}$ and a shear
$\gamma_{II}$ by setting

${\sf B}=\left(\begin{array}[]{cc}1&\gamma_{II}\\\
0&1+\mathrm{S}\end{array}\right)\cdot{\sf{F_{0}}},$

where ${\sf F_{0}}\propto\mathrm{I}$ is a macroscopic stretch or compression
and $\mathrm{S}=\gamma_{II}=0.03$.

Figure 7a shows the convergence of the true error $\|\nabla u_{h}-\nabla
u\|_{L^{2}(\Omega)}$ with respect to the number of degrees of freedom $N$ and
Figure 7b shows the efficiency factors (which is defined by the ratios of the
error estimators and the true error) for various a posteriori error estimators
and mesh structures.

(a) Convergence

(b) Efficiency factors

Figure 7: The convergences of the error and the efficiency factors for
different error estimators with respect to the number of degrees of freedom
for the micro-crack.

First of all, Figure 7a clearly verifies the necessity of the prudence of
treating the modeling residual estimator around the interface. In terms of the
error versus the number of degrees of freedom, all adaptive computations
achieve the same optimal convergence rate compared with the a prior graded
mesh and the original residual based a posteriori error estimator and are
barely distinguishable if we fix the width of the buffer region to be 3, i.e.,
$R_{\rm buf}=3$. Totally removing the modeling residual from the a posteriori
error estimator (represented by the green line with round markers) or crudely
approximating it near the interface (represented by the light blue line with
six-pointed star markers for $R_{\rm buf}=1$) apparently result in a lower
convergence rate. However, we may not need to be too cautious either, since
the influence of the modeling error decays rapidly [46, 47, 48]. This is
demonstrated by the dark blue line with star markers for $R_{\rm buf}=10$,
where apparently unnecessary degrees of freedom are placed around the
interface that causes a even worse rate of convergence even compared with the
crude approximation.

Figure 7b plots the efficiency factors for various error estimators and mesh
structures. The factors are uniformly moderate for all the estimators with
some oscillation which we believe has little effect for the adaptivity. We
note that rigorous efficiency analysis, which proves that the residual based
error estimator provides both upper and lower bounds for the true error up to
some constants, requires a substantial amount of additional technicality and
should be included in a separate work. We refer to [49] for a theoretical
consideration of the efficiency for the GRAC method in one dimension that may
be the only existing literature in this field to the best knowledge of the
authors.

(a) CPU times

(b) Time constitution

Figure 8: The CPU times for each step in the adaptive computations for the
micro-crack.

The most interesting observation that illustrates the significance of the
current work can be concluded in Figure 8. Figure 8a plots the total CPU times
versus the numbers of degrees of freedom $N$ for various error estimators and
mesh structures while Figure 8b plots those for computing the error estimators
and solving the a/c coupling problems separately. We see clearly in Figure 8a
that the CPU time of adaptively solving the a/c problem using the original a
posteriori error estimator scales as $O(N^{2})$ whereas those using all other
estimators scales as $O(N)$, except for the case of $R_{{\rm bld}}=6$ (a
slightly large blending region). This coincides with our estimate of the
computational cost in § 3.3 and evidently shows that the original a posteriori
error estimator (represented by the blue line with asterisk markers) may not
be employed for the purpose of practical adaptive computations. It is also
interesting to see that although we observe in Figure 7a an optimal
convergence rate of the error using the blended approximated estimator, Figure
8a shows that such construction may gild the lily, even when the blending
region is just slightly large (represented by the purple line with triangle
markers for $R_{\rm bld}=6$). It may be explained by the fact that all the
computations for the modeling residual, including computing the atomistic
stress tensor and finding the geometric relation between $\mathcal{T}_{\rm a}$
and $\mathcal{T}_{h}$, are still carried out in the blending region which
increase the computational cost rapidly with the growing width of the blending
region. The dark blue line with star markers, which represents the mesh
structure having a large buffer region ($R_{\rm buf}=10$), is noticed to be
slightly higher than other “$O(N)$” curves. This shows that computing the
modeling residual is still more expensive than the coarsening residual for the
reason that the modeling residual is ultimately nonlocal. This indicates that
we may avoid computing the atomistic stress tensor as much as possible if the
interface effect is properly captured for the purpose of efficiency.

A more straightforward comparison of the original error estimator and the
modified one is given in Figure 8b. It shows that the additional work purely
comes from the computing the original estimator but the cost of computing the
modified estimator is somehow marginal compared with that of solving the a/c
problems themselves.

###  5.2 Anti-plane Screw Dislocation

The second defect for which we conduct the adaptive computation using
Algorithm 2 is the anti-plane screw dislocation. Following [8], we restrict
ourselves to an anti-plane shear motion which is illustrated in Figure 1b. In
this case, the Burgers vector is set to be $b=(0,0,1)^{T}$ and the center of
the dislocation core is chosen to be $\hat{x}:=\frac{1}{2}(1,1,\sqrt{3})^{T}$,
and we assume that no additional shear deformation is applied. Thus the
unknown for this dislocation model is the displacement in the
$e_{3}$-direction. As we mentioned at the beginning of § 2.1, we choose
$u_{0}\in\mathscr{U}$ given by the solution of a linear elasticity model as a
far-field crystalline environment or a predictor in the language of numerical
computation and we need to compute $u$ which can be considered as a corrector.
For the sake of completeness, we review the detailed derivation of $u_{0}$ in
Appendix A.

Similar results are observed for the anti-plane screw dislocation as those in
the case of micro-crack with half of the convergence rate as predicted in [8].
It is interesting to see that the adaptive computation using the pure
coarsening residual (represented by the green line with round markers) again
results in half of the optimal convergence rate (this time being DOF to the
power of $0.25$ rather than to the power of $0.5$). Whether it is a
coincidence needs a further investigation.

(a) Convergence

(b) Efficiency factors

Figure 9: The convergences of the error and the efficiency factors for
different error estimators with respect to the number of degrees of freedom
for the anti-plane screw dislocation.

(a) CPU times

(b) Time constitution

Figure 10: The CPU times for each step in the adaptive computations for the
anti-plane screw dislocation.

###  5.3 Multiple Vacancies

The last type of defect we consider is the case of multiple vacancies. We
select three vacancy sites which are separated away from each other by a
distance of $40$ (see Figure 1c for an illustration). We then apply the same
isotropic stretch $\mathrm{S}$ and shear $\gamma_{II}$ to this defected system
as that in the micro-crack. This is a very simple model for the interactions
of defects and we choose this case to see if the adaptive a/c coupling method
can be potentially applied to study such problems.

We show in Figure 11 the evolution of the atomistic and continuum partitioning
in the adaptive process. Unlike the micro-crack and the anti-plane screw
dislocation where only one atomistic region exists throughout, we start with
three disjoint atomistic regions initially which contain the three separated
vacancies. The atomistic regions then grow and get very close as the
adaptivity goes and finally merge into one “big” atomistic region. We see in
the fourth subplot Figure 11 that the sizes of the mesh between the atomistic
regions also shrink so small that we should adaptively shift the continuum
model to the atomistic model.

Figure 11: Illustration of the evolution of atomistic and continuum partition
in the adaptive process for the multiple vacancies

We observe in Figure 12a a sudden drop of the error in each line which
corresponds to the merge of the disjoint atomistic regions though the lines
all possess their own “correct” convergence rate. We put the term “correct” in
a quotation mark because there is actually no predicted convergence rate for
this case before the merge of the atomistic regions. This is also the reason
why there is no line representing the graded mesh, whose construction is often
based on such prediction, is plotted in Figure 12a as opposed in Figure 7a and
Figure 9a. However, this may be a more common situation we need to face when
dealing with realistic material systems. It is interesting to discuss the
issue of how we might adjust the adaptive strategy so that the time spot of
the merge of the atomistic is (quasi-)optimal though it is beyond the scope of
the current work and may be included in a future numerical study of the
adaptive a/c method.

(a) Convergence

(b) Efficiency factors

Figure 12: The convergences of the error and the efficiency factors for
different error estimators with respect to the number of degrees of freedom
for the multiple vacancies.

(a) CPU times

(b) Time constitution

Figure 13: The CPU times for each step in the adaptive computations for the
multiple vacancies.

The final figure in this paper once again shows the advantage of using the
properly modified a posteriori error estimator. We observe a slightly worse
rate for the blended estimator (the purple line with triangle markers with the
rate of $1.6$). This might be caused by the complex geometry between the
atomistic regions, which confirms our belief that any unnecessary search for
the geometric relationship between $\mathcal{T}_{\rm a}$ and $\mathcal{T}_{h}$
should be avoided.

##  6 Conclusion and Outlook

We develop the efficient adaptive algorithm for the geometric reconstruction
based consistent atomistic/continuum (GRAC) coupling scheme based on a
modified residual based a posteriori error estimator for various types of
crystalline defects. The key for achieving such efficiency but maintaining
comparable accuracy is the strategy of constructing an a/c coupling mesh with
a small “interface buffer” region around the coupling interface in which the
elementwise modeling residual is computed exactly but outside which the
modeling residual is replaced by a theoretically and numerically justified
approximation using its coarsening counterpart. Although we believe such
strategy is generally applicable for a/c coupling schemes and more complex
crystalline defects, the current research still raises a few open problems
which deserve further mathematical analysis and algorithmic developments.

The first one is the extension to general short range interactions [34]. To
generalize our strategy to that setting, a detailed study of the width of
“interface buffer” region with respect to the interaction range is needed.

The second one is the problem of adaptive error control for energy blended a/c
coupling schemes, such as energy-based blended quasi-continuum (BQCE) method
and blended ghost force correction (BGFC) method. The a priori analysis in [9,
10, 20] provide an analytical framework and the stress tensor based
formulation plays a key role in the analysis. By inheriting such framework and
using the modified error estimator proposed in this paper, the adaptivity for
those coupling schemes should be somehow straightforward for the case of fixed
width of the blending region. However, a considerable amount of effort may be
required to problem of the adaptive choice of the width of the blending
region.

The third one is the extension to three dimensional problems. We note that
constructing a consistent a/c coupling scheme with a sharp interface in three
dimensions still remains an open problem. Thus the development of adaptive
energy blended methods mentioned in the previous paragraph may be the
appropriate direction.

The last but may be the most important problem to consider is the adaptive
computation for complex crystalline defects. The simulation of the
interactions of defects, which is briefly touched by the third numerical
experiment, and the dislocation nucleation have already attracted considerable
attention. It is very difficult (if not impossible) to give rigorous a priori
analysis for such problems, including the proper boundary conditions and
complicated geometry at the interface. However, it should be the place where
the adaptive a/c coupling methods reveal its most advantage and power and we
believe that the idea of the current research makes it possible to develop
efficient and robust adaptive algorithms for those problems with practical
interests.

##  Appendix Appendix A Far-field boundary predictor for anti-plane screw
dislocations

We model dislocations by following the setting used in [8]. We consider a
model for straight dislocations obtained by projecting a 3D crystal into 2D.
Let $B\in\mathbb{R}^{3\times 3}$ be a nonsingular matrix. Given a Bravais
lattice $B\mathbb{Z}^{3}$ with dislocation direction parallel to $e_{3}$ and
Burgers vector ${\sf b}=({\sf b}_{1},0,{\sf b}_{3})$, we consider
displacements $W:B\mathbb{Z}^{3}\rightarrow\mathbb{R}^{3}$ that are periodic
in the direction of the dislocation direction of $e_{3}$. Thus, we choose a
projected reference lattice
$\Lambda:=A\mathbb{Z}^{2}:=\\{(\ell_{1},\ell_{2})~{}|~{}\ell=(\ell_{1},\ell_{2},\ell_{3})\in
B\mathbb{Z}^{3}\\}$. We also introduce the projection operator

$P(\ell_{1},\ell_{2})=(\ell_{1},\ell_{2},\ell_{3})\quad\text{for}~{}\ell\in
B\mathbb{Z}^{3}.$ (52)

It can be readily checked that this projection is again a Bravais lattice. For
anti-plane screw dislocation, $\Lambda$ is obtained as projection of a 3D
Bravais lattice along the screw dislocation direction (and the direction of
slip) and we restrict the displacements of the form $u=(0,0,u_{3})$.

We follow the constructions in [8, 50] for modeling dislocations and prescribe
$u_{0}$ as follows. Let $\Lambda\subset\mathbb{R}^{2}$,
$\hat{x}\in\mathbb{R}^{2}$ be the position of the dislocation core and
$\Upsilon:=\\{x\in\mathbb{R}^{2}~{}|~{}x_{2}=\hat{x}_{2},~{}x_{1}\geq\hat{x}_{1}\\}$
be the “branch cut”, with $\hat{x}$ chosen such that
$\Upsilon\cap\Lambda=\emptyset$.

We define the far-field predictor $u_{0}$ by solving the continuum linear
elasticity (CLE)

$\displaystyle\mathbb{C}^{j\beta}_{i\alpha}\frac{\partial^{2}u^{\rm
lin}_{i}}{\partial x_{\alpha}\partial x_{\beta}}$ $\displaystyle=$
$\displaystyle 0\qquad\text{in}~{}~{}\mathbb{R}^{2}\setminus\Upsilon,$
$\displaystyle u^{\rm lin}(x+)-u^{\rm lin}(x-)$ $\displaystyle=$
$\displaystyle-{\sf
b}\qquad\text{for}~{}~{}x\in\Upsilon\setminus\\{\hat{x}\\},$ (53)
$\displaystyle\nabla_{e_{2}}u^{\rm lin}(x+)-\nabla_{e_{2}}u^{\rm lin}(x-)$
$\displaystyle=$ $\displaystyle
0\qquad\text{for}~{}~{}x\in\Upsilon\setminus\\{\hat{x}\\},$

where the forth-order tensor $\mathbb{C}$ is the linearised Cauchy-Born tensor
(derived from the potential $V$, see [8, § 7] for more detail).

We mention that for the anti-plane screw dislocation, under the proper
assumptions on the interaction range $\mathcal{R}$ and the potential $V$, the
first equation in (Appendix A) simply becomes to $\Delta u^{\rm lin}=0$ [50].
The system (Appendix A) then has the well-known solution

$\displaystyle u_{0}(x):=u^{\rm lin}(x)=\frac{{\sf b}}{2\pi}\arg(x-\hat{x}),$
(54)

where we identify $\mathbb{R}^{2}\cong\mathbb{C}$ and use $\Upsilon-\hat{x}$
as the branch cut for arg.

Note that for the purpose of analysis, we have $\nabla u_{0}\in
C^{\infty}(\mathbb{R}^{2}\setminus\\{0\\})$ and $|\nabla^{j}u_{0}|\leq
C|x|^{-j}$ for all $j\geq 0$ and $x\neq 0$.

## References

  * [1] E. V. D. Giessen, P. A. Schultz, N. Bertin, V. V. Bulatov, W. Cai, G. Csányi, S. M. Foiles, M. G. Geers, C. González, M. Hütter, et al., Roadmap on multiscale materials modeling, Model. Simul. Mat. Sci. Eng. 28 (4) (2020) 043001.
  * [2] E. B. Tadmor, R. E. Miller, Modeling Materials: Continuum, Atomistic and Multiscale Techniques, Cambridge University Press, 2012.
  * [3] M. Luskin, C. Ortner, Atomistic-to-continuum coupling, Acta Numer. 22 (2013) 397\.
  * [4] M. Ortiz, R. Phillips, E. B. Tadmor, Quasicontinuum analysis of defects in solids, Philos. Mag. A 73 (6) (1996) 1529–1563.
  * [5] R. Miller, E. Tadmor, A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods, Modelling Simul. Mater. Sci. Eng. 17.
  * [6] D. Olson, P. Bochev, M. Luskin, A. Shapeev, An optimization-based atomistic-to-continuum coupling method, SIAM J. Numer. Anal. 52 (4) (2014) 2183–2204.
  * [7] D. Olson, A. Shapeev, P. Bochev, M. Luskin, Analysis of an optimization-based atomistic-to-continuum coupling method for point defects, ESAIM: Math. Model. Numer. Anal. 50 (1) (2016) 1–41.
  * [8] V. Ehrlacher, C. Ortner, A. V. Shapeev, Analysis of boundary conditions for crystal defect atomistic simulations, Arch. Rat. Mech. Anal. 222 (3) (2016) 1217–1268.
  * [9] X. H. Li, C. Ortner, A. Shapeev, B. V. Koten, Analysis of blended atomistic/continuum hybrid methods, Numer. Math. 134 (2016) 275–326.
  * [10] C. Ortner, L. Zhang, Atomistic/continuum blending with ghost force correction, SIAM J. Sci. Comput. 38 (1) (2016) A346–A375.
  * [11] H. B. Dhia, L. Chamoin, J. T. Oden, S. Prudhomme, A new adaptive modeling strategy based on optimal control for atomic-to-continuum coupling simulations, Comput. Methods Appl. Mech. Engrg. 200 (37-40) (2011) 2675–2696.
  * [12] A. S. Davis, J. T. Lloyd, V. Agrawal, Moving window techniques to model shock wave propagation using the concurrent atomistic–continuum method, Comput. Methods Appl. Mech. Engrg. 389 (2022) 114360.
  * [13] S. Chakraborty, S. Ghosh, Crystal plasticity phase-field model with crack tip enhancement through a concurrent atomistic-continuum model, J. Mech. Phys. Solids 155 (2021) 104563.
  * [14] J. Tinsley Oden, S. Prudhomme, A. Romkes, P. T. Bauman, Multiscale modeling of physical phenomena: Adaptive control of models, SIAM J. Sci. Comput. 28 (6) (2006) 2359–2389.
  * [15] J. T. Oden, Adaptive multiscale predictive modelling, Acta Numerica 27 (2018) 353–450.
  * [16] S. Chakraborty, S. Ghosh, A concurrent atomistic-crystal plasticity multiscale model for crack propagation in crystalline metallic materials, Comput. Methods Appl. Mech. Engrg. 379 (2021) 113748.
  * [17] P. Gupta, M. Ortiz, D. M. Kochmann, Nonequilibrium thermomechanics of gaussian phase packet crystals: Application to the quasistatic quasicontinuum method, J. Mech. Phys. Solids 153 (2021) 104495.
  * [18] P. Ming, J. Z. Yang, Analysis of a one-dimensional nonlocal quasi-continuum method, Multiscale Model. Simul. 7 (4) (2009) 1838–1875.
  * [19] C. Ortner, H. Wang, A priori error estimates for energy-based quasicontinuum approximations of a periodic chain, Math. Models Methods Appl. Sc. 21 (2011) 2491–2521.
  * [20] J. Lu, P. Ming, Convergence of a force-based hybrid method for atomistic and continuum models in three dimension, Comm. Pure Appl. Math. 66 (2013) 83–108.
  * [21] C. Ortner, L. Zhang, Construction and sharp consistency estimates for atomistic/continuum coupling methods with general interfaces: a 2D model problem, SIAM J. Numer. Anal. 50 (2012) 2940–2965.
  * [22] P. Lin, Convergence analysis of a quasi-continuum approximation for a two-dimensional material without defects, SIAM J. Numer. Anal. 45 (1) (2007) 313–332.
  * [23] J. T. Oden, S. Prudhomme, Estimation of modeling error in computational mechanics, J. Comput. Phys. 182 (2) (2002) 496–515.
  * [24] S. Prudhomme, P. T. Bauman, J. T. Oden, Error control for molecular statics problems, Int. J. Multiscale Comput. Eng. 4 (5-6) (2006) 647–662.
  * [25] M. Arndt, M. Luskin, Goal-oriented atomistic-continuum adaptivity for the quasicontinuum approximation, Int. J. Multiscale Comput. Eng. 5 (5) (2007) 407–415.
  * [26] M. Arndt, M. Luskin, Error estimation and atomistic-continuum adaptivity for the quasicontinuum approximation of a Frenkel-Kontorova model, Multiscale Model. Simul. 7 (2008) 147–170.
  * [27] M. Arndt, M. Luskin, Goal-oriented adaptive mesh refinement for the quasicontinuum approximation of a Frenkel-Kontorova model, Comput. Methods Appl. Mech. Engrg. 197 (2008) 4298–4306.
  * [28] C. Ortner, A. V. Shapeev, Analysis of an energy-based atomistic/continuum coupling approximation of a vacancy in the 2d triangular lattice, Math. Comp. 82 (2013) 2191–2236.
  * [29] S. Prudhomme, L. Chamoin, H. B. Dhia, P. T. Bauman, An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations, Comput. Methods Appl. Mech. Engrg. 198 (21-26) (2009) 1887–1901.
  * [30] H. B. Dhia, G. Rateau, The arlequin method as a flexible engineering design tool, Int. J. Numer. Methods Eng. 62 (11) (2005) 1442–1462.
  * [31] C. Ortner, A priori and a posteriori analysis of the quasi-nonlocal quasicontinuum method in 1D, Math. Comp. 80 (275) (2011) 1265–1285.
  * [32] C. Ortner, H. Wang, A posteriori error control for a quasi-continuum approximation of a periodic chain, IMA J. Numer. Anal. 34 (3) (2013) 977–1001.
  * [33] H. Wang, M. Liao, P. Lin, L. Zhang, A posteriori error estimation and adaptive algorithm for atomistic/continuum coupling in two dimensions, SIAM J. Sci. Comput. 40 (4) (2018) A2087–A2119.
  * [34] M. Liao, P. Lin, L. Zhang, A posteriori error estimate and adaptive mesh refinement algorithm for atomistic/continuum coupling with finite range interactions in two dimensions, Commun. Comput. Phys. 27 (1) (2020) 198–226.
  * [35] C. Ortner, L. Zhang, Energy-based atomisitic-to-continuum coupling without ghost forces, Comput. Methods Appl. Mech. Engrg. 279 (1) (2014) 29–45.
  * [36] J. Braun, T. Hudson, C. Ortner, Asymptotic expansion of the elastic far-field of a crystalline defect, ArXiv e-prints 2108.04765.
  * [37] H. Chen, C. Ortner, Y. Wang, QM/MM methods for crystalline defects. part 3: Machine-learned interatomic potentials, ArXiv e-prints 2106.14559.
  * [38] M. S. Daw, M. I. Baskes, Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals, Phys.Rev.B 29 (12) (1984) 6443–6453.
  * [39] M. W. Finnis, J. E. Sinclair, A simple empirical n-body potential for transition-metals, Philos. Mag. A 50 (1) (1984) 45–55.
  * [40] W. E, P. Ming, Cauchy-Born rule and the stability of crystalline solids: static problems, Arch. Ration. Mech. Anal. 183 (2) (2007) 241–297.
  * [41] C. Ortner, F. Theil, Justification of the Cauchy-Born approximation of elastodynamics, Arch. Ration. Mech. Anal. 207 (3) (2013) 1025–1073.
  * [42] C. Ortner, The role of the patch test in 2D atomistic-to-continuum coupling methods, ESAIM Math. Model. Numer. Anal. 46 (2012) 1275–1319.
  * [43] P. Clément, Approximation by finite element functions using local regularization, ESAIM: Math. Model. Numer. Anal. 9 (R2) (1975) 77–84.
  * [44] R. Verfürth, Error estimates for some quasi-interpolation operators, ESAIM: Math. Model. Numer. Anal. 33 (4) (1999) 695–713.
  * [45] Y. Wang, H. Chen, M. Liao, C. Ortner, H. Wang, L. Zhang, A posteriori error estimates for adaptive qm/mm coupling methods, SIAM J. Sci. Comput. 43 (4) (2021) A2785–A2808.
  * [46] M. Dobson, M. Luskin, An analysis of the effect of ghost force oscillation on quasicontinuum error, ESAIM: Math. Model. Numer. Anal. 43 (3) (2009) 591–604.
  * [47] J. Chen, P. Ming, Ghost force influence of a quasicontinuum method in two dimension, J. Comput. Math. (2012) 657–683.
  * [48] L. Cui, P. Ming, The effect of ghost forces for a quasicontinuum method in three dimension, Sci. China Math. 56 (12) (2013) 2571–2589.
  * [49] H. Wang, S. Yang, Analysis of the residual-type and the gradient recovery-type a posteriori error estimators for a consistent atomistic-to-continuum coupling method in one-dimension, Multiscale Model. Simul. 16 (2) (2018) 679–709.
  * [50] J. Braun, M. Buze, C. Ortner, The effect of crystal symmetries on the locality of screw dislocation cores, SIAM J. Math. Anal. 51 (2019) 1108–1136.