# Adaptive Multiscale Coupling Methods of Molecular Mechanics based on a
Unified Framework of a Posteriori Error Estimates

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

###### Abstract

Multiscale coupling methods are significant methodologies for the modeling and
simulation of materials with defects, intending to achieve the (quasi-)optimal
balance of accuracy and efficiency. The a posteriori analysis and
corresponding adaptive algorithms play a crucial role in the efficient
implementation of multiscale coupling methods. This paper proposes a unified
framework for residual-based a posteriori error estimates that can be applied
to general consistent multiscale coupling methods. In particular, we prove
that the error estimator based on the residual force can provide the upper
bound of the true approximation error. As prototypical examples, we present a
variety of adaptive computations based on this reliable error estimator for
the blended atomistic-to-continuum (a/c) coupling methods, including the
energy-based blended quasi-continuum (BQCE), the force-based blended quasi-
continuum (BQCF) and the recently developed blended ghost force correction
(BGFC) methods. We develop a coarse-grained technique for the efficient
evaluation of the error estimator. A robust adaptive algorithm is therefore
proposed and validated with different types of crystalline defects, some of
which are not considered in previous related literature on the adaptive a/c
coupling methods. The results demonstrate that the adaptive algorithm leads to
the same optimal convergence rate of the error as the a priori error estimate,
but with considerable computational efficiency. This study provides valuable
insights into the design and implementation of adaptive multiscale methods,
and represents a significant contribution to the literature on a/c coupling
methods.

###### keywords:

atomistic-to-continuum coupling , a posteriori error estimate , adaptive
algorithm , crystalline defects

##  1 Introduction

In the past two decades, multiscale coupling methods have attracted great
attention from various academic communities, including those focused on
biochemistry, engineering, and mathematics [1, 2, 3, 4, 5, 6]. Atomistic-to-
continuum (a/c) coupling methods are a typical class of concurrent multiscale
schemes aiming to achieve the (quasi-)optimal balance of accuracy and
efficiency for modeling crystalline defects [7, 8, 9, 10, 11, 12]. The
fundamental concept underlying a/c coupling methods is to employ the more
precise atomistic model in the immediate vicinity of localized defects, while
the continuum model (such as the Cauchy-Born rule) is utilized in regions away
from the defect cores.

The a/c coupling schemes can be broadly categorised into sharp-interface
coupling and blended coupling methods. Each of these categories can further be
divided into energy-based (conservative) and force-based (non-conservative)
a/c couplings. The modeling and a priori analysis of a/c coupling methods have
received extensive and comprehensive investigations, for example, the QNL
(quasi-nonlocal quasicontinuum) method [13, 14], the BQCE (blended energy-
based quasi-continuum) method [12, 15], the BQCF (blended force-based quasi-
continuum) method [16, 17], the GRAC (geometric reconstruction based
atomistic/continuum coupling) method [18, 19] and the recently developed BGFC
(atomistic/continuum blending with ghost force correction) method [10, 20]. We
refer to [1, 21] for the extensive overview and benchmark of a/c coupling
methods and [7] for their rigorous a priori error analysis.

One of the primary obstacles facing a/c coupling methods is the optimal
allocation of atomistic and continuum regions, as well as the determination of
appropriate mesh structures, to achieve the (quasi-)optimal balance between
accuracy and efficiency. While feasible, the a priori choices generally result
in sub-optimal distribution of computational resources. Therefore, the a
posteriori analysis and corresponding adaptive algorithms are essential for
the efficient implementation and simulation of a/c coupling methods in real-
world material systems.

Several approaches have been proposed for adaptive a/c coupling methods, and
these methods vary in terms of how they construct the error estimators. The
first approach, known as the goal-oriented a posteriori error estimate, was
developed from a heuristic and engineered perspective [22, 23, 24, 25, 11].
The primary disadvantage of this approach is that the reliability of the error
estimator can not be guaranteed and more importantly, all these works utilize
the original energy-based atomistic-to-continuum (a/c) method as the
underlying model, which is known to be inconsistent due to the presence of the
so-called ghost force [7, 14, 26]. The second approach, referred to as the
residual-based a posteriori error estimate, adopts the idea of classic
adaptive finite element methods by providing a quantitative estimate of the
error in a particular energy-norm. This approach was first analyzed in one
dimension for consistent atomistic-to-continuum (a/c) coupling methods in [27,
28]. The extension to two dimensions for the GRAC method with nearest neighbor
interaction was conducted in [6] and further generalized to the case of finite
range interactions in [29]. However, these works suffered from the
inefficiency of the resulting adaptive algorithms due to the high
computational cost associated with evaluating the modeling residual. To
address this issue, our recent work [30] developed a theory-based
approximation for the residual-based a posteriori error estimator, resulting
in significant improvements in efficiency for the adaptivity computations.

Nonetheless, the limitations of the existing residual based a posteriori error
estimators for a/c coupling methods are significant. The complicated
construction of the error estimator based on the atomistic stress and the need
for additional steps such as the stress tensor correction [6] make the
approach less efficient and less tractable. Moreover, these estimators are
limited to the GRAC method, and there is a lack of transferability to other
methods. Finally, the absence of a rigorous a posteriori error estimate for
blended a/c coupling methods is a significant gap in the current research, as
such methods are commonly used in practice. These limitations highlight the
need for a new approach to a posteriori error estimation in a/c coupling
methods that can address these issues and provide more efficient and reliable
adaptive algorithms.

The purpose of the current work is to propose a unified framework of the a
posteriori error estimates for consistent atomistic-to-continuum coupling
methods which are used in designing the corresponding adaptive algorithms, and
implement the adaptive a/c methods for the simulations of several two
dimensional crystalline defects with practical importance. To be precise, our
contribution lies in the following three aspects.

First, we provide a framework to establish that the residual in the dual norm
of a given energy-norm is equivalent to the true error and to derive an error
estimator based on the residual force that provides an upper bound of the
residual. We show that such framework is unified such that it only depends on
the consistency condition that the solution of the coupling model converges to
the atomistic model but does not depend on the specific coupling scheme we
employ. The framework is pioneered by our previous research in the context of
QM/MM coupling methods [2] but only for point defects. In particular, the
proof of the equivalence of the approximation error and the residual in the
current work is adapted to a/c methods and is more quantitative and clear.
Moreover, we provide detailed estimates for the truncation error by analyzing
the decay of the residual that are crucial for determining the size of the
computational domain in the adaptive simulation of screw dislocation and crack
whose behaviors are substantially more complicated than point defects. The
efficient evaluation of local error contributions based on certain sampling
scheme is also presented.

Second, we design novel adaptive algorithms for several prototypical
atomistic-to-continuum coupling methods ranging from energy-based to force-
based, blended type to sharp interface, which could automatically adjust the
atomistic region and the continuous mesh structures as well as the blending
region, if applicable, on the fly. As opposed to the residual-stress based a
posteriori error estimates where the derivation for individual residual-stress
based estimator for each coupling scheme is needed [6, 29, 28, 31, 32], all
the adaptive algorithms in the current work are essentially based on the same
residual-force based estimator just given and only algorithmic adjustment is
needed. In addition, the adaptive algorithms for blended type of a/c methods
which are easily implemented and are thus widely adopted in the community of
mechanics and engineering [15, 12, 33], are in fact developed for the first
time in more than one dimension. We note that the residual-stress based a
posteriori error estimates for blended type of a/c methods are nontrivial even
in one dimension [32]. This clearly reveals the advantage of our unified
framework which provides considerable flexibility for the development of
adaptive algorithms of different coupling schemes.

Third, we provide numerical validation of our proposed adaptive algorithms for
a set of a/c coupling methods in the presence of various types of crystalline
defects, including crack which, to the best knowledge of the authors, is
implemented for the first time in the literature of a/c. The adaptive
simulations demonstrate that the proposed algorithms produce optimal
convergence rates of the error and (quasi-)optimal decomposition of the
domain. In addition, a thorough comparison between the adaptive algorithm
designed in the current work and a recently developed one based on a modified
residual-based error estimator, which has significant improvement in
computational cost but only limited to a/c methods with sharp interface [30],
is given. The result suggests that the approach proposed in the current work
achieves the same accuracy and efficiency but with much better effective
implementation again thanks to the identical formulation of the residual-force
based a posteriori error estimator. The systematic study of adaptive a/c
methods for cracks represents a pioneering effort, leveraging the inherent
attributes of the proposed residual-force based a posteriori error estimator.
Conversely, the a priori analysis poses challenges that require advanced
techniques to be addressed effectively. The inherent effectiveness and
remarkable flexibility of this approach not only broaden the horizons but also
provide a robust framework for investigating practical crystalline defects
such as grain boundaries through real-world simulations.

To ensure clarity of presentation, we focus on the atomistic system with
finite range interaction in two dimensions. However, the proposed unified
framework has the potential to be extended to other consistent multiscale
coupling methods and three-dimensional problems. Further investigations
including more efficient strategies for other coupling schemes and the
extension to realistic crystalline defects in three dimensions, such as
partial dislocations and grain boundaries, are also discussed at the end of
the paper.

###  1.1 Outline

The paper is structured as follows. In Section 2, we introduce the concept of
crystalline defects and the general formulation of the atomistic and the a/c
coupling models. Section 3 begins with a review of the a priori error
estimates for the blended a/c coupling methods previously presented in [15,
10], laying the foundation for the analysis in this work. We then provide a
rigorous residual-based a posteriori error estimate (Theorem 3.1) and
introduce the corresponding residual-force-based error estimator, which
provides an upper bound of the true approximation error (Theorem 3.2). In
Section 4, we propose an efficient evaluation of the local error contributions
and develop the corresponding adaptive algorithm for the blended a/c coupling
methods, in accordance with the theoretical results presented in the previous
section. In Section 5, we present numerical results for the adaptive
computations on the crystalline defects we consider and provide a
comprehensive discussion and explanation of our findings. We conclude by
summarizing our results and discussing possible future directions in Section
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$. For a finite set $A$, we will use $\\#A$ to denote the
cardinality of $A$. For second order tensors A and B, we denote
${\textsf{A}}:\textsf{B}=\sum_{i,j}{\textsf{A}}_{ij}\textsf{B}_{ij}$. The
closed ball with radius $r$ and center $x$ is denoted by $B_{r}(x)$, or
$B_{r}$ if the center is the origin.

##  2 The Atomistic Model and A/C Coupling Methods

In this section, we provide an exposition of the atomistic model and a general
formulation of a/c coupling methods. Several variations of a/c couplings [7,
15] can be accommodated within this form. Our motivation for presenting this
general form is to establish the principal contribution of our research: a
unified framework for a posteriori error analysis that can be extended to a
wide range of consistent a/c coupling methods. We defer the discussion of
specific a/c coupling methods to Section 4.

For the sake of simplicity, we consider the single-species Bravais lattices,
and we note that both the analysis and the algorithms discussed in this work
can be applied to multilattice crystals [34] with suitable minor
modifications. We keep the presentation as concise as possible since much of
the details can be found in various earlier works [7, 15, 10, 12, 20].

###  2.1 Atomistic model

Let ${\Lambda^{\rm hom}}={\textsf{A}}\mathbb{Z}^{2}$ with some non-singular
matrix ${\textsf{A}}\in\mathbb{R}^{2\times 2}$ be a perfect single lattice
possessing no defects and $\Lambda\subset\mathbb{R}^{2}$ be the corresponding
single lattice with certain local defects. The mismatch between $\Lambda$ and
${\Lambda^{\rm hom}}$ represents possible defects that are often contained in
some localized defect cores. To be precise, we assume that the defect is
contained in $\Omega^{\rm DEF}$ and $\Lambda\backslash\Omega^{\rm
DEF}={\Lambda^{\rm hom}}\backslash\Omega^{\rm DEF}$.

The deformed configuration of the infinite lattice $\Lambda$ is a map
$y\in\mathscr{U}$ with $\mathscr{U}:=\\{v:\Lambda\to\mathbb{R}^{2}\\}$, and it
can be decomposed as

$\displaystyle
y(\ell)=F\ell+u_{0}(\ell)+u(\ell),\qquad\forall~{}\ell\in\Lambda,$ (2.1)

where $F$ is a macroscopic deformation gradient, $u_{0}\in\mathscr{U}$ is a
far-field predictor resulting from the presence of the defect and
$u\in\mathscr{U}$ is a corrector. For point defects, we simply take $u_{0}=0$.
For anti-plane screw dislocation and anti-plane crack, $u_{0}$ can be derived
by solving a continuum linearized elasticity (CLE) equation [4] and they are
briefly reviewed in Appendix A.

We define the interaction neighbourhood for each $\ell\in\Lambda$ by
$\mathcal{N}_{\ell}:=\\{\ell^{\prime}\in\Lambda~{}|~{}0<|\ell^{\prime}-\ell|\leq
r_{\textrm{cut}}\\}$ with a given cut-off radius $r_{\textrm{cut}}$. We also
denote the interaction range
$\mathscr{R}_{\ell}:=\\{\ell^{\prime}-\ell~{}|~{}\ell^{\prime}\in\mathcal{N}_{\ell}\\}$.
For each atom $\ell\in\Lambda$, we define the finite difference stencil for
$u\in\mathscr{U}$

$\displaystyle
Du(\ell):=\\{D_{\rho}u(\ell)\\}_{\rho\in\mathscr{R}_{\ell}}:=\\{u(\ell+\rho)-u(\ell)\\}_{\rho\in\mathscr{R}_{\ell}}.$

To measure the local “regularity” of a displacement $u\in\mathscr{U}$, it is
convenient to use a background mesh, for example, a canonical triangulation
$\mathcal{T}_{\Lambda}$ of $\mathbb{R}^{2}$ into triangles whose nodes are the
reference lattice sites in $\Lambda$, see [4, Figure 1] for an illustration.
We define $I_{\textrm{a}}u$ as the standard piecewise affine interpolation of
$u$ with respect to $\mathcal{T}_{\Lambda}$. Notice that $\Lambda$ may
contains vacancies, and we can construct the interpolant with respect to
${\Lambda^{\rm hom}}$ by extending $u$ to vacancy sites [6, Appendix 1]. When
no confusion arises, we identify $u$ with interpolation $I_{\rm a}u$ and
denote the piecewise constant gradient $\nabla u=\nabla I_{\rm a}u$. We then
introduce the functional space of finite-energy displacements

$\displaystyle{\mathscr{U}}^{1,2}(\Lambda):=\big{\\{}u:\Lambda\rightarrow\mathbb{R}^{d}~{}\big{\lvert}~{}\|\nabla
u\|_{L^{2}}<\infty\big{\\}}$ (2.2)

with the associated norm $\|v\|_{{\mathscr{U}}^{1,2}}:=\|\nabla v\|_{L^{2}}$
for $v\in{\mathscr{U}}^{1,2}(\Lambda)$. We also define the following subspace
of compact displacements

$\displaystyle\mathscr{U}^{\rm
c}(\Lambda):=\big{\\{}u:\Lambda\rightarrow\mathbb{R}^{d}~{}\big{\lvert}~{}\textrm{supp}(\nabla
u)~{}~{}\textrm{is compact}\big{\\}}.$ (2.3)

It was shown in [7, Proposition 3.3] that $\mathscr{U}^{\rm c}$ is dense in
${\mathscr{U}}^{1,2}$.

The site potential is a collection of mappings
$V_{\ell}:(\mathbb{R}^{2})^{\mathscr{R}_{\ell}}\rightarrow\mathbb{R}$, which
represents the energy distributed to each atomic site. We refer to [7, Sectino
2] for a detailed discussion on the assumptions of general site potentials. We
can then formally define the energy-difference functional of the atomistic
model

$\displaystyle\mathcal{E}^{\textrm{a}}(u)$
$\displaystyle=~{}\sum_{\ell\in\Lambda}\Big{(}V_{\ell}\big{(}Du_{0}(\ell)+Du(\ell)\big{)}-V_{\ell}\big{(}Du_{0}(\ell)\big{)}\Big{)}=:\sum_{\ell\in\Lambda}V^{\prime}_{\ell}\big{(}Du(\ell)\big{)},$
(2.4)

where $V^{\prime}_{\ell}$ is the first renormalization of the potential
$V_{\ell}$ [10, Section 2.1].

The atomistic problem we would like to solve is formulated as

$u^{\rm a}\in\arg\min\big{\\{}\mathcal{E}^{\rm
a}(u)~{}|~{}u\in\mathscr{U}^{1,2}(\Lambda)\big{\\}},$ (2.5)

where “$\arg\min$” is understood as the set of local minimisers. It was shown
in [7, Lemma 1] that $\mathcal{E}^{\textrm{a}}(u)$ is well-defined on the
space ${\mathscr{U}}^{1,2}$, namely, the solutions to (2.5) exist under
suitable assumptions (cf. [7, Section 2.1]). A local minimizer $u$ satisfies
the following first and second order optimality conditions

$\displaystyle\big{\langle}\delta\mathcal{E}^{\textrm{a}}(u),v\big{\rangle}=0,\qquad\big{\langle}\delta^{2}\mathcal{E}^{\textrm{a}}(u)v,v\big{\rangle}\geq
0,\qquad\forall~{}v\in{\mathscr{U}}^{1,2}(\Lambda).$ (2.6)

In order to present the error analysis in Section 3, a stronger stability
assumption is required [2, 7]:

######  Assumption 2.1.

$\exists~{}\gamma>0$, such that for all $v\in{\mathscr{U}}^{1,2}(\Lambda)$,
$\big{\langle}\delta^{2}\mathcal{E}^{\rm
a}(u)v,v\big{\rangle}\geq\gamma\|v\|^{2}_{{\mathscr{U}}^{1,2}}$.

###  2.2 The general formulation of the a/c coupling methods

Since the atomistic problem is defined on an infinite domain and considers
every atom as a degree of freedom, it is practically unsolvable. To overcome
this challenge, and in line with the principle of the first order optimality
(2.6), the a/c coupling methods  [7, 15, 10] aim to solve an approximated
variational nonlinear system

$\displaystyle\langle\mathcal{F}_{h}^{\rm
ac}(u_{h}),v_{h}\rangle_{\textrm{ac}}=0,\quad\forall v_{h}\in\mathscr{U}_{h},$
(2.7)

where $\mathscr{U}_{h}$ denotes an approximation of the solution space of the
displacement ${\mathscr{U}}^{1,2}(\Lambda)$. The nonlinear map
$\mathcal{F}^{\rm ac}_{h}:\mathscr{U}_{h}\rightarrow(\mathscr{U}_{h})^{*}$
serves as an approximation of the first variation of the atomistic energy
$\delta\mathcal{E}^{\textrm{a}}(u)$. The duality
$\langle\cdot,\cdot\rangle_{\textrm{ac}}$ is induced by $\mathscr{U}_{h}$. It
is crucial to construct $\mathscr{U}_{h}$ and $\mathcal{F}^{\rm ac}_{h}$
effectively to ensure the well-posedness and solvability of (2.7) with
considerably fewer degrees of freedom.

The general formulation (2.7) fits most (if not all) of the a/c methods though
certain principles of approximation need to follow which we will explain
immediately. However, since our subsequent a posterior error estimates in
Section 3 do not depend on the individual formulation of a coupling scheme, we
keep the approximation principles as general as possible and concrete examples
will be specified later in Section 4.

The first approximation is the restriction of the computational domain to a
simply connected polygon $\Omega\subset\mathbb{R}^{d}$ and the error
introduced by this restriction is called the truncation error.

The second approximation is essentially a model reduction so that the nonlocal
atomistic model is approximated by a local continuum model. Recall that the
displacement is smooth in the region away from the defect core. A sufficiently
accurate approximated solution may be produced if we replace the atomistic
model by a continuum approximation in the far field. A typical choice is the
Cauchy-Born nonlinear elasticity model [35, 36], where the strain energy
function $W:\mathbb{R}^{2}\rightarrow\mathbb{R}^{2}$ measures the energy per
unit volume under a global homogeneous deformation, that is,

$W(\textsf{F}):=\det({\textsf{A}}^{-1})\cdot V(\textsf{F}\mathscr{R}).$ (2.8)

We define the first renormalization of $W$ as
$W^{\prime}(\textsf{F}):=W(\textsf{F})-W({\bf 0})$.

The third approximation is a domain decomposition with further coarse-graining
by the finite element discretization. In this step, we decompose the
computational domain $\Omega$ into the atomistic region $\Omega_{\textrm{a}}$,
the interface region $\Omega_{\rm i}$ and the continuum region
$\Omega_{\textrm{c}}$, so that $\Omega=\Omega_{\textrm{a}}\cup\Omega_{\rm
i}\cup\Omega_{\textrm{c}}\subset\mathbb{R}^{2}$. We assume that the defect
core is contained in the atomistic region, i.e., $\Omega^{\rm
DEF}\subset\Omega_{\textrm{a}}$, the displacement $u$ is sufficiently smooth
in $\Omega_{\textrm{c}}$ and $\Omega_{\rm i}$ is the region between $\Omega$
and $\Omega_{\textrm{c}}$ where the transition of different models happens. To
accommodate our presentation, we define the set of core atoms
$\Lambda_{\textrm{a}}:=\Lambda\cap\Omega_{\textrm{a}}$, the set of the
interface atoms $\Lambda_{\rm i}:=\Lambda\cap\Omega_{\rm i}$ and the set of
atoms in the continuum region $\Lambda_{\rm c}:=\Lambda\cap\Omega_{\rm i}$.
Let $\mathcal{T}_{\rm a}$ be the canonical triangulation induced by
$\Lambda_{\textrm{a}}\cup\Lambda_{\rm i}$ and $\mathcal{T}_{\rm c}$ be a
shape-regular simplicial partition of the continuum region. We can define
$\mathcal{T}_{h}=\mathcal{T}_{\textrm{a}}\cup\mathcal{T}_{\textrm{c}}$ as the
triangulation of the computational domain for the a/c coupling schemes. Let
$\Omega_{h}=\bigcup_{T\in\mathcal{T}_{h}}T$. We define the space of the
coarse-grained displacements by

$\displaystyle\mathscr{U}_{h}:=\big{\\{}u_{h}:\Omega_{h}\to\mathbb{R}^{2}~{}\big{|}~{}\text{$u_{h}\in\mathcal{P}_{1}(\mathcal{T}_{h})$,
$u_{h}=0$ in $\mathbb{R}^{d}\setminus\Omega_{h}$}\big{\\}},$ (2.9)

where a $\mathcal{P}_{1}$ finite element discretization is utilized in the
continuum region.

We illustrate the feasibility of the general formulation (2.7) and the
principles of approximation by the following prototypical examples.

The first example is the so-called energy-based a/c methods [12, 15, 18, 10].
For $u_{h}\in\mathscr{U}_{h}$, we first define the a/c coupling energy
functional

$\displaystyle\mathcal{E}^{\rm ac}_{h}(u_{h}):=$
$\displaystyle\sum_{\ell\in\Lambda_{\textrm{a}}}V^{\prime}_{\ell}(Du_{h})+\sum_{\ell\in\Lambda_{\rm
i}}V^{\rm
I}_{\ell}(Du_{h})+\sum_{T\in\mathcal{T}_{h}}\omega_{T}\big{(}W^{\prime}(\nabla
u_{h}+\nabla u_{0})-W^{\prime}(\nabla u_{0})\big{)},$ (2.10)

where $V^{\rm I}_{\ell}\in C^{k}((\mathbb{R}^{2})^{\mathcal{R}})$ is an
interface potential and $\omega_{T}$ is the effective volume of $T$ which
specify individual coupling schemes and are subtly defined to deal with the
spurious forces (commonly referred as “ghost force”) around the interface. In
this context, we have $\mathcal{F}_{h}^{\rm ac}:=\delta\mathcal{E}^{\rm
ac}_{h}$ in (2.7) which is the first variation of the coupling energy
functional $\mathcal{E}^{\rm ac}_{h}$ with respect to the space
$\mathscr{U}_{h}$.

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

$u^{\rm ac}_{h}\in\arg\min\big{\\{}\mathcal{E}^{\rm
ac}_{h}(u_{h}),~{}u_{h}\in\mathscr{U}_{h}\big{\\}}.$ (2.11)

Similar to the optimality conditions given in (2.6), the solution to the
energy minimization problem in (2.11) also satisfies the variational problem
in (2.7). However, it is important to note that the reverse is only true when
the corresponding coupling scheme is stable [19, 7].

The second example is the so-called force-based a/c coupling methods [16, 17,
32]. Instead of having a defined energy, the force-based methods directly
solve an equilibrium equation defined by

$\displaystyle\langle\mathcal{F}_{h}^{\rm
ac}(u_{h}),v_{h}\rangle_{\textrm{ac}}:=\langle\delta\mathcal{E}^{\rm
a}(u_{h}),(1-\beta)v_{h}\rangle+\langle\delta\mathcal{E}^{\rm
cb}_{h}(u_{h}),\beta v_{h}\rangle=0,$ (2.12)

where the function $\beta:\Omega\rightarrow\mathbb{R}$ specifies the force-
based coupling schemes and $\mathcal{E}^{\rm cb}_{h}$ is the Cauchy-Born
finite element functional whose concrete formulation is given in (4.55) in
Section 4. We note that the function $\beta$ often satisfies that $\beta=0$ in
$\Omega_{a}$ and $\beta=1$ in $\Omega_{c}$ but takes values between $0$ and
$1$ in $\Omega_{i}$ to generate the transition of the models. The force-based
coupling offers the notable benefit of eliminating spurious interface forces
commonly present in energy-based schemes. Nonetheless, this advantage is often
accompanied by a non-conservative force field, resulting in the inability to
conserve the moment in the simulations using molecular dynamics.

We note again that since our a posterior error estimates in Section 3 do not
depend on the formulation of a coupling scheme, we postpone the construction
detail of specific a/c schemes to Section 4 where adaptive algorithms, which
do depend on the corresponding coupling scheme, are developed. What
essentially makes our a posterior error estimates possible is the consistency
of the a/c method which we introduce in the following section.

###  2.3 The consistency of a/c coupling methods

Let $N$ denote the total number of degrees of freedom of the a/c coupling
method. We say that the method is consistent if the approximation error
decreases as $N$ increases. To be more precise, under additional assumptions
(e.g., blending function and mesh qualities [20, Assumption 1]), for $N$
sufficiently large, we have

$\displaystyle\|u-u^{\rm ac}_{h}\|_{{\mathscr{U}}^{1,2}}\leq CN^{-k}(\log
N)^{s},$ (2.13)

where $C$ is a constant independent of any model parameters, and $s,k\geq 0$
which are related to the types of defects and specify the convergence rate of
the approximation error with respect to $N$. The rate of convergence $k$ is
primarily determined by the regularity of the defect equilibrium and the
specific coupling schemes employed. We will provide a summary that outlines
the convergence rates for various cases considered in this work in Section 4.

It is worth mentioning that the a priori error estimate (2.13) provides a
solid theoretical foundation of the current work, which plays a key role in
proving the following residual based a posteriori error estimate (cf. Theorem
3.1).

##  3 The a Posteriori Error Estimates based on Residual Forces

In this section, we propose a unified framework of the a posteriori error
estimate that is generally applicable to any consistent atomistic-to-continuum
(a/c) coupling methods satisfying (2.7) which provides the groundwork for the
adaptive algorithms for a range of a/c coupling methods in the following
section. To be precise, we adapt the analytical framework originally
introduced by our previous research of adaptive QM/MM methods for simple point
defects [2] to the setting of a/c and extend it to more complicated defects.
Such framework can be summarized in four steps: the proof of the equivalence
of the true approximation error and the residual in certain dual norm, the
derivation of the a posteriori error estimator based on the residual forces
which gives upper bound of the residual, the estimate of the error committed
by the truncation of the computational domain for different types of defects
and an efficient sampling of the residual forces which forms the local error
contributions. We explain these four steps in detail in the following
subsections.

###  3.1 Equivalence of the residual and the approximation error

For the simplicity of presentation, we denote the solution to the atomistic
problem (2.5) as $u$ and that to the a/c problem (2.5) or (2.12) as $u_{h}$.
Recall that $\mathscr{U}^{\rm c}$ is dense in ${\mathscr{U}}^{1,2}$, we have
$I_{\textrm{a}}u_{h}\in{\mathscr{U}}^{1,2}$ where $I_{\textrm{a}}$ is the
standard piecewise affine interpolation operator introduced in Section 2.1.
The residual $\textsf{R}(I_{\rm a}u_{h})$ is then defined as an operator on
$\mathscr{U}^{1,2}$,

$\textsf{R}(I_{\rm
a}u_{h})[v]:=\langle\delta\mathcal{E}^{\textrm{a}}(I_{\textrm{a}}u_{h}),v\rangle,\quad\forall
v\in\mathscr{U}^{1,2}.$ (3.14)

The theorem below shows that the dual norm of the residual $\textsf{R}(I_{\rm
a}u_{h})$ gives both the upper and the lower bounds for the true approximation
error $\|u-I_{\rm a}u_{h}\|_{\mathscr{U}^{1,2}}$.

######  Theorem 3.1.

Let $u$ be a strongly stable solution of (2.5). Suppose that
$\delta\mathcal{E}^{\rm a}$ and $\delta^{2}\mathcal{E}^{\rm a}$ are Lipschitz
continuous in $B_{r}(u)$ with uniform constants $L_{1}$ and $L_{2}$, the a/c
coupling method is consistent in the sense of (2.13) and that $u_{h}$
satisfies the same decay as $u$. Then for $N$ sufficiently large, there exists
constants $c$ and $C$ independent of the approximation parameters such that

$\displaystyle c\|u-I_{\rm
a}u_{h}\|_{{\mathscr{U}}^{1,2}}\leq\|\mathsf{R}(I_{\rm
a}u_{h})\|_{({\mathscr{U}}^{1,2})^{*}}\leq C\|u-I_{\rm
a}u_{h}\|_{{\mathscr{U}}^{1,2}},$ (3.15)

where the residual $\mathsf{R}(I_{\rm a}u_{h})$ is defined by (3.14).

###### Proof.

Let $r>0$ satisfies $B_{r}(u)\subset{\mathscr{U}}^{1,2}$. By the Lipschitz
continuity of $\delta\mathcal{E}^{\textrm{a}}$ and
$\delta^{2}\mathcal{E}^{\textrm{a}}$ in $B_{r}(u)$ with uniform constants
$L_{1}$ and $L_{2}$, for any $w\in B_{r}(u)$, we have

$\displaystyle\|\delta\mathcal{E}^{\textrm{a}}(u)-\delta\mathcal{E}^{\textrm{a}}(w)\|\leq
L_{1}\|u-w\|_{{\mathscr{U}}^{1,2}},$
$\displaystyle\|\delta^{2}\mathcal{E}^{\textrm{a}}(u)-\delta^{2}\mathcal{E}^{\textrm{a}}(w)\|\leq
L_{2}\|u-w\|_{{\mathscr{U}}^{1,2}}.$

Using the fact that $\langle\delta\mathcal{E}^{\textrm{a}}(u),v\rangle=0$, we
can obtain for $N$ sufficiently large

$\textsf{R}(I_{\rm
a}u_{h})[v]=\langle\delta\mathcal{E}^{\textrm{a}}(I_{\textrm{a}}u_{h})-\delta\mathcal{E}^{\textrm{a}}(u),v\rangle\leq
L_{1}\|u-I_{\rm
a}u_{h}\|_{{\mathscr{U}}^{1,2}}\|v\|_{{\mathscr{U}}^{1,2}},\quad\forall
v\in\mathscr{U}^{1,2}.$

To estimate the term $\|u-I_{\rm a}u_{h}\|_{{\mathscr{U}}^{1,2}}$, from the a
priori estimate (2.13), we have

$\displaystyle\|u-I_{\rm a}u_{h}\|_{{\mathscr{U}}^{1,2}}$
$\displaystyle\leq\|u-u_{h}\|_{{\mathscr{U}}^{1,2}}+\|u_{h}-I_{\rm
a}u_{h}\|_{{\mathscr{U}}^{1,2}}\lesssim N^{-k}(\log N)^{s},$ (3.16)

where the last inequality follows from the assumption on $u_{h}$ that it
satisfies the same decay estimate as $u$.

The above result leads to the lower bound of the true approximation error

$\displaystyle\|\textsf{R}(I_{\rm a}u_{h})\|_{({\mathscr{U}}^{1,2})^{*}}\leq
C\|u-I_{\rm a}u_{h}\|_{{\mathscr{U}}^{1,2}},$ (3.17)

where the constant $C:=L_{1}$.

To get the upper bound, by the definition of $\textsf{R}(I_{\rm a}u_{h})$ in
(3.14) and the Galerkin orthogonality resulted from the first optimality
condition of the atomistic problem (2.6), it is straightforward to see that

$\displaystyle\|\textsf{R}(I_{\rm
a}u_{h})\|_{({\mathscr{U}}^{1,2})^{*}}\|u-I_{\rm
a}u_{h}\|_{{\mathscr{U}}^{1,2}}\geq\langle\delta\mathcal{E}^{\textrm{a}}(I_{\textrm{a}}u_{h}),u-I_{\textrm{a}}u_{h}\rangle=\langle\delta\mathcal{E}^{\textrm{a}}(I_{\textrm{a}}u_{h})-\delta\mathcal{E}^{\textrm{a}}(u),u-I_{\textrm{a}}u_{h}\rangle.$
(3.18)

A further application of the intermediate value theorem yields

$\displaystyle\langle\delta\mathcal{E}^{\textrm{a}}(I_{\textrm{a}}u_{h})-\delta\mathcal{E}^{\textrm{a}}(u),u-I_{\textrm{a}}u_{h}\rangle=\langle\delta^{2}\mathcal{E}^{\textrm{a}}(w)(u-I_{\textrm{a}}u_{h}),u-I_{\textrm{a}}u_{h}\rangle,$
(3.19)

where $w=tu+(1-t)I_{\textrm{a}}u_{h}$ for some $t\in(0,1)$.

To analyze the second variation of $\mathcal{E}^{\textrm{a}}$ at $w$, we add
and subtract the same term and separate it into two parts

$\displaystyle\langle\delta^{2}\mathcal{E}^{\textrm{a}}(w)v,v\rangle$
$\displaystyle=\langle\delta^{2}\mathcal{E}^{\textrm{a}}(w)v,v\rangle-\langle\delta^{2}\mathcal{E}^{\textrm{a}}(u)v,v\rangle+\langle\delta^{2}\mathcal{E}^{\textrm{a}}(u)v,v\rangle$
$\displaystyle=:S_{1}+S_{2}.$ (3.20)

For $N$ sufficiently large, e.g. $N\geq(\frac{\gamma}{2}L_{2})^{k}$, we can
bound $S_{1}$ by the Lipschitz continuity of
$\delta^{2}\mathcal{E}^{\textrm{a}}$ such that

$\displaystyle\big{|}S_{1}\big{|}\leq
L_{2}\|u-w\|_{{\mathscr{U}}^{1,2}}\|v\|^{2}_{{\mathscr{U}}^{1,2}}\leq\frac{\gamma}{2}\|v\|^{2}_{{\mathscr{U}}^{1,2}},$
(3.21)

where the second inequality holds by the fact that
$\|u-w\|_{{\mathscr{U}}^{1,2}}\leq\|u-I_{\textrm{a}}u_{h}\|_{{\mathscr{U}}^{1,2}}$
and (3.16). Since $u$ is a strongly stable solution of (2.5) satisfying
Assumption 2.1, we obtain the following stability estimate by (3.1) and (3.21)
that

$\displaystyle\langle\delta^{2}\mathcal{E}^{\textrm{a}}(w)v,v\rangle\geq\frac{\gamma}{2}\|v\|^{2}_{{\mathscr{U}}^{1,2}}.$
(3.22)

Let $v=u-I_{\rm a}u_{h}$ in (3.22), we have

$\displaystyle\langle\delta^{2}\mathcal{E}^{\textrm{a}}(w)(u-I_{\textrm{a}}u_{h}),u-I_{\textrm{a}}u_{h}\rangle\geq\frac{\gamma}{2}\|u-I_{\rm
a}u_{h}\|^{2}_{{\mathscr{U}}^{1,2}},$ (3.23)

Combining (3.18), (3.19) and (3.23) and dividing both sides by $\|u-I_{\rm
a}u_{h}\|_{{\mathscr{U}}^{1,2}}$ leads to the upper bound of the true
approximation error

$\displaystyle\|\textsf{R}(I_{\rm a}u_{h})\|_{({\mathscr{U}}^{1,2})^{*}}\geq
c\|u-I_{\rm a}u_{h}\|_{{\mathscr{U}}^{1,2}},$ (3.24)

where $c:=\gamma/2$. The stated results are obtained by (3.17) and (3.24). ∎

Theorem 3.1 establishes the equivalence of the residual $\textsf{R}(I_{\rm
a}u_{h})$ in the dual norm $\|\cdot\|_{({\mathscr{U}}^{1,2})^{*}}$ and true
approximation error $\|u-I_{\rm a}u_{h}\|_{\mathscr{U}^{1,2}}$. More
importantly, we note that such equivalence is independent of the specific
formulation of the a/c coupling schemes. The conditions on which the result
does rely are essentially the following three: the Lipschitz continuity of the
derivatives or the variations of $\mathcal{E}^{\textrm{a}}$, the consistency
of the coupling method defined by (2.13) and the assumption that $u_{h}$
shares the same decay estimates as $u$. The first condition is often
guaranteed by the regularity of the (empirical) interatomic potential (c.f.
[4, Lemma 2.1]). The second condition is the focus of the a priori analysis of
any a/c coupling method [7, 15]. The third condition is in fact not easy to
prove rigorously and is often proposed as an assumption [30]. One possible
approach to address this difficulty is to utilize the technique for
constructing the effective Green’s functions of a/c coupling methods which is
initially introduced in [4]. We also refer the interested reader to the
analysis in [5, Section 5.1] for further information. In the current work, we
numerically verify this assumption in Appendix C (cf. Figures A.1, A.2, and
A.3 for micro-crack, anti-plane screw dislocation, and anti-plane crack,
respectively).

From Theorem 3.1, we see that the ideal a posteriori error estimator is

$\eta^{\rm ideal}(u_{h}):=\|\textsf{R}(I_{\rm
a}u_{h})\|_{({\mathscr{U}}^{1,2})^{*}}.$

However, such quantity is not computable in practice since it is defined in a
dual norm. Hence, we derive a computable estimator in the following section,
which provides the upper bound of the true approximation error.

###  3.2 a Posteriori error estimator based on residual forces

Before we derive the computable a posteriori error estimator, we give a brief
review of the common practice in previous researches in the adaptive a/c
methods.

Since $\mathcal{E}^{\textrm{a}}(\cdot)$ is an energy functional on
${\mathscr{U}^{1,2}}$, the first variation of $\mathcal{E}^{\textrm{a}}$
naturally defines a stress-strain relation in the language of continuum
mechanics. For example, the direct calculation of the first variation of
$\mathcal{E}^{\textrm{a}}$ at $I_{\textrm{a}}u_{h}$ often gives

$\textsf{R}(I_{\textrm{a}}u_{h})[v]=\langle\delta\mathcal{E}^{\textrm{a}}(I_{\rm
a}u_{h}),v\rangle=\sum_{T\in\mathcal{T}_{\Lambda}}|T|\sigma^{\textrm{a}}(I_{\rm
a}u_{h}):\nabla_{T}v,$ (3.25)

where we can define $\sigma^{\textrm{a}}$ to be the atomistic stress tensor.
We refer to [36] for the detailed formulations of $\sigma^{\textrm{a}}$. One
way to give a quantitative and computable estimate of $\|\textsf{R}(I_{\rm
a}u_{h})\|_{({\mathscr{U}}^{1,2})^{*}}$ is to introduce the corresponding
stress-strain formulation of the a/c coupling method we are interested, which
is often obtained from (2.7), and add it to (3.25) so that the following
representation of the residual is obtained

$\textsf{R}(I_{\textrm{a}}u_{h})[v]=\sum_{T\in\mathcal{T}_{\Lambda}}|T|\sigma^{\textrm{a}}(I_{\rm
a}u_{h}):\nabla_{T}v-\sum_{T\in\mathcal{T}_{h}}|T|\sigma^{\textrm{ac}}(u_{h};T):\nabla_{T}v_{h}$
(3.26)

if $\mathcal{P}_{1}$ finite element discretization is used as we did in (2.9),
where $\sigma^{\textrm{ac}}$ is the coupling stress tensor and $v_{h}$ is a
certain interpolation of $v$ in $\mathscr{U}_{h}$. This often serves as the
starting point of the so-called residual-stress based a posteriori error
estimates for a/c methods. However, it is worth noting that there are two
major disadvantages brought by the introduction of the variational formulation
of the coupling scheme as given, for example, in (3.26). The first one is the
particularity of the a/c method so that the estimate of (3.26) should be
derived for each coupling scheme. The second one is the complicated
formulation of $\sigma^{\textrm{ac}}$ and the discrepancy between the test
function spaces ${\mathscr{U}}^{1,2}$ and $\mathscr{U}_{h}$ which result in
the involved subsequent manipulation. We refer to [6] for detailed examples
and a more recent work in [30] to illustrate that such complexity may
significantly bring down the efficiency of the adaptive algorithms.

On the other hand, the residual $\delta\mathcal{E}^{\textrm{a}}(I_{\rm
a}u_{h})$ can also be expressed in terms of the residual force for any
$v\in{\mathscr{U}}^{1,2}$ such that

$\textsf{R}(I_{\textrm{a}}u_{h})[v]=\langle\delta\mathcal{E}^{\textrm{a}}(I_{\rm
a}u_{h}),v\rangle=\sum_{\ell\in\Lambda}\mathscr{F}^{\rm a}_{\ell}(I_{\rm
a}u_{h})\cdot v(\ell),$ (3.27)

where the residual force is defined as $\mathscr{F}^{\rm a}_{\ell}(I_{\rm
a}u_{h}):=\partial_{u(\ell)}\mathcal{E}^{\textrm{a}}(u)\big{|}_{u={I_{\rm
a}u_{h}}}$. We note that such expression is again independent of the specific
formulation of the a/c coupling schemes and allows us to avoid the involved
stress-strain or the “weak” representations of the residual
$\textsf{R}(I_{\textrm{a}}u_{h})$. We thus adopt the residual-force based
representation of the residual given by (3.27) and consequently, the a
posteriori error estimate is given by the following lemma and theorem.

######  Lemma 3.1.

Let $d=2$. Then for any $\ell\in\Lambda$ there exists a constant $C^{\rm
d2}>0$ such that

$|v(\ell)-v(0)|\leq C^{\rm
d2}\log(2+|\ell|)\cdot\|v\|_{{\mathscr{U}}^{1,2}}\qquad\forall
v\in{\mathscr{U}}^{1,2}.$

Let $d=3$. Then there exists a constant $C^{\rm d3}>0$ such that, for each
$v\in{\mathscr{U}}^{1,2}$ there exists $v_{\infty}\in\mathbb{R}^{3}$ the
following estimate holds

$\|v-v_{\infty}\|_{\ell^{6}}\leq C^{\rm d3}\|v\|_{{\mathscr{U}}^{1,2}}.$

######  Theorem 3.2.

Let the residual forces $\mathscr{F}^{\rm a}_{\ell}(u_{h})$ be defined by
(3.27). Under the conditions of Theorem 3.1, there exists a constant $C^{\rm
resF}$ such that

$\displaystyle\|u-I_{\rm a}u_{h}\|_{\mathscr{U}^{1,2}}\leq
C\|\mathsf{R}(I_{\rm a}u_{h})\|_{(\mathscr{U}^{1,2})^{*}}\leq C^{\rm
resF}\eta(u_{h}),$ (3.28)

where

$\displaystyle\eta(u_{h}):=\left\\{\begin{array}[]{ll}\big{\|}\log(2+|\ell|)\cdot\mathscr{F}^{\rm
a}_{\ell}(I_{\rm a}u_{h})\big{\|}_{\ell^{1}},&d=2\\\\[4.30554pt]
\|\mathscr{F}^{\rm a}_{\ell}(I_{\rm
a}u_{h})\|_{\ell^{\frac{6}{5}}},&d=3\end{array}\right..$ (3.31)

We note that Lemma 3.1 serves as the Poincaré inequality for
${\mathscr{U}}^{1,2}$ and Theorem 3.2 gives an upper bound of
$\mathsf{R}(I_{\rm a}u_{h})$ using the residual force $\mathscr{F}^{\rm
a}_{\ell}(I_{\rm a}u_{h})$ on each lattice point $\ell\in\Lambda$. The proof
of Lemma 3.1 can be found in [37, Lemma A.2] while the proof of Theorem 3.2
follows the same line as that for a posteriori error estimate in QM/MM
coupling methods (cf. [2, Theorem 3.1]) and is thus omitted. However, a few
remarks should be made as below.

######  Remark 3.1 (Stress based vs. force based).

The residual-stress based estimate is often adopted in the a posteriori error
analysis for the finite element method (FEM) solving elliptic PDEs which, to
some extend, are similar to our models. There are two reasons for this. The
first is that the weak/variational formulation, which admits possible
singularity in the weak solution, is the first step of FEM as well as its
analysis. The one is that stress, rather than force, is the quantity related
to the geometric change in continuum mechanics which are modeled by elliptic
PDEs [38, 39]. In contrast, there are no singularities in the microscopic
model we are interested due to the discrete setting so that the strong and the
weak formulations, which are defined by (3.27) and (3.25) respectively, are
exactly the same (which may not be the case in PDEs). In addition, forces on
atoms are responsible for the geometric change in the microscopic scale [7,
40]. Therefore, the residual-force based estimate as presented in (3.2) is a
more natural choice for the problem we concern. However, we have to emphasize
that, by properly introducing additional terms, the stress based formulation
may help to separate the error into different parts, which are the modeling
error, the coarsening error and the truncation error respectively in the case
of a/c methods [6]. Such separation is crucial particularly in the a priori
analysis to obtain the convergence result of the methods [7, 4].

######  Remark 3.2 (Three dimensions).

Although the a posteriori error estimates are provided for both two and three
dimensions in Theorem 3.2, we focus on the implementation in two dimensions.
The adaptive a/c methods in three dimensions require a reliable three-
dimensional mesh generator and adaptation techniques, which is out of the
scope of the current work. However, we note that the framework and results can
be easily extended to three dimensions with corresponding modifications, e.g.
estimates of the truncation error which will be discussed in the next section.
We refer the interested readers to the recent work by Fu et al. [41] in this
direction.

######  Remark 3.3 (Lower bounds).

We note that the lower bound of either the ideal a posteriori error estimator
$\|\mathsf{R}(I_{\rm a}u_{h})\|_{(\mathscr{U}^{1,2})^{*}}$ or the true
approximation error $\|u-I_{\rm a}u_{h}\|_{\mathscr{U}^{1,2}}$ is yet to be
proved which is often technically involved. We refer to [31] in this direction
for the residual-stress based a posteriori error estimator and a more recent
progress in [5] on the Riesz representation of the residual force
$\mathscr{F}^{\rm a}_{\ell}(u_{h})$ in the context of QM/MM coupling methods
which, in principle, could be extended to the current work. However, the loss
of the lower bound does not have a significantly impact in our adaptive
simulations, as it primarily affects the efficiency of the error estimator.
Further discussions are made in Section 5.

The evaluation of $\eta(u_{h})$ requires the computation of the residual force
at each lattice point of the infinite lattice $\Lambda$ and is thus not
feasible in practice. In the subsequent sections, we develop an approximation
of $\eta(u_{h})$ so that it is only computed in the finite domain $\Omega$ and
is efficiently assembled and allocated to each element to guide the
adaptivity.

###  3.3 Estimates of the truncation error

In this section, we consider the approximation by the finite computational
domain.

We first define the coupling error estimator by

$\displaystyle\eta^{\textrm{ac}}(u_{h}):=\sum_{\ell\in\Lambda\cap\Omega}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm
a}_{\ell}(I_{\rm a}u_{h})\big{|},$ (3.32)

which is a weighted summation of the residual force $\mathscr{F}^{\rm
a}_{\ell}(I_{\rm a}u_{h})$ for $\ell\in\Lambda\cap\Omega$. To estimate the
error originated from the truncation of finite computational domain, we
introduce the truncation error estimator as follows. Recall that $r_{\rm cut}$
is the maximum radius of the interaction range introduced in Section 2.1, we
denote the extended domain of $\Omega$ as

$\Omega_{\rm ext}:=\bigcup_{\ell\in\Lambda\cap\Omega}B_{r_{\rm cut}+1}(\ell).$

The truncation error estimator is then defined by

$\displaystyle\rho^{\rm tr}(u_{h}):=\sum_{\ell\in\Lambda\cap(\Omega_{\rm
ext}\setminus\Omega)}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm
a}_{\ell}(I_{\rm a}u_{h})\big{|}.$ (3.33)

The following Lemma shows that the global error estimator $\eta(u_{h})$ can be
bounded by the global coupling error estimator $\eta^{\rm ac}(u_{h})$ and the
truncation indicator $\eta^{\rm tr}(u_{h})$.

######  Lemma 3.2.

Suppose that $\eta(u_{h})$ and $\eta^{\rm ac}(u_{h})$ are defined by (3.31)
and (3.32) respectively, then we have

$\eta(u_{h})\leq\eta^{\rm ac}(u_{h})+\eta^{\rm tr}(u_{h}),$

where the truncation indicator $\eta^{\rm tr}(u_{h})$ is defined by

$\displaystyle\eta^{\rm tr}(u_{h}):=\left\\{\begin{array}[]{ll}\rho^{\rm
tr}(u_{h}),&{\rm Point~{}defects}\\\\[4.30554pt] \rho^{\rm tr}(u_{h})+C^{\rm
disloc}\big{(}3+\log(R_{\rm ext})\big{)}\cdot R_{\rm ext}^{-1},&{\rm
Dislocations}\end{array}\right.,$ (3.36)

with $R_{\rm ext}:=R+r_{\rm cut}+1$, $\rho^{\rm tr}(u_{h})$ defined by (3.33)
and $C^{\rm disloc}$ a constant.

###### Proof.

Let $\Lambda_{\rm ext}:=\Lambda\cap\Omega_{\rm ext}$. From the definition of
$\eta(u_{h})$, $\eta^{\textrm{ac}}(u_{h})$ and $\rho^{\rm tr}(u_{h})$ by
(3.31), (3.32) and (3.33) respectively. We obtain

$\displaystyle\eta(u_{h})=\eta^{\textrm{ac}}(u_{h})+\rho^{\rm
tr}(u_{h})+\sum_{\ell\in\Lambda\setminus\Lambda_{\rm
ext}}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm a}_{\ell}({\bf 0})\big{|}.$
(3.37)

For the case of point defects, the far-field predictor is $u_{0}=0$. Hence,
for $\ell\in\Lambda\setminus\Lambda_{\rm est}$, we have $\mathscr{F}^{\rm
a}_{\ell}({\bf 0})=0$, which leads to

$\eta(u_{h})=\eta^{\textrm{ac}}(u_{h})+\rho^{\rm tr}(u_{h}).$

For general straight dislocations, the far-field predictor is not zero.
Therefore the above equation is no longer held. However, according to the
sharp decay estimates of the residual forces $|\mathscr{F}^{\rm a}_{\ell}({\bf
0})|\leq C^{\rm disloc}\cdot|\ell|^{-3}$ (c.f. [4, Lemma 5.8]), we have

$\displaystyle\sum_{\ell\in\Lambda\setminus\Lambda_{\rm
ext}}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm a}_{\ell}({\bf 0})\big{|}$
$\displaystyle\leq~{}C^{\rm disloc}\int^{\infty}_{R_{\rm ext}}\log(2+r)\cdot
r^{-2}\,{\textrm{d}}r$ $\displaystyle\leq~{}C^{\rm disloc}\big{(}3+\log(R_{\rm
ext})\big{)}\cdot R_{\rm ext}^{-1}.$

Hence, taking into account this with (3.37) yields the stated results for
dislocations. ∎

######  Remark 3.4.

The sharp decay estimate of the residual force for cracks has not been
rigorously established yet. However, it is reasonable to expect that
$|\mathscr{F}^{\rm a}_{\ell}({\bf 0})|\leq C^{\rm crack}\cdot|\ell|^{-2.5}$ by
exploiting the fact that the decay of the far-field predictor for crack is
$|Du_{0}(\ell)|\lesssim|\ell|^{-0.5}$ whereas
$|Du_{0}(\ell)|\lesssim|\ell|^{-1}$ for straight dislocation. This speculation
will be numerically verified in Figure 10. Hence, for anti-plane crack we
consider in this work, the truncation indicator can be defined by

$\eta^{\rm tr}(u_{h}):=\rho^{\rm tr}(u_{h})+C^{\rm crack}\big{(}3+\log(R_{\rm
ext})\big{)}\cdot R_{\rm ext}^{-0.5}.$

The rigorous proof of the elastic far-field behavior for cracks is currently
not available and beyond the scope of current work. Readers who are interested
in this direction may refer to the analysis presented in [42, 43].

The truncation indicator $\rho^{\rm tr}(u_{h})$ serves as a measure of the
error that results from the finite size of the computational domain. The
constants $C^{\rm disloc}$ and $C^{\rm crack}$ in the definition of $\rho^{\rm
tr}(u_{h})$ can be empirically determined by fitting the residual forces, as
shown in Figure 5 and Figure 10 for anti-plane screw dislocation and anti-
plane crack, respectively. Therefore, the truncation indicator $\rho^{\rm
tr}(u_{h})$ can be computed for all types of defects studied in this work, and
its value is checked against a threshold during the main adaptive procedure.
Further details will be discussed in Section 4.2, Algorithm 2.

###  3.4 Local error contribution by sampling

The coupling a posteriori error estimator $\eta^{\textrm{ac}}(u_{h})$ is
defined globally. In order to determine the optimal locations for refining the
mesh in the continuum region or expanding the atomistic/blended regions, it is
necessary to distribute the global residual estimator
$\eta^{\textrm{ac}}(u_{h})$ into element-wise local contributions.

Recall that $\mathcal{T}_{h}$ is the finite element partition of the
computational domain $\Omega$ defined in Section 2.2, intuitively we can
rewrite $\eta^{\textrm{ac}}(u_{h})$ in terms of the element
$T\in\mathcal{T}_{h}$ as

$\displaystyle\eta^{\textrm{ac}}(u_{h})$
$\displaystyle=\sum_{\ell\in\Lambda\cap\Omega}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm
a}_{\ell}(I_{\rm a}u_{h})\big{|}$
$\displaystyle=~{}\sum_{T\in\mathcal{T}_{h}}\sum_{\ell\in
T}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm a}_{\ell}(I_{\rm a}u_{h})\big{|}$
$\displaystyle=:\sum_{T\in\mathcal{T}_{h}}\eta^{\textrm{ac}}_{T}(u_{h}).$
(3.38)

It is important to note that although the derivation of the elementwise error
estimator (3.4) is conceptually straightforward, it can be computationally
expensive. In particular, evaluating (3.4) requires determining the atoms in
$T$ for each $T\in\mathcal{T}_{h}$, which require a careful consideration of
the geometric relationship between $\ell$ and $T$. The approximate
computational cost for this process would be estimated as ${\rm DoF}^{2}$,
where ${\rm DoF}$ denotes the total number of degrees of freedom of a/c
coupling methods. See our recent work [30] for more detailed discussions.
Moreover, since $u_{h}$ is smooth in regions far from the defect core, the
change residual force $\mathscr{F}^{\rm a}_{\ell}(I_{\rm a}u_{h})$ is also
relatively smooth in such region. Consequently, there is no need to evaluate
the residual force for every atom $\ell\in T$ for $T\in\mathcal{T}_{h}$ that
is far from the defect core during adaptive computations.

Motivated by the sampling strategy given in [2, Section 3.3], we propose the
following approximation of the coupling error estimator
$\eta^{\textrm{ac}}(u_{h})$ by

$\displaystyle\tilde{\eta}^{\textrm{ac}}(u_{h})$
$\displaystyle:=\sum_{T\in\mathcal{T}_{h}}\omega(T)\log(2+|\tilde{{\ell}}(T)|)\cdot\big{|}\mathscr{F}^{\rm
a}_{\tilde{{\ell}}(T)}(I_{\rm
a}u_{h})\big{|}=:\sum_{T\in\mathcal{T}_{h}}\tilde{\eta}^{\textrm{ac}}_{T}(u_{h}),$
(3.39)

with the corresponding approximated elementwise coupling error estimator

$\displaystyle\tilde{\eta}^{\textrm{ac}}_{T}(u_{h}):=\omega(T)\log(2+|\tilde{\ell}(T)|)\cdot\big{|}\mathscr{F}^{\rm
a}_{\tilde{\ell}(T)}(I_{\rm a}u_{h})\big{|},$ (3.40)

where $\tilde{\ell}(T)$ is the “repatom” of $T$ and $w(T)$ gives the weight of
the element $T\in\mathcal{T}_{h}$ (e.g., one can take $w(T)$ to be the area of
$T$). Multiple repatoms within one element $T$ is also allowed. For example,
$\tilde{\ell}(T)$ can be the set of nodes of element $T$ with
$\omega(T):=|T|/3$ for triangular mesh in two dimensions. We will use this
strategy throughout our numerical experiments in Section 5.

The following theorem shows that the approximated estimator (3.40) provides an
accurate approximation, given that the residual force is adequately smooth
when the atomistic region is sufficiently large. The proof is mainly dependent
on the interpolation error analysis in relation to the definition provided in
equation (3.40), which will be given in Appendix B.

######  Theorem 3.3.

Let $\tilde{\eta}^{\textrm{ac}}(u_{h})$ and $\eta^{\textrm{ac}}(u_{h})$ be
defined by (3.4) and (3.39), respectively. Denote $R_{\Omega}$ as the radius
of $\Omega$ and let $\widetilde{\mathscr{F}}^{\rm a}$ be a
$C^{2}$-interpolation of $\mathscr{F}^{\rm a}$ in $\Omega$. For sufficiently
large $R_{\rm a}$, we have

$\big{|}\tilde{\eta}^{\textrm{ac}}(u_{h})-\eta^{\textrm{ac}}(u_{h})\big{|}\lesssim\log(R_{\Omega})\cdot\|\nabla^{2}\widetilde{\mathscr{F}}^{\rm
a}(I_{\rm a}u_{h})\|_{L^{2}(\Omega)}.$ (3.41)

The evaluation of (3.39) incurs a significantly lower computational cost as it
avoids the need to determine geometric relationships and requires only a few
evaluation of the residual forces which is proportion to the number of
elements. This approach is similar to that proposed in our recent work [30].
We will provide numerical verification of the computational cost reduction in
Section 5.

##  4 Adaptive Algorithms for Various A/C Methods

In this section, we develop the adaptive algorithms for the a posteriori error
control problem for various a/c coupling methods based on our residual-force
based error estimator. We concentrate on two major classes of a/c methods,
namely the a/c methods with sharp interface and the blended type a/c methods,
which are the methods most often encountered. We note that such classification
is somehow different from the one we indicated in Section 2 where the methods
are categorized as energy-based and force-based. The reason for this
difference is that the energy vs. force-based categorization stems from a
modeling perspective while the sharp interface vs. blended categorization is
more suitable in view of the adaptive algorithms that will be developed
immediately. As we mentioned previously, because of the unified framework for
the a posteriori error estimates, only algorithmic differences exist for
different a/c schemes. In addition, a unified adaptive algorithm can be
developed for each class of methods which further simplify the problem.

For each class of a/c methods, we will first briefly introduce the detailed
formulations of typical coupling schemes and then develop the corresponding
adaptive algorithms. Potential applications of our unified a posteriori error
estimates to other classes of a/c methods are also commented at the end of
this section.

###  4.1 The a/c methods with sharp interface and the adaptive algorithm

We first consider the a/c methods with sharp interface and corresponding
adaptive algorithm. The character of this class of methods is to have a very
narrow transition region from the atomistic to the continuum regions (often a
few layers of atoms depending on the interaction range of the system we
consider), which is achieved either by introducing an interface energy
functional (often identified by $V^{\rm I}$ in (2.10)) or using a sharp
characteristic function instead of a blending function (often identified by a
$0-1$ function $\beta$ in (2.12)) for energy-based and force-based schemes
respectively. We also note that most of the existing literature on the a
posteriori error control and adaptivity of the a/c methods focus on such class
of methods, which provides us a good point to illustrate the advantage of our
unified framework.

####  4.1.1 The formulations of the a/c methods with sharp interface

Here we introduce two consistent a/c methods with sharp interface.

GRAC method

The geometric reconstruction based consistent a/c (GRAC) coupling method was
initially proposed in [18] and urther developed in [44]. The GRAC method match
exactly the formulation as shown in (2.10), that is,

$\displaystyle\mathcal{E}^{\rm grac}_{h}(u_{h}):=$
$\displaystyle\sum_{\ell\in\Lambda_{\textrm{a}}}V^{\prime}_{\ell}(Du_{h})+\sum_{\ell\in\Lambda_{\textrm{i}}}\omega_{\ell}V^{\textrm{i}}_{\ell}(Du_{h})+\sum_{T\in\mathcal{T}_{h}}\omega_{T}W(\nabla
y_{h}|_{T}),$ (4.42)

where $V_{\ell}^{\textrm{i}}$ is a modified interface site potential and
$\omega_{\ell}$ and $\omega_{T}$ are certain coefficients which are usually
called the effective volumes of lattice sites or elements (see [18, Section
2.2] for a detailed discussion). The essential idea for GRAC method is to
introduce the geometric reconstruction parameters $C_{\ell;\rho,\varsigma}$ so
that for each
$\ell\in\Lambda^{\textrm{i}},\rho,\varsigma\in\mathscr{R}_{\ell}$ the
interface potential is redefined by

$V_{\ell}^{\textrm{i}}(Du_{h}):=V\Big{(}\big{(}{\textstyle\sum_{\varsigma\in\mathscr{R}_{\ell}}C_{\ell;\rho,\varsigma}D_{\varsigma}u_{h}(\ell)}\big{)}_{\rho\in\mathscr{R}_{\ell}}\Big{)},$
(4.43)

and parameters are then determined by solving the following equations:

$V_{\ell}^{\textrm{i}}({\bf 0})=V({\bf
0}),\quad\forall\ell\in\Lambda_{\textrm{i}}\quad\text{and}\quad\langle\delta\mathcal{E}^{\rm
grac}_{h}({\bf 0}),v\rangle=0\quad\forall v\in\mathscr{U}_{h},$ (4.44)

which are termed as the energy and force patch test consistency [45].

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

$u^{\rm
grac}_{h}\in\arg\min\big{\\{}\mathcal{E}_{h}(u_{h}),~{}u_{h}\in\mathscr{U}_{h}\big{\\}}.$
(4.45)

QCF method

The QCF method is the force-based a/c coupling method without blending, that
is,

$\displaystyle\langle\mathcal{F}_{h}^{\rm
qcf}(u_{h}),v_{h}\rangle:=\langle\delta\mathcal{E}^{\rm
a}(u_{h}),\chi_{\Lambda_{\textrm{a}}}v_{h}\rangle+\langle\delta\mathcal{E}^{\rm
cb}_{h}(u_{h}),(1-\chi_{\Lambda_{\textrm{a}}})v_{h}\rangle,$ (4.46)

where $\chi_{\Lambda_{\textrm{a}}}$ is the characteristic function of the
atomistic region $\Lambda_{\textrm{a}}$.

In the QCF method we solve the following variational nonlinear system

$\displaystyle\langle\mathscr{F}_{h}^{\rm qcf}(u^{\rm
qcf}_{h}),v_{h}\rangle=0,\quad\forall v_{h}\in\mathscr{U}_{h}.$ (4.47)

####  4.1.2 Consistency of the a/c methods with sharp interface

The consistency of the a/c methods with sharp interface including the above
mentioned GRAC and QCF methods have been studied in detail in the previous
works [18, 7]. Here we only give a brief summary. Under the assumptions on
mesh qualities [20, Assumption 1]), for $N$ sufficiently large, we have

$\displaystyle\|u-u^{\rm grac,qcf}_{h}\|_{{\mathscr{U}}^{1,2}}\leq
C\left\\{\begin{array}[]{lll}N^{-1},&\textrm{Point defects}\\\\[4.30554pt]
N^{-1}(\log N)^{0.5},&\textrm{Dislocations}\\\\[4.30554pt] N^{-0.25}(\log
N)^{0.5},&\textrm{Cracks}\end{array}\right..$ (4.51)

Note that the rigorous a priori error estimates of the a/c methods with sharp
interface for dislocations and cracks are still lacking. The corresponding
rates shown above are speculated based on the methodologies employed in [4],
which are beyond the scope of this work but deserve a future study.

####  4.1.3 Adaptive algorithms for the a/c methods with sharp interface

Existing literatures on the a posteriori error estimate for sharp interface
a/c methods rely on the stress-based approach which we briefly touched at the
beginning of Section 3.2. The disadvantage of such approach is that the a
posteriori error estimator should be derived for each method and the detailed
formulation of the estimator may be involved [6, 32, 29] which may increase
the complexity of the implementation. Moreover, it was first pointed out in
[29] that the computational cost of the residual-stress based a posteriori
error estimator for the GRAC method is even larger than that of solving the
a/c problem. It was later shown that the inefficiency may be improved by a
modified a posteriori error estimator together with an “interface buffer”
region which lies in several layers of atoms around the a/c interface where
the finite element mesh coincides with the crystal lattice. However, such
improvement is only demonstrated for systems with nearest-neighbor
interactions and the “interface buffer” region may be large for the system
with finite range interactions such that unnecessary degrees of freedoms may
be introduced to raise the computational cost again.

To avoid the difficulty just mentioned, we adopt the residual-force based
error estimator just derived to develop our adaptive algorithm given in
Algorithm 1 which can be applied to general sharp interface a/c methods. We
note that the algorithm essentially follows the one previously presented in
[6, 30] with only the replacement of the error estimators. The accuracy and
computational efficiency are later demonstrated in Section 5.

Algorithm 1 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}$, $K$ and $R_{\max}$.

  2. Step 1

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

  3. Step 2

Estimate: Compute the local error estimator $\tilde{\eta}^{\textrm{ac}}_{T}$
by (3.40) for each $T\in\mathcal{T}_{h}$, and the truncation indicator
$\rho^{\rm tr}$ by (3.36). Compute the degrees of freedom $N$ and
$\tilde{\eta}^{\textrm{ac}}:=\sum_{T}\tilde{\eta}^{\textrm{ac}}_{T}$. If
$\rho^{\rm tr}>\tau_{1}\tilde{\eta}^{\textrm{ac}}$, enlarge the computational
domain. Stop if $N>N_{\max}$ or $\tilde{\eta}^{\textrm{ac}}<\eta_{\rm tol}$.

  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}}\tilde{\eta}^{\textrm{ac}}_{T}\geq\frac{1}{2}\tilde{\eta}^{\textrm{ac}}.$

    2. Step 3.2

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

$\sum_{T\in\mathcal{M}^{k}_{\textrm{i}}}\tilde{\eta}^{\textrm{ac}}_{T}\geq\tau_{2}\sum_{T\in\mathcal{M}}\tilde{\eta}^{\textrm{ac}}_{T}.$
(4.52)

Let $\mathcal{M}=\mathcal{M}\setminus\mathcal{M}^{k}_{\textrm{i}}$. Expand the
interface $\Lambda^{\textrm{i}}$ outward by $k$ layers.

  5. Step 4

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

###  4.2 The blended a/c methods and their adaptive algorithms

We then consider the blended a/c coupling methods and the corresponding
adaptive algorithm. The primary advantage of the blended methods is their
simplicity of construction and implementation, which make these methods
attractive for the simulations of complex crystalline defects in real-world
applications.

####  4.2.1 The formulations of the blended a/c methods

All the blended a/c methods depend on a so-called blending function $\beta\in
C^{2,1}(\mathbb{R}^{d})$ satisfying $\beta=0$ in $\Omega_{\textrm{a}}$,
$\beta=1$ in $\Omega_{\textrm{c}}$ and ${\rm
supp}(\nabla\beta)\subset\Omega_{\rm b}$ which characterizes the transition of
the atomistic model in the defect core to the continuum model in the far
field. To present the detailed formulations, we also define the piecewise
constant mid-point interpolant $Q_{h}v\in\mathcal{P}_{0}(\mathcal{T}_{h})$ for
a function $v:\Omega_{h}\rightarrow\mathbb{R}$ such that $Q_{h}v(x):=v(x_{T})$
for $x\in T\in\mathcal{T}_{h}$, where $x_{T}:=\frac{1}{|T|}\int_{T}x\rm dx$.
We then introduce the detailed formulations of three typical blended a/c
methods.

BQCE method

The BQCE method is an energy-based method [12, 15]. For
$u_{h}\in\mathscr{U}_{h}$, the BQCE energy functional is obtained by combining
atomistic and continuum energy functionals via a blending function $\beta$,
that is,

$\mathcal{E}_{h}^{\rm
bqce}(u_{h}):=\sum_{\ell\in\Lambda\cap\Omega_{h}}\big{(}1-\beta(\ell)\big{)}V^{\prime}_{\ell}(Du_{h})+\int_{\Omega_{h}}Q_{h}\beta(x)\big{(}W^{\prime}(\nabla
u_{h}+\nabla u_{0})-W^{\prime}(\nabla u_{0})\big{)}\rm dx,$ (4.53)

where $V^{\prime}_{\ell}$ and $W^{\prime}$ are defined by (2.4) and (2.8)
respectively. The BQCE problem we would like to solve is to find

$u^{\rm bqce}_{h}\in\arg\min\big{\\{}\mathcal{E}^{\rm
bqce}_{h}(u_{h}),~{}u_{h}\in\mathscr{U}_{h}\big{\\}}.$ (4.54)

BQCF method

While the BQCE method (4.53) blends the atomistic and the continuum energies,
the BQCF method [15, 46] blends atomistic and continuum forces and is thus a
force-based method. We first define the pure Cauchy-Born finite element
functional [15], for $u_{h}\in\mathscr{U}_{h}$,

$\displaystyle\mathcal{E}^{\rm
cb}_{h}(u_{h}):=\int_{\Omega_{h}}Q_{h}\big{[}W^{\prime}(\nabla u_{h}+\nabla
u_{0})-W^{\prime}(\nabla u_{0})\big{]}\rm dx.$ (4.55)

Recall that $\beta\in C^{2,1}(\mathbb{R}^{2})$ is a blending function, then
the BQCF operator is the nonlinear map $\mathscr{F}^{\rm
bqcf}_{h}:\mathscr{U}_{h}\rightarrow(\mathscr{U}_{h})^{*}$, defined by

$\displaystyle\langle\mathscr{F}^{\rm
bqcf}_{h}(u_{h}),v_{h}\rangle:=\langle\delta\mathcal{E}^{\rm
a}(u_{h}),(1-\beta)v_{h}\rangle+\langle\delta\mathcal{E}^{\rm
cb}_{h}(u_{h}),\beta v_{h}\rangle.$ (4.56)

In the BQCF method we solve the following variational nonlinear system

$\displaystyle\langle\mathscr{F}_{h}^{\rm bqcf}(u^{\rm
bqcf}_{h}),v_{h}\rangle=0,\quad\forall v_{h}\in\mathscr{U}_{h}.$ (4.57)

BGFC method

The BGFC method is first proposed and analyzed in [10] combining the benefits
of blending [15] and ghost force correction [47], which aims to achieve the
optimal rate of convergence in terms of the degrees of freedom. The detailed
derivation of BGFC method is given in [10] and here we introduce an
alternative but simplified approach considering the “ghost force removal
formulation” [47] to derive the an equivalent formulation.

Let $\mathcal{E}^{\rm bqce}_{\rm hom}(\hat{u}_{0})$ be the BQCE energy
functional defined on the homogeneous lattice at a suitable “predictor”
$\hat{u}_{0}$, then the energy funtional of BGFC method can be written as

$\displaystyle\mathcal{E}^{\rm bgfc}(u_{h})=\mathcal{E}^{\rm
bqce}(u_{h})-\big{\langle}\delta\mathcal{E}^{\rm bqce}_{\rm
hom}(\hat{u}_{0}),u_{h}\big{\rangle}.$ (4.58)

In our implementation of the BGFC method, motivated by [10, Section 2.7 and
Section 4.2], we use (4.58) with $\hat{u}_{0}=0$ for all types of crystalline
defects considered in this paper. It is worth mentioning that the choice of
$\hat{u}_{0}$ could have a significant impact on the accuracy of the method,
especially in applications involving dislocations and cracks. We plan to
investigate this alternative point of view rigorously in future work.

The BGFC problem we would like to solve is to find

$u^{\rm bgfc}_{h}\in\arg\min\big{\\{}\mathcal{E}^{\rm
bgfc}_{h}(u_{h}),~{}u_{h}\in\mathscr{U}_{h}\big{\\}}.$ (4.59)

####  4.2.2 Consistency of the blended a/c methods

The selection of the parameter $\beta$ significantly impacts the accuracy of
the error estimates for the corresponding blended a/c coupling methods.
Interested readers may refer to previous works such as [7, 12, 15] for a
comprehensive analysis and discussion on this topic. In this work, we employ
the same blending function developed and analyzed in [10, 12], where $\beta$
is obtained in a preprocessing step by approximately minimizing
$\|\nabla^{2}\beta\|_{L^{2}}$ [10].

Angolous to the a/c methods with sharp interface, the consistency of the
blended a/c methods for point defects can also be found in [7]. The
consistency results for anti-plane dislocation essentially follows from the
earlier works [10, 15] but are adapted to the setting of the current work,
while the error estimate for anti-plane crack is a new speculated result that
will be analyzed rigorously in our future work.

Under the assumptions on blending function and mesh qualities [20, Assumption
1]), for $N$ sufficiently large, we have

$\displaystyle\|u-u^{\rm bqce}_{h}\|_{{\mathscr{U}}^{1,2}}\leq$ $\displaystyle
C\left\\{\begin{array}[]{lll}N^{-0.5},&\textrm{Point defects}\\\\[4.30554pt]
N^{-1}(\log N)^{0.5},&\textrm{Dislocations}\\\\[4.30554pt] N^{-0.25}(\log
N)^{0.5},&\textrm{Cracks}\end{array}\right..$ (4.63) $\displaystyle\|u-u^{\rm
bqcf,bgfc}_{h}\|_{{\mathscr{U}}^{1,2}}\leq$ $\displaystyle
C\left\\{\begin{array}[]{lll}N^{-1},&\textrm{Point defects}\\\\[4.30554pt]
N^{-1}(\log N)^{0.5},&\textrm{Dislocations}\\\\[4.30554pt] N^{-0.25}(\log
N)^{0.5},&\textrm{Cracks}\end{array}\right..$ (4.67)

It is worthwhile mentioning that to obtain the above mentioned convergence
rates, we exploit the a priori knowledge that $R_{a}=R_{b}$, i.e., the radius
of the atomistic region is the same as the width of the blending region. In
Section 5 we will observe the same relationship on the fly during the adaptive
procedure for all types of crystalline defects considered in this paper.

####  4.2.3 Adaptive algorithms for the blended a/c methods

It could be easily anticipated from, e.g. (4.53), (4.55) and (4.58), that the
residual-stress based a posteriori error estimates for the blended a/c methods
are even more involved than those for the sharp interface ones. Indeed, the
only existing work for blended type of methods using the residual-stress based
approach is [32] where the a posteriori error estimator for BQCF method in 1D
is derived which already has a complicated expression (c.f. [32, Theorem
3.2]). It is also worth noting that the adaptive algorithm developed there
uses the a priori result that $R_{a}\approx R_{b}$. However, since the
blending width plays a crucial role in the performance of the blending type
methods, a fully adaptive algorithm should dynamically determine the length of
the blending region as well as expand the computational domain, enlarge the
atomistic region and refine the elements in the continuum region.

We note that Algorithm 2 is the main adaptive algorithm for the blended a/c
coupling methods which follows a similar line of Algorithm 1 and Algorithm 3
is in charge of the marking step which is explained in detail below.

Algorithm 2 Adaptive blended a/c coupling based on residual forces

  1. Step 0

Prescible $\Omega$, $\mathcal{T}_{h}$, $N_{\max}$, $\eta_{\rm tol}$, $K$,
$\tau_{1}$ and $\tau_{2}$.

  2. Step 1

Solve: Solve the BQCE (BQCF or BGFC) solution $u_{\textrm{h}}$ of (4.54)
((4.57) or (4.59)) on the current mesh $\mathcal{T}_{h}$.

  3. Step 2

Estimate: Compute the local error estimator $\tilde{\eta}^{\textrm{ac}}_{T}$
by (3.40) for each $T\in\mathcal{T}_{h}$, and the truncation indicator
$\rho^{\rm tr}$ by (3.36). Compute the degrees of freedom $N$ and
$\tilde{\eta}^{\textrm{ac}}:=\sum_{T}\tilde{\eta}^{\textrm{ac}}_{T}$. If
$\rho^{\rm tr}>\tau_{1}\tilde{\eta}^{\textrm{ac}}$, enlarge the computational
domain. Stop if $N>N_{\max}$ or $\tilde{\eta}^{\textrm{ac}}<\eta_{\rm tol}$.

  4. Step 3

Mark: Apply Algorithm 3 to construct the refinement set $\mathcal{M}$, the
total number of atomistic layers to be expanded for both atomistic and
blending regions $k$ and the ratio $\alpha$ representing the error
contribution of atomistic region.

  5. Step 4

Refine: Expand the atomistic region outward by $[\alpha k]$ layers and the
blending region outward by $k-[\alpha k]$ layers. Bisect all elements
$T\in\mathcal{M}$. Go to Step 1.

Algorithm 3 Mark step

  1. Step 1

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

$\sum_{T\in\mathcal{M}}\tilde{\eta}^{\textrm{ac}}_{T}\geq\frac{1}{2}\tilde{\eta}^{\textrm{ac}}.$

  2. Step 2

: We can find the interface elements which are within $k$ layers of atomistic
distance, $\mathcal{M}^{k}:=\\{T\in\mathcal{M}\bigcap(\mathcal{T}_{\rm
b}\cup\mathcal{T}_{\textrm{c}}):{\rm{dist}}(T,\Lambda_{\textrm{a}})\leq k\\}$.
Find the first $k\leq K$ such that

$\sum_{T\in\mathcal{M}^{k}}\tilde{\eta}^{\textrm{ac}}_{T}\geq\tau_{2}\sum_{T\in\mathcal{M}}\tilde{\eta}^{\textrm{ac}}_{T}.$
(4.68)

  3. Step 3

: Construct the set $\mathcal{M}^{k}_{\textrm{a}}:=\\{T\in\mathcal{M}^{k}:{\rm
dist}(T,\Lambda_{\textrm{a}})\leq{\rm dist}(T,\Lambda_{\textrm{c}})\\}$.
Compute the ratio

$\displaystyle\alpha:=\frac{\sum_{T\in\mathcal{M}^{k}_{\textrm{a}}}\tilde{\eta}^{\textrm{ac}}_{T}}{\sum_{T\in\mathcal{M}^{k}}\tilde{\eta}^{\textrm{ac}}_{T}}.$
(4.69)

Let $\mathcal{M}:=\mathcal{M}\setminus\mathcal{M}^{k}$.

The distance between $T$ and $\Lambda_{\textrm{a}}$ is defined as ${\rm
dist}(T,\Lambda_{\textrm{a}}):=\inf\\{|\ell-
x_{T}|,\forall\ell\in\Lambda_{\textrm{a}}\\}$ and the barycenter of $T$ is
defined to be $x_{T}$ . Similarly we can define ${\rm
dist}(T,\Lambda_{\textrm{c}})$. The Mark step mainly applies the standard
Dörfler strategy [48] to choose a set $\mathcal{M}$ for model refinement. More
precisely, we construct a set $\mathcal{M}^{k}\subset\mathcal{M}$ such that it
contains the marked elements that are within $k$ layers of atomistic distance
to $\Lambda_{\rm a}$. This $k$ layers are the total layers of atomistic and
blending regions to be expanded. The set $\mathcal{M}\setminus\mathcal{M}^{k}$
contains the elements to be refined in continuum region. Next, we design a
strategy demonstrating how to enlarge $\Lambda_{\textrm{a}}$ and $\Lambda_{\rm
b}$ based on the error contributions in $\mathcal{M}^{k}$ to these two
regions. For $T\in\mathcal{M}^{k}$, if ${\rm
dist}(T,\Lambda_{\textrm{a}})\leq{\rm dist}(T,\Lambda_{\textrm{c}})$, we
assign it to constitute a new subset $\mathcal{M}^{k}_{\textrm{a}}$. We
compute the ratio of the summation of the error estimator for
$T\in\mathcal{M}^{k}_{\textrm{a}}$ and that for $T\in\mathcal{M}^{k}$ (cf.
(4.69)). According to this ratio $\alpha$, we extend the atomistic and
blending regions outwards by $[\alpha k]$ and $k-[\alpha k]$ layers
respectively.

###  4.3 Discussion and possible extension to other a/c coupling methods

We would like to highlight that, apart from the blended and shape-interface
coupling methods discussed in the previous sections, there are various other
a/c coupling methods that are used in both engineering and mathematical
communities. To demonstrate the versatility of the error estimator proposed in
this work, we present two examples as follows.

The first one is the so-called flexible boundary condition (FBC) method [49,
50], which applies the continuum solutions as the boundary conditions of the
atomistic problem in the core region. The FBC method can be formulated as:
find $u^{\rm fbc}_{h}:=\\{u^{\rm a},u^{\rm c}\\}$ such that

$\displaystyle({\bf P}^{\rm a})\;\left\\{\begin{aligned} \;\mathcal{L}[u^{\rm
a}]&=0&\qquad&\text{in}\;\Lambda_{\rm a},\\\ \;u&=u^{\rm
c}&&\text{in}\;\Lambda_{\rm i},\end{aligned}\right.$ $\displaystyle({\bf
P}^{\rm c})\;\left\\{\begin{aligned} \;\mathcal{L}_{\rm cb}[u^{\rm
c}]&=0&\qquad&\text{in}\;\Lambda_{\rm c},\\\ \;u&=u^{\rm
a}&&\text{on}\;\Lambda_{\rm i}.\end{aligned}\right.,$ (4.70)

where $\mathcal{L}$ and $\mathcal{L}_{\rm cb}$ are the force operator of
atomistic and Cauchy-Born models, respectively. It is shown in [49] that the
consistency of the FBC method is the same as that of the a/c methods with
sharp interface (cf. (4.51)).

The second one is the optimization-based a/c methods, which require the
atomistic and continuum subdomains with an overlap region
$\Omega_{\textrm{i}}:=\Omega_{\textrm{a}}\cap\Omega_{\textrm{c}}$. This is an
alternative ghost-force-free method. The optimization-based a/c methods are to
solve the constrained minimization problem: find $u^{\rm opt}_{h}:=\\{u^{\rm
a},u^{\rm c}\\}$ such that $\|\nabla u^{\textrm{a}}-\nabla
I_{\textrm{a}}u^{\textrm{c}}\|_{L^{2}(\Omega_{\rm i})}$ is minimized subject
to

$\displaystyle\left\\{\begin{array}[]{ll}\langle\delta\mathcal{E}^{\textrm{a}}(u^{\textrm{a}}),v\rangle=0&\forall
v\in\mathscr{U}^{\textrm{a}}_{0}\\\\[4.30554pt]
\langle\delta\mathcal{E}^{\textrm{c}}(u^{\textrm{c}}),v\rangle=0&\forall
v\in\mathscr{U}^{\textrm{c}}_{0}\end{array}\right.\quad\textrm{and}\quad\int_{\Omega_{\textrm{i}}}u^{\textrm{a}}-I_{\textrm{a}}u^{\textrm{c}}\rm
dx=0.$ (4.73)

The objective ensures that the mismatch between $\bar{u}^{\textrm{a}}$ and
$\bar{u}^{\textrm{c}}$ over $\Omega_{\textrm{i}}$ is as small as possible.
Again, the consistency of this method is the same as that of a/c methods with
sharp interface [51].

Both of these aforementioned methods hold substantial significance within the
engineering community, yet a noticeable gap exists in the availability of
corresponding adaptive algorithms. In light of the framework shown in Section
3, the proposed residual force-based error estimator can in principle be
applied to both two methods. This primarily attributed to the fact that it is
independent from the detailed formulations of the a/c methods. Notably, within
the framework of the FBC method, the residual force-based error estimator can
be interpreted as an extended speed for an interface motion problem. This
perspective enables its resolution through the utilization of the well-
established fast marching method [52]. Conversely, in the case of
optimization-based adaptive/corrective methods, the integration of the
adaptive algorithm necessitates its incorporation within the interface
minimization problem. A thorough exploration of this promising direction will
be a focal point of our future work.

##  5 Numerical Experiments

In this section, we perform adaptive simulations for three prototypical types
of defects: micro-crack, anti-plane screw dislocation, and anti-plane crack,
using our adaptive algorithms developed in Algorithm 2. The geometries of the
defects are illustrated in Figure 1.

(a) Micro-crack

(b) Anti-plane screw dislocation

(c) Anti-plane crack

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

We focus on the implementation for a two dimensional triangular lattice
${\Lambda^{\rm hom}}=\mathsf{A}\mathbb{Z}^{2}$ defined by

${\sf A}=\left(\begin{array}[]{cc}1&\cos(\pi/3)\\\
0&\sin(\pi/3)\end{array}\right),$

where ${\Lambda^{\rm hom}}$ is in fact the projection of a BCC crystal along
the (111) direction upon rotating and possibly dilating [4].

We adopt the well-known EAM potential as the site potential $V_{\ell}$ all
through our simulation 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\mathscr{R}_{\ell}}\phi\big{(}|D_{\rho}y(\ell)|\big{)}+F\Big{(}{\textstyle\sum_{\rho\in\mathscr{R}_{\ell}}}\psi\big{(}|D_{\rho}y(\ell)|\big{)}\Big{)},$
(5.74)

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)$
$F(\tilde{\rho})=C\left[(\tilde{\rho}-\tilde{\rho_{0}})^{2}+(\tilde{\rho}-\tilde{\rho_{0}})^{4}\right]$

with the parameters $a=4,b=3,c=10$ and $\tilde{\rho_{0}}=6\exp(-0.9b)$, which
are the same as the numerical experiments presented in [6, 29, 30]. The radius
of the computational domain $\Omega$ is set to be 300 initially and the
adaptive processes start with the initial configuration with $R_{\rm a}=3$.
The adaptive parameters in Algorithm 1 and Algorithm 2 are fixed to be
$\tau_{1}=1.0$ and $\tau_{2}=0.7$.

###  5.1 Micro-crack

The first defect we consider is the micro-crack, which is a prototypical
example of point defects which 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
simulation, we set $k=11$ and we apply an isotropic stretch $\mathrm{S}$ and
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$.

####  5.1.1 Adaptive a/c methods with sharp interface for micro-crack

In this section, we specifically consider the nearest neighbor setting for the
GRAC method.This choice is made due to the increased complexity and potential
instability associated with generalizing the method to finite interaction
ranges, as discussed later in this section. Moreover, by focusing on the
nearest neighbor setting, we can directly compare our current work with our
previous study conducted in Section 5 as presented in [30]. The adaptive
results for the QCF method are expected to exhibit qualitative similarity to
those of the GRAC method. Consequently, we omit them here for the sake of
simplicity of presentation.

According to [18, Proposition 3.7], under the setting of the current work, for
$\ell\in\Lambda^{\textrm{i}}$, the parameters $C_{\ell;\rho,\zeta}$ in (4.43)
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^{\textrm{i}}$; (d) $C_{\ell,\rho,\rho}=1$ for
$\ell+\rho\in\Lambda^{\textrm{a}}$; (e) $C_{\ell,\rho,\rho}=2/3$ for
$\ell+\rho\in\Lambda^{\textrm{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 [44].

We first consider the adaptive simulations of the micro-crack by adaptive a/c
methods with sharp interface. In particular, we test the adaptive GRAC method
based on three different a posteriori error estimators which are the original
residual-stress based error estimator developed in [6], the modified residual-
stress based error estimator developed in [30] and the residual-force based
error estimator given in (3.4). The corresponding adaptive algorithms are [6,
Algorithm 3], [30, Algorithm 2] and Algorithm 1 in the current work. We
restrict ourselves to a system with nearest neighbor interactions since the
construction and the implementation of the GRAC method and its adaptivity for
finite range interactions are very much involved. We will comment on this
issue later at the end of this section.

(a) Convergence

(b) Efficiency factors

(c) Estimation time

Figure 2: The convergence of the error, the CPU time of evaluating the
estimators and solving the a/c problem and the efficiency factors of the
estimators for adaptive GRAC method for the micro-crack.

Figure 2a shows the convergence of the true error
$\|u-I_{\textrm{a}}u_{h}\|_{\mathscr{U}^{1,2}}$ with respect to the number of
degrees of freedom (DoF). For the purpose of comparison, we also plot the
convergence of the error based on an a priori graded mesh. It is clearly seen
that all three adaptive algorithms achieve the almost the same optimal
convergence rate and the lines of convergence are barely distinguishable.

Figure 2b shows the efficiency factors (which is defined to be the ratios of
the error estimators and the true error) of the three different error
estimators. As we expect, the residual-force based error estimator developed
in the current work only provides an upper bound of the true error and the
efficiency factor of which is slightly inferior to those of the stress-based
estimators developed in [6] and [30]. However, the over estimate still
moderate since we consider a two dimensional problem and it is observed that
the residual-force based estimator possesses certain asymptotic exactness that
deserves a further investigation.

Figure 2c shows the CPU time of evaluating the a posteriori error estimators.
As we pointed out at the beginning of Section 3.2 and in [30] in more detail,
the original residual-stress based a posteriori error estimator is highly
inefficient as the time of evaluating the error estimator exceeds that of
solving the problem itself. On the contrary, the costs of evaluations of the
modified residual-stress based error estimator and the residual-force based
error estimator are comparable and are marginal compared with those for
original residual-stress based error estimator, which indicates that both
estimators may be adopted in practice.

We have to comment that the true advantages of the residual-force based error
estimator are not demonstrated by the above figures. They lie in the following
two aspects. First, the significant reduction of complexity in the
implementation. The development and computation of the residual-force based
error estimator introduced in this study exhibit a notably higher degree of
simplicity compared to its counterpart based on residual-stress. As
demonstrated in [6], the residual-stress based error estimator necessitates a
laborious computation involving the so-called “stress tensor correction”
technique. This requirement stems from the inherent lack of uniqueness in
stress formulations. Additionally, in order to further mitigate computational
expenses, a manual selection of an ”interface buffer” region is required [30].
Second, the possible extension to the systems with finite range interactions.
The residual-stress based error estimator for the GRAC method with finite
range interactions has been discussed in [29]. However, its formulation
necessitates significant approximations. Notably, when finite range
interactions are taken into account, the atomistic stress loses its locality,
prompting the introduction of ad hoc approximations to render the residual
amenable to element-wise summation. It is important to emphasize that these
approximations are effective only for simple point defects; for more intricate
defects such as dislocations and cracks, which are the focus of our present
investigation, this approach becomes untenable. In contrast, the residual-
force based error estimator proposed in this study offers a versatile
solution, readily extendable to finite range interactions due to its
straightforward construction. The adaptive simulation for systems of finite
range interactions are given in the numerical experiments for blended a/c
methods immediately in the next section.

####  5.1.2 Adaptive blended a/c methods for micro-crack

We then conduct the adaptive simulations of the micro-crack by blended a/c
methods. In particular, we test the adaptive BQCE, BQCF and BGFC methods which
share the same adaptive algorithms proposed in Algorithm 2 and 3. We
incorporate nearest and next-nearest neighbor interactions in our simulations
to demonstrate the capability of our a posteriori error estimate and adaptive
algorithms for dealing with systems with finite range interactions which is
one of the significant contributions of this study.

The implementations of all three blended a/c coupling (BQCE, BQCF and BGFC)
methods considered in this study remain consistent with those detailed in a
prior work [10]. Specifically, the blending function is derived during a
preprocessing stage by seeking an approximate minimization of
$\|\nabla^{2}\beta\|_{L^{2}}$, as elaborated upon comprehensively in [12]. For
the BGFC method, we employ the equivalent ”ghost force removal formulation”
(4.58), selecting $\hat{u}_{0}=0$ as the predictor for the sake of simplicity.

(a) Convergence

(b) Efficiency factors

Figure 3: The convergences of the error and the efficiency factors for both
rigorous and approximated error estimators with respect to the number of
degrees of freedom for the micro-crack.

In Figure 3a, we present the convergence of the true error
$\|u-I_{\textrm{a}}u_{h}\|_{\mathscr{U}^{1,2}}$ for the adaptive BQCE, BQCF,
and BGFC methods. Again we observe that all adaptive computations achieve the
theoretically verified optimal convergence rate ($N^{-0.5}$ for BQCE [12] and
$N^{-1.0}$ for both BQCF [46] and BGFC [10] respectively) compared with the a
priori results (cf. (4.63)) generated by a graded mesh.

Figure 3b displays the efficiency factors the residual-force based error
estimator for the three adaptive blended a/c methods. The constant prefactors
in the estimators overestimate the true approximation error at the beginning
of the adaptive simulations but become moderate asymptotically which have
similar behavior as that for the adaptive a/c method with sharp interface.

Figure 3b presents the efficiency factors for the rigorous and approximated
error estimators. As discussed in Remark 3.3, the a posteriori error estimator
based on a weighted $\ell^{2}$-norm of residual force is challenging to ensure
a lower bound on the true approximation error. The efficiency factors for all
three coupling methods shown in Figure 3b indicate that the constant
prefactors in the estimators overestimate the true approximation error
significantly at the beginning of the adaptive computations. Subsequently, the
prefactors become more moderate with slight oscillation, which has a
negligible effect on the adaptivity as the number of degrees of freedom
increases.

(a) CPU times

(b) $R_{\textrm{a}}/R_{\rm b}$

Figure 4: The CPU times for each steps (left) and the ratio between the radius
of the atomistic region $R_{\textrm{a}}$ and the width of the blending region
$R_{\rm b}$ (right) for three blended a/c coupling methods in the adaptive
computations for the micro-crack.

Figure 4a illustrates the CPU times required for computing the error
estimators and solving the corresponding blended a/c coupling problems,
plotted against the number of degrees of freedom $N$. Notably, it is evident
that the CPU times for error estimator computation exhibit a linear scaling,
denoted as $O(N)$, while the CPU times for solving all three blended a/c
methods demonstrate an approximately quadratic scaling behavior, described as
$O(N^{2})$. Furthermore, we observe that solving the BQCF coupling method is
more computationally expensive than solving the BQCE or BGFC methods with the
same number of degrees of freedom since the nonlinear system needs to be
solved for the force-based coupling scheme.

We plot the ratio between $R_{\textrm{a}}$ and $R_{\rm b}$ in the adaptive
computations for the micro-crack in Figure 4b, where $R_{\textrm{a}}$ is the
radius of the atomistic region while $R_{\rm b}$ represent the width of the
blending region. We observe that our adaptive algorithm can achieve nearly
optimal computational balance, i.e., $R_{\textrm{a}}\approx R_{\rm b}$, which
also verifies the a priori assumption on the relationship presented in [15,
Eq. (19)]). Furthermore, the ratio is consistently greater than one,
indicating that more atoms are allocated to the atomistic region during
adaptive computations. This observation substantiates the effectiveness of our
adaptive algorithm, as it consistently strives to reduce errors in the
direction of increased atomistic resolution.

###  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 [4], we restrict
the discussion and implementation to anti-plane shear motion which is
illustrated in Figure 1b. In this case, the Burgers vector is $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 there is no
additional shear deformation applied. The unknown for the anti-plane model is
the displacement in $e_{3}$-direction. The derivation of the far-field
predictor $u_{0}$ for anti-plane screw dislocation is reviewed in the Appendix
A.1.

Replicating the setting from [4], we use a simplified EAM-type interatomic
potential for the anti-plane case, given by

$V_{\ell}(y):=G\Big{(}\sum_{\rho\in\mathscr{R}_{\ell}}\phi\big{(}y(\ell+\rho)-y(\ell)\big{)}\Big{)},$

where

$G(s)=1+0.5s^{2}\quad\textrm{and}\quad\phi(r)=\sin^{2}(\pi r).$

Note that in this case, the BQCE and BGFC methods (cf. (4.53) and (4.58)) are
in fact identical since $\delta\mathcal{E}^{\rm bqce}_{\rm hom}(0)=0$ in
(4.58), which is an artefact of the anti-plane setting [10].

Figure 5 plots the decay of the residual force
$\mathcal{F}^{\textrm{a}}_{\ell}({\boldsymbol{0}})$ for the anti-plane screw
dislocation, which verifies the sharp estimates in [4, Lemma 5.8]. The black
dashed line serves as a validation of our selection of the constant ($C^{\rm
disloc}\approx 0.30$) utilized in defining the truncation indicator by
equation (3.36). This choice is justified by the observation that all the
atoms (sites) exhibit residual forces below this line. This represents a
pivotal distinction from the scenarios involving point defects, where the
residual forces are precisely zero. We leverage this distinction to construct
the truncation indicator, enhancing the robustness and accuracy of the
corresponding adaptive algorithm.

Figure 5: Decay of $\mathcal{F}^{\textrm{a}}_{\ell}({\boldsymbol{0}})$ for the
anti-plane screw dislocation.

####  5.2.1 Adaptive a/c methods with sharp interface for screw dislocation

In this section, we consider the adaptive simulations of screw dislocation by
adaptive a/c methods with sharp interface. Similarly as presented in Section
5.1.1, we consider the nearest neighbor setting for the GRAC method by testing
three different a posteriori error estimators.

(a) Convergence

(b) Efficiency factors

(c) Estimation time

Figure 6: The convergence of the error, the CPU time of evaluating the
estimators and solving the a/c problem and the efficiency factors of the
estimators for adaptive GRAC method for screw dislocation.

Much like the results for the micro-crack case presented in Section 5.1.1,
Figure 6a presents the convergence behavior of the true error
$\|u-I_{\textrm{a}}u_{h}\|_{\mathscr{U}^{1,2}}$ concerning the number of
degrees of freedom. Remarkably, all three adaptive algorithms exhibit nearly
identical optimal convergence rates, with their convergence curves being
indistinguishable from one another. Figure 6b offers insights into the
efficiency factors. The efficiency factor of the residual-force based error
estimator displays a marginal underestimation in comparison to the stress-
based estimators. However, it is crucial to emphasize that this slight
underestimation remains well within acceptable limits. Figure 6c presents the
CPU time necessary for evaluating the a posteriori error estimators. The
computational costs associated with assessing the modified residual-stress
based error estimator and the residual-force based error estimator are closely
aligned and are significantly lower than those associated with the original
residual-stress based error estimator.

####  5.2.2 Adaptive blended a/c methods for screw dislocation

We then conduct the adaptive simulations of the screw dislocation by blended
a/c methods. First of all, the corresponding computational mesh and the
atomistic region used in the construction of the blended a/c coupling methods
are illustrated in Figure 7.

Figure 7: The illustration of the computational mesh and the atomistic region
as used in the construction of the blended a/c coupling methods for anti-plane
screw dislocation.

We observe in Figure 8a that all adaptive computations achieve the same
optimal convergence rate compared with the a priori graded mesh given by
(4.63). We note again that the BGFC method is identity to the BQCE method due
to the anti-plane setting. Figure 8b plots the efficiency factors for both
rigorous and approximated error estimators for the blended a/c coupling
methods, which are still moderate.

(a) Convergence

(b) Efficiency factors

Figure 8: 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. In this case BGFC is identity to BQCE
due to the artefact of the anti-plane setting.

Figure 9a plots the total CPU times versus the numbers of degrees of freedom
$N$ for both rigorous and approximated error estimators. Figure 9b visualizes
the relationship between $R_{\textrm{a}}$ and $R_{\rm b}$ during the adaptive
computations for the anti-plane screw dislocation. Similar results are
observed for the anti-plane screw dislocation as those in the case of micro-
crack, which demonstrates not only the accuracy but also the efficiency of the
proposed residual-force error estimator. Our adaptive algorithm consistently
attains the optimal relationship between $R_{\textrm{a}}$ and $R_{\rm b}$,
thus reinforcing its effectiveness even in the context of anti-plane screw
dislocation.

(a) CPU times

(b) $R_{\textrm{a}}/R_{\rm b}$

Figure 9: The CPU times for each steps (left) and the ratio between the radius
of the atomistic region $R_{\textrm{a}}$ and the width of the blending region
$R_{\rm b}$ (right) for three blended a/c coupling methods in the adaptive
computations for the anti-plane screw dislocation.

###  5.3 Anti-plane crack

The last type of defect we consider is the case of anti-plane crack, which is
of much more practical interest than the previous two cases conducted above
(cf. the micro-crack in Section 5.1 and the anti-plane screw dislocation in
Section 5.2). It has not been considered in the literature of adaptive a/c
coupling methods to the best knowledge of the authors. We apply the same anti-
plane setting as that for screw dislocation presented in the last example. The
only difference is that the far-field predictor $u_{0}$ for anti-plane crack
is more complicated, which is briefly reviewed in the Appendix A.2.

Figure 10 shows the decay of the residual force
$\mathcal{F}^{\textrm{a}}_{\ell}({\boldsymbol{0}})$ for the anti-plane crack.
We observe that the residual force of the surface atoms lying around crack tip
only decays like $|\ell|^{-1.5}$ and the others decay as $|\ell|^{-2.5}$. The
black dashed and the black dotted lines validate our choice of the constants
($C^{\rm crack}_{\rm surf}$ for the surface atoms and $C^{\rm crack}_{\rm
oth}$ for others) as the residual force of all the atoms (sites) are under
these two lines. According to the discussions given in Remark 3.4, we choose
$C^{\rm crack}=C^{\rm crack}_{\rm surf}+C^{\rm crack}_{\rm oth}=3.35$ in our
adaptive computations.

Figure 10: Decay of $\mathcal{F}^{\textrm{a}}_{\ell}({\boldsymbol{0}})$ for
the anti-plane crack.

The significance and motivation behind the study of cracks, distinct from
micro-cracks and screw dislocations discussed in previous sections, encompass
three key aspects. Firstly, the implementation of the GRAC method for cracks
presents considerable challenges and appears impractical. In accordance with
the fundamental principles of the GRAC method [18], the determination of
parameters $C_{\ell;\rho,\zeta}$ in (4.43) around the interface for cracks is
hindered by the presence of crack tips. Consequently, the blended a/c coupling
method appears to be a more suitable choice for simulating cracks. Secondly,
conducting rigorous a priori analysis of a/c coupling methods proves
considerably more challenging. Within the framework proposed in [4], the
critical components of a priori analysis entail the establishment of rigorous
equilibrium decay estimates. This encompasses an intricate analysis of the
lattice Green’s function, along with the theoretical proof of residual-force
estimates, as numerically illustrated in Figure 11 above. However, when it
comes to a posteriori error estimates, we can somewhat circumvent these
intricate technical nuances, highlighting the novelty and significance of our
proposed methodology, especially if it can be effectively applied to crack
simulations. Thirdly, constructing a residual-stress based error estimator for
cracks remains ambiguous and seemingly implausible. This is due to the fact
that such construction relies heavily on the geometry of the defect and the
formulation of specific a/c coupling methods. In light of these
considerations, we proceed with numerical tests employing the proposed
residual-force based error estimator for all three blended a/c coupling
methods.

The corresponding computational mesh and the atomistic region used in the
construction of the blended a/c coupling methods are then illustrated in
Figure 11. Here we need to point out that though the same anti-plane setting
is applied as that for anti-plane screw dislocation presented in the last
section, the BQCE and BGFC methods are not identical any more for anti-plane
crack since a row of mesh around crack tip is removed as shown in Figure 11.

Figure 11: The illustration of the computational mesh and the atomistic region
as used in the construction of the blended a/c coupling methods for anti-plane
crack. A row of mesh around crack tip is removed to simulate the anti-plane
crack system [42].

Figure 3a plots that all adaptive computations achieve the same optimal
convergence rate ($N^{-0.25}$ for all three blended methods) compared with the
a priori results (cf. (4.63)) generated by a graded mesh. As we discussed, the
rigorous a priori error estimate for this case deserves further investigation.

(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 anti-plane crack.

Figure 12b plots the efficiency factors for the blended a/c coupling methods.
Figure 13a plots the total CPU times versus the numbers of degrees of freedom
$N$ for both rigorous and approximated error estimators. Similar results are
observed for the anti-plane crack as those in the case of micro-crack and
anti-plane screw dislocation, which demonstrates not only the accuracy but
also the efficiency of the residual-force based error estimator.

(a) CPU times

(b) $R_{\textrm{a}}/R_{\rm b}$

Figure 13: The CPU times for each steps (left) and the ratio between the
radius of the atomistic region $R_{\textrm{a}}$ and the width of the blending
region $R_{\rm b}$ (right) for three blended a/c coupling methods in the
adaptive computations for the anti-plane crack.

In Figure 13b, we depict the correlation between $R_{\textrm{a}}$ and $R_{\rm
b}$ for the adaptive computations performed on the anti-plane crack. Still,
our adaptive algorithm successfully achieves the optimal relationship between
$R_{\textrm{a}}$ and $R_{\rm b}$ for the anti-plane crack scenario.

##  6 Conclusion

In this work we propose a unified framework of the residual-based a posteriori
error estimates and design the corresponding adaptive algorithms that are
essentially applicable to any consistent multiscale coupling methods. We prove
that the error estimator based on the residual forces can provide the upper
bound of the true approximation error. As prototypical examples, we present a
range of adaptive computations based on this reliable error estimator for the
blended atomistic-to-continuum (a/c) coupling methods including the energy-
based blended quasi-continuum (BQCE), the force-based blended quasi-continuum
(BQCF) and the recently developed blended ghost force correction (BGFC)
methods. We develop coarse-grained techniques to efficiently evaluate the
error estimator and test them with different types of crystalline defects,
some of which have not been previously considered in related literature on
adaptive a/c coupling methods. The various numerical results demonstrate that,
compared to a priori error estimates, the adaptive algorithm provides the same
optimal convergence rate of the error with significant computational
efficiency. We note that the techniques and strategies presented in this study
do not require a specific choice of the multiscale coupling scheme as long as
the method is consistent. The present investigation offers valuable insights
into the development and application of adaptive multiscale techniques, and
constitutes a noteworthy addition to the existing body of literature on
atomistic-to-continuum coupling methods.

Although we believe that the adaptive algorithm proposed in this paper is
generally applicable for other common multiscale coupling schemes and more
complex crystalline defects, this research still raise a few open problems
which deserve further mathematical analysis and algorithmic developments.

  * •

A priori error estimate of BGFC method for straight edge dislocation and anti-
plane crack: The a priori results for edge dislocations and anti-plane cracks
are not currently available. To address this issue, a possible solution is to
investigate an equivalent ghost force removal formulation (4.58), where a
suitable “predictor” $\hat{u}_{0}$ needs to be constructed. In future work, we
plan to explore this alternative approach thoroughly, particularly with
regards to selecting nontrivial $\hat{u}_{0}$ in applications that involve
cracks and dislocations.

  * •

Three dimensional model problems: It is anticipated that the current work will
be extended to three dimensions in future. This, however, presents significant
challenges in mesh generation and adaptation in implementation. In comparison
to the two-dimensional case, the difficulties in three dimensions are
threefold: firstly, the atomistic region is not guaranteed to be convex;
secondly, a smooth transition region may be required for complex defects with
large distortions; and thirdly, a robust and efficient mesh adaptation is
essential. Recent advancements in this area [41] should provide a solid
foundation for further research.

  * •

More complex crystalline defects: Practical crystalline defects such as
partial dislocations connected by a stacking fault, dislocation nucleation,
and grain boundaries have already attracted significant attention. However,
constructing and implementing appropriate a/c coupling methods for these
problems remains a major challenge. Rigorous a priori analysis for such
problems, including proper boundary conditions and complex interface geometry,
is very difficult if possible. However, these are precisely the scenarios
where adaptive a/c coupling methods are expected to shine, demonstrating their
greatest advantage and potential. We believe that the approach presented in
this research provides a foundation for developing efficient and robust
adaptive algorithms for these practical and important problems.

Both the theoretical and practical aspects discussed above will be explored in
future work.

##  Appendix Appendix A Far-field Predictors

###  Appendix A.1 Dislocations

We model dislocations by following the setting used in [4]. 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}.$ (A.1)

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 [53, 4] 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}\\},$ (A.2)
$\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 [4, § 7] for more detail).

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

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

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$.

###  Appendix A.2 Cracks

We present the setting of cracks by following [54], which stems from the
limitation of the continuum elasticity approaches to static crack problems.
Similar with the discussions of dislocations, we introduce the following CLE

$\displaystyle-{\rm div}~{}(\mathbb{C}:\nabla u)=0\qquad$
$\displaystyle\text{in}~{}~{}\mathbb{R}^{2}\setminus\Gamma,$
$\displaystyle(\mathbb{C}:\nabla u)\nu=0\qquad$
$\displaystyle\text{on}~{}~{}\Gamma,$ (A.4)

supplied with a suitable boundary condition coupling to the bulk [55]. It is
well-known that near the crack tip, the gradients of solutions to (Appendix
A.2) exhibit a persistent $1/\sqrt{r}$ behaviour, where $r$ is the distance
from the crack tip (cf. [56]).

For Mode III (anti-plane) cracks we consider in the numerics, as discussed in
[53], the PDE (Appendix A.2) then reduces to a Poisson equation, which has a
canonical solution, given by

$u^{\rm lin}_{k}(x)=k\sqrt{r}\sin{\frac{\theta}{2}},$

where $(r,\theta)$ representing standard cylindrical polar coordinates centred
at the crack tip. The scalar parameter $k$ corresponds to the (rescaled)
stress intensity factor (SIF) [57].

##  Appendix Appendix B Proof of Theorem 3.3

###### Proof.

From the definitions (3.4) and (3.39), for sufficiently large $R_{\rm a}$, we
have

$\displaystyle\big{|}\tilde{\eta}^{\textrm{ac}}(u_{h})-\eta^{\textrm{ac}}(u_{h})\big{|}$
$\displaystyle=\sum_{T\in\mathcal{T}_{h}}\Big{(}\omega(T)\log(2+|\tilde{{\ell}}(T)|)\cdot\big{|}\mathscr{F}^{\rm
a}_{\tilde{{\ell}}(T)}(I_{\rm a}u_{h})\big{|}-\sum_{\ell\in
T}\log(2+|\ell|)\cdot\big{|}\mathscr{F}^{\rm a}_{\ell}(I_{\rm
a}u_{h})\big{|}\Big{)}$
$\displaystyle\lesssim\sum_{T\in\mathcal{T}_{h}}(2+|\tilde{{\ell}}(T)|)^{-2}\cdot\|\nabla^{2}\widetilde{\mathscr{F}}^{\rm
a}(I_{\rm a}u_{h})\|_{L^{2}(T)}$
$\displaystyle\lesssim\log(R_{\Omega})\cdot\|\nabla^{2}\widetilde{\mathscr{F}}^{\rm
a}(I_{\rm a}u_{h})\|_{L^{2}(\Omega)},$

where the first inequality follows from the standard interpolation error
estimate. ∎

##  Appendix Appendix C Numerical Supplements

To validate the central assumption of Theorem 3.1, i.e., that $u_{h}$ exhibits
the same decay estimates as $u$, we present numerical results for all blended
a/c coupling methods used to model various defect cases in this study. The
setting details for the numerical experiments can be found in Section 5. The
decay results for all defect cases are plotted in the following figures, where
the $x$-axis represents the distance from the ”center” of the defects. In all
cases, we observe that $u_{h}$ exhibits decay rates that are consistent with
those of $u$, as assumed in Theorem 3.1.

(a) BQCE

(b) BQCF

(c) BGFC

Figure A.1: Numerical verification of the decay rate of $u$ and $u_{h}$ for
micro-crack.

(a) BQCE/BGFC

(b) BQCF

Figure A.2: Numerical verification of the decay rate of $u$ and $u_{h}$ for
anti-plane screw dislocation. In this case BGFC is identity to BQCE due to the
artefact of the anti-plane setting.

(a) BQCE

(b) BQCF

(c) BGFC

Figure A.3: Numerical verification of the decay rate of $u$ and $u_{h}$ for
anti-plane crack.

## References

  * [1] E. Van Der 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, Modelling and Simulation in Materials Science and Engineering 28 (4) (2020) 043001.
  * [2] H. Chen, M. Liao, H. Wang, Y. Wang, L. Zhang, Adaptive qm/mm coupling for crystalline defects, Comput. Methods Appl. Mech. Engrg. 354 (2019) 351–368.
  * [3] H. Chen, C. Ortner, Y. Wang, Qm/mm methods for crystalline defects. part 3: Machine-learned mm models, Multiscale Modeling & Simulation 20 (4) (2022) 1490–1518.  doi:10.1137/21M1441122.
  * [4] V. Ehrlacher, C. Ortner, A. V. Shapeev, Analysis of boundary conditions for crystal defect atomistic simulations, Arch. Ration. Mech. Anal. 222 (2016) 1217–1268.
  * [5] 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.
  * [6] 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 (2018) A2087–A2119.
  * [7] M. Luskin, C. Ortner, Atomistic-to-continuum coupling, Acta Numerica 22 (2013) 397–508.
  * [8] D. Olson, P. B. Bochev, M. Luskin, A. V. Shapeev, An optimization-based atomistic-to-continuum coupling method, SIAM J. Numer. Anal. 52 (2014) 2183–2204.
  * [9] I. Tembhekar, J. S. Amelang, L. Munk, D. M. Kochmann, Automatic adaptivity in the fully-nonlocal quasicontinuum method for coarse-grained atomistic simulations, Int. J. Numer. Meth. Engng. 110 (2016) 878–900.
  * [10] C. Ortner, L. Zhang, Atomistic/continuum blending with ghost force correction, SIAM J. Sci. Comput. 38 (2016) A346–A375.
  * [11] 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.
  * [12] M. Luskin, C. Ortner, B. V. Koten, Formulation and optimization of the energy-based blended quasicontinuum method, Comput. Methods Appl. Mech. Engrg. 253 (2013) 160–168.
  * [13] P. Ming, J. Z. Yang, Analysis of a one-dimensional nonlocal quasi-contunuum method, Multiscale Model. Simul. 47 (4) (2009) 1838–1875.
  * [14] C. Ortner, H. Wang, A priori error estimates for energy-based quasicontinuum approximations of a periodic chain, Math. Models Methods Appl. Sc. 21 (12) (2011) 2491–2521.
  * [15] X. H. Li, C. Ortner, A. V. Shapeev, B. V. Koten, Analysis of blended atomistic/continuum hybrid methods, Numer. Math. 134 (2016) 275–326.
  * [16] J. Lu, P. Ming, Convergence of a force-based hybrid method in three dimensions, Comm. Pure Appl. Math. 66 (2012) 83–108.
  * [17] M. Dobson, M. Luskin, Analysis of a force-based quasicontinuum approximation, ESAIM:M2AN 42 (2008) 113–139.
  * [18] 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.
  * [19] C. Ortner, A. V. Shapeev, L. Zhang, (in-)stability and stabilization of qnl-type atomistic-to-continuum coupling methods, Multiscale Model. Simul. 12 (2014) 1258–1293.
  * [20] L. Fang, L. Zhang, Blended ghost force correction method for 3d crystalline defects, COMMUNICATIONS IN COMPUTATIONAL PHYSICS 29 (4) (2021) 1246–1272.  doi:10.4208/cicp.OA-2020-0107.
  * [21] R. E. Miller, E. B. Tadmor, A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods, Modelling Simul. Mater. Sci. Eng. 17 (2009) 1–51.
  * [22] J. T. Oden, S. Prudhomme, Estimation of modeling error in computational mechanics, J. Comput. Phys. 182 (2002) 496–515.
  * [23] S. Prudhomme, P. T. Bauman, J. T. Oden, Error control for molecular statics problems, International Journal for Multiscale Computational Engineering 4 (2006) 647–662.
  * [24] J. T. 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.
  * [25] M. Arndt, M. Luskin, Goal-oriented atomistic-continuum adaptivity for the quasicontinuum approximation, International Journal for Multiscale Computational Engineering 5 (2007) 407–415.
  * [26] C. Ortner, A. V. Shapeev, Analysis of an energy-based atomistic/continuum coupling approximationof a vacancy in the 2d triangular lattice, Math. Comp. 82 (2013) 2191–2236.
  * [27] C. Ortner, A priori and a posteriori analysis of the quasinonlocal quasicontinuum method in 1d, Math. Comp. 275 (2011) 1265–1285.
  * [28] C. Ortner, H. Wang, A posteriori error control for a quasi-continuum approximation of a periodic chain, IMA Journal of Numerical Analysis 34 (2014) 977–1001.
  * [29] 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.
  * [30] Y. Wang, H. Wang, Efficient a posteriori error control of a consistent atomistic/continuum coupling method for two dimensional crystalline defects, ArXiv e-prints arXiv:2211.14143.
  * [31] 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 (2018) 679–709.
  * [32] H. Wang, S. Liu, F. Yang, A posteriori error control for three typical force-based atomistic-to-continuum coupling methods for an atomistic chain, Numer. Math. Theor. Meth. Appl. 12 (2019) 233–264.
  * [33] D. Olson, X. Li, C. Ortner, B. V. Koten, Force-based atomistic/continuum blending for multilattices, Numerische Mathematik 140 (2018) 703–754.
  * [34] D. Olson, C. Ortner, Y. Wang, L. Zhang, Theoretical study of elastic far-field decay from dislocations in multilattices, ArXiv e-prints 1910.12269.
  * [35] W. E, P. Ming, Cauchy-born rule and the stability of crystalline solids: dynamic problems, Acta Math. Appl. Sin. Engl. Ser. 23 (4) (2007) 529–550.  doi:10.1007/s10255-007-0393.
URL http://dx.doi.org/10.1007/s10255-007-0393

  * [36] C. Ortner, F. Theil, Justification of the cauchy–born approximation of elastodynamics, Arch. Rat. Mech. Anal. 207 (2013) 1025–1073.
  * [37] H. Chen, C. Ortner, Qm/mm methods for crystalline defects. part 2: Consistent energy and force-mixing, Multiscale Model. Simul. 15 (2016) 184–214.
  * [38] R. Verfürth, A review of a posteriori error estimation and adaptive mesh-refinement techniques, J. Comput. Appl. Math 50 (1996) 67–83.
  * [39] D. Braess, Finite elements: Theory, fast solvers, and applications in solid mechanics, Cambridge University Press, 2007.
  * [40] E. B. Tadmor, R. E. Miller, Modeling Materials: Continuum, Atomistic and Multiscale Techniques, Cambridge University Press, 2012.
  * [41] K. Fu, M. Liao, Y. Wang, J. Chen, L. Zhang, Adaptive multigrid strategy for geometry optimization of large-scale three dimensional molecular mechanics, ArXiv e-prints arXiv:2105.02672.
  * [42] M. Buze, T. Hudson, C. Ortner, Analysis of an atomistic model for anti-plane fracture, Math. Models Methods Appl. Sc. 29 (13) (2019) 2469–2521.
  * [43] J. Braun, T. Hudson, C. Ortner, Asymptotic expansion of the elastic far-field of a crystalline defect, ArXiv e-prints 2108.04765.
  * [44] C. Ortner, L. Zhang, Energy-based atomistic-to-continuum coupling without ghost forces, Comput. Methods Appl. Mech. Engrg. 279 (2014) 29–45.
  * [45] C. Ortner, The role of the patch test in 2d atomistic-to-continuum coupling methods, ESAIM: Math. Model. Numer. Anal. 46 (2012) 1275–1319.
  * [46] X. H. Li, M. Luskin, C. Ortner, A. V. Shapeev, Theory-based benchmarking of the blended force-based quasicontinuum method, Comput. Methods Appl. Mech. Engrg. 268 (2014) 763–781.
  * [47] V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, M. Ortiz, An adaptive finite element approach to atomic-scale mechanics–the quasicontinuum method, Journal of the Mechanics and Physics of Solids 47 (1999) 611–642.
  * [48] W. Dorfler, A convergent adaptive algorithm for poisson’s equation, SIAM JOURNAL ON NUMERICAL ANALYSIS 33 (3) (1996) 1106–1124.  doi:10.1137/0733054.
  * [49] M. Hodapp, Analysis of a sinclair-type domain decomposition solver for atomistic/continuum coupling, Multiscale Modeling & Simulation 19 (4) (2021) 1499–1537.
  * [50] M. Hodapp, Efficient flexible boundary conditions for long dislocations, arXiv preprint arXiv:2105.08798.
  * [51] D. Olson, A. V. Shapeev, P. B. Bochev, M. Luskin, Analysis of an optimization-based atomistic-to-continuum coupling method for point defects, ESAIM:M2AN 50 (2016) 1–41.
  * [52] J. Sethian, Fast marching methods, SIAM review 41 (2) (1999) 199–235.
  * [53] 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.
  * [54] M. Buze, T. Hudson, C. Ortner, Analysis of cell size effects in atomistic crack propagation, ESAIM: Math. Model. Numer. Anal. 54 (2020) 1821–1847.
  * [55] L. B. Freund, Dynamic fracture mechanics, Cambridge university press, 1998.
  * [56] J. R. Rice, et al., Mathematical analysis in the mechanics of fracture, Fracture: an advanced treatise 2 (1968) 191–311.
  * [57] C.-T. Sun, Z. Jin, Fracture mechanics, Academic Press, 2011.