Theory-based Benchmarking of the Blended Force-Based Quasicontinuum Method
Xingjie Helen Li, Mitchell Luskin, Christoph Ortner, Alexander V. Shapeev
TTHEORY-BASED BENCHMARKING OF THE BLENDED FORCE-BASEDQUASICONTINUUM METHOD
XINGJIE LI, MITCHELL LUSKIN, CHRISTOPH ORTNER, AND ALEXANDER V. SHAPEEV
Abstract.
We formulate an atomistic-to-continuum coupling method based on blending atomisticand continuum forces. Our precise choice of blending mechanism is informed by theoretical predic-tions. We present a range of numerical experiments studying the accuracy of the scheme, focusingin particular on its stability. These experiments confirm and extend the theoretical predictions,and demonstrate a superior accuracy of B-QCF over energy-based blending schemes. Introduction
Atomistic-to-continuum coupling methods (a/c methods) have been proposed to increase thecomputational efficiency of atomistic computations involving the interaction between local crystaldefects with long-range elastic fields [6, 7, 15, 18, 22, 29, 30, 40]; see [26] for a recent review of a/ccoupling methods and their numerical analysis. Energy-based methods in this class, such as thequasicontinuum model (denoted QCE [41]), exhibit spurious interfacial forces (“ghost forces”) evenunder uniform strain [8,39]. The effect of the ghost force on the error in computing the deformationand the lattice stability by the QCE approximation has been analyzed in [8–10, 31], where latticestability refers to the positive definiteness of the Hessian matrix of the total potential energy. Thedevelopment of more accurate energy-based a/c methods is an ongoing process [5,15,20,34,37,38,40].An alternative approach to a/c coupling is the force-based quasicontinuum (QCF) approxima-tion [7, 11, 12, 25, 29], but the non-conservative and indefinite equilibrium equations make the iter-ative solution and the determination of lattice stability more challenging [12–14]. Indeed, it is anopen problem whether the (sharp-interface) QCF method is stable in dimension greater than one.Although some recent results in this direction exist [24], it is still unclear to what extent they canbe extended for general atomistic domains and in the presence of defects.Many blended a/c coupling methods have been proposed in the literature, e.g., [1–4, 16, 23, 35,36, 42]. In [21], we formulated a blended force-based quasicontinuum (B-QCF) method, similar tothe method proposed in [25], which smoothly blends the forces of the atomistic and continuum
Xingjie Li, 182 George St., Providence, RI 02912, USA, xingjie [email protected]. Luskin (Corresponding Author), 127 Vincent Hall, 206 Church St. SE, Minneapolis, MN 55455,USA, [email protected], Phone: 612-625-6565, FAX 612-626-2017C. Ortner, Mathematics Institute, Zeeman Building, University of Warwick, Coventry CV4 7AL,UK, [email protected]. V. Shapeev, 127 Vincent Hall, 206 Church St. SE, Minneapolis, MN 55455, USA,[email protected]
Date : November 9, 2018.2000
Mathematics Subject Classification.
Key words and phrases. quasicontinuum, error analysis, atomistic to continuum, embedded atom model, quasi-nonlocal.This work was supported in part by the NSF PIRE Grant OISE-0967140, DOE Award de-sc0002085, and AFOSRAward FA9550-12-1-0187. CO was supported by EPSRC grant EP/H003096 “Analysis of atomistic-to-continuumcoupling methods.” a r X i v : . [ m a t h . NA ] S e p HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 2 model instead of the sharp transition in the QCF method. Under the simplifying assumptionthat deformation is homogeneous, we established sharp conditions under which a linearized B-QCF operator is positive definite, which effectively guarantees stability of the numerical scheme.Surprisingly, the required blending width to ensure positive definiteness of the linearized B-QCFoperator is asymptotically small (however typical prefactors in the relative size of the blendingregion are not predicted by the theory). The one-dimensional theory developed in [21] is completeand agrees with the numerical experiments. However, the two-dimensional theory was based on aconjecture that has been proved only in a particular case (see Remark 3.1 for more details) andtherefore requires numerical validation.In the present paper, we present focused numerical experiments to validate and extend thetheoretical conclusions in [19, 21]. In particular, we study (i) whether stability of the B-QCFmethod in 2D can be systematically improved with increasing the blending width, (ii) whethera relatively narrow blending, as suggested by the theory, is enough in practice, and (iii) whetherusing the quintic spline (that has the regularity assumed in the theory) has advantages over thecubic spline. In addition we provide accuracy benchmarks similar to those in [27]. Our numericalbenchmarks demonstrate that the B-QCF scheme is a practical a/c coupling mechanism withperformance (accuracy versus computational cost) superior to energy-based blending schemes.1.1.
Summary.
In section 2, we introduce the B-QCF model for a 1D atomistic chain. We statethe asymptotically optimal condition on the blending size in Theorem 2.1 and apply a uniformexpansion to the atomistic chain in subsection 2.2. The critical strain errors between the atomisticand B-QCF models with different blending size are computed in this subsection. The numericalresults perfectly match the analytic prediction, that is, the errors decay polynomially in terms ofthe blending size.In section 3, we establish the B-QCF model for a 2D hexagonal lattice. We state sufficient andnecessary conditions on the blending width under which the B-QCF operator is positive definite.To numerically investigate the positive-definiteness of the B-QCF operators in 2D, we apply threedifferent classes of deformations to the perfect lattice, which are the uniform expansion, two typesof shear deformation, and a general class of homogeneous deformations. The results of 2D uniformexpansion are similar to those of the 1D example, and they agree with the theoretical conclusionswell.The stability regions of the different models under homogeneous deformations are consistentwith the analytic prediction. By using a small blending region, the 2D B-QCF operator becomesalmost as stable as the atomistic model, compared to the fact that the stability region of the force-based quasicontinuum (QCF) method, i.e., the B-QCF method without blending region, is a propersubset of the fully atomistic model [12–14]. However, the stability error under shear deformationfor the B-QCF operator seems to only depend linearly on the system size, which is observed fromthe numerical experiments.In section 4, we implement the B-QCF method from a practical point of view. We briefly reviewthe accuracy results in terms of computational cost, i.e., the total number of degrees of freedom DoF,and then include some numerical experiments for a di-vacancy and microcrack to demonstrate thesuperior accuracy of B-QCF over other a/c coupling schemes that we have investigated previouslyin [27].
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 3 The B-QCF Operator in D. Notation.
We denote the scaled reference lattice by (cid:15) Z := { (cid:15)(cid:96) : (cid:96) ∈ Z } . We apply a macro-scopic strain F > y F := F (cid:15) Z = ( F (cid:15)(cid:96) ) (cid:96) ∈ Z . The space U of 2 N -periodic zero mean displacements u = ( u (cid:96) ) (cid:96) ∈ Z from y F is given by U := (cid:26) u : u (cid:96) +2 N = u (cid:96) for (cid:96) ∈ Z , and (cid:80) N(cid:96) = − N +1 u (cid:96) = 0 (cid:27) , and we thus admit deformations y from the space Y F := { y : y = y F + u for some u ∈ U } . We set (cid:15) = 1 /N throughout so that the reference length of the computational cell remains fixed.We define the discrete differentiation operator, D u , on periodic displacements by( D u ) (cid:96) := u (cid:96) − u (cid:96) − (cid:15) , −∞ < (cid:96) < ∞ . We note that ( D u ) (cid:96) is also 2 N -periodic in (cid:96) and satisfies the zero mean condition. We will oftendenote ( D u ) (cid:96) by Du (cid:96) . We then define (cid:0) D (2) u (cid:1) (cid:96) and (cid:0) D (3) u (cid:1) (cid:96) for −∞ < (cid:96) < ∞ by (cid:16) D (2) u (cid:17) (cid:96) := Du (cid:96) +1 − Du (cid:96) (cid:15) ; (cid:16) D (3) u (cid:17) (cid:96) := D (2) u (cid:96) − D (2) u (cid:96) − (cid:15) . To make the formulas more concise, we sometimes denote Du (cid:96) by u (cid:48) (cid:96) , D (2) u (cid:96) by u (cid:48)(cid:48) (cid:96) , etc., when thereis no confusion in the expressions.For a displacement u ∈ U and its discrete derivatives, we employ the weighted discrete (cid:96) p(cid:15) and (cid:96) ∞ (cid:15) norms by (cid:107) u (cid:107) (cid:96) p(cid:15) := (cid:32) (cid:15) N (cid:88) (cid:96) = − N +1 | u (cid:96) | p (cid:33) /p for 1 ≤ p < ∞ , (cid:107) u (cid:107) (cid:96) ∞ (cid:15) := max − N +1 ≤ (cid:96) ≤ N | u (cid:96) | , and the weighted inner product for (cid:96) (cid:15) is (cid:104) u , w (cid:105) := N (cid:88) (cid:96) = − N +1 (cid:15)u (cid:96) w (cid:96) . The B-QCF Operator.
We consider a one-dimensional (1D) atomistic chain with periodicity2 N , denoted y ∈ Y F , under second-neighbor pair interaction. The total atomistic energy per periodof y is given by E a ( y ) − (cid:15) (cid:80) N(cid:96) = − N +1 f (cid:96) y (cid:96) , where E a ( y ) = (cid:15) N (cid:88) (cid:96) = − N +1 (cid:2) φ ( y (cid:48) (cid:96) ) + φ ( y (cid:48) (cid:96) + y (cid:48) (cid:96) − ) (cid:3) (2.1)for external forces f (cid:96) and a two-body potential φ ∈ C (0 , + ∞ ) such as the Morse potential givenby (2.12). Implicitly we also assume that φ ( r ) , φ (cid:48) ( r ) and φ (cid:48)(cid:48) ( r ) decay rapidly as r increases, so thatwe only have to take into account first and second neighbors.The equilibrium equations are given by the force balance at each atom: F a(cid:96) + f (cid:96) = 0 where F a (cid:96) ( y ) := − (cid:15) ∂ E a ( y ) ∂y (cid:96) = 1 (cid:15) (cid:110) (cid:2) φ (cid:48) ( y (cid:48) (cid:96) +1 ) + φ (cid:48) ( y (cid:48) (cid:96) +2 + y (cid:48) (cid:96) +1 ) (cid:3) − (cid:2) φ (cid:48) ( y (cid:48) (cid:96) ) + φ (cid:48) ( y (cid:48) (cid:96) + y (cid:48) (cid:96) − ) (cid:3) (cid:111) . (2.2) HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 4
The linearized equilibrium equations about y F are( L a u a ) (cid:96) = f (cid:96) , for (cid:96) = − N + 1 , . . . , N, where ( L a v ) for a displacement v ∈ U is given by( L a v ) (cid:96) := φ (cid:48)(cid:48) F ( − v (cid:96) +1 + 2 v (cid:96) − v (cid:96) − ) (cid:15) + φ (cid:48)(cid:48) F ( − v (cid:96) +2 + 2 v (cid:96) − v (cid:96) − ) (cid:15) . Here and throughout we use the notation φ (cid:48)(cid:48) F := φ (cid:48)(cid:48) ( F ) and φ (cid:48)(cid:48) F := φ (cid:48)(cid:48) (2 F ), where φ is the potentialin (2.1). We assume that φ (cid:48)(cid:48) F >
0, which holds for typical pair potentials such as the Lennard-Jonespotential under physically relevant deformations. Appropriate extensions of the stability resultsin this paper can likely be obtained for more general smooth deformations by utilizing the moretechnical formalism developed, for example, in [18, 32, 33].The local QC (or Cauchy-Born) approximation (QCL) uses the Cauchy-Born extrapolation rule[40, 41], that is, approximating y (cid:48) (cid:96) + y (cid:48) (cid:96) − in (2.1) by 2 y (cid:48) (cid:96) in our context. Thus, the QCL energy isgiven by E qcl ( y ) = (cid:15) N (cid:88) (cid:96) = − N +1 (cid:2) φ ( y (cid:48) (cid:96) ) + φ (2 y (cid:48) (cid:96) ) (cid:3) . (2.3)Then the local continuum forces F qcl ( y ) are F qcl (cid:96) ( y ) := − (cid:15) ∂ E qcl ( y ) ∂y (cid:96) = 1 (cid:15) (cid:110) (cid:2) φ (cid:48) ( y (cid:48) (cid:96) +1 ) + 2 φ (cid:48) (2 y (cid:48) (cid:96) +1 ) (cid:3) − (cid:2) φ (cid:48) ( y (cid:48) (cid:96) ) + 2 φ (cid:48) (2 y (cid:48) (cid:96) ) (cid:3) (cid:111) . We can similarly obtain the linearized QCL equilibrium equations about the uniform deformation (cid:16) L qcl u qcl (cid:17) (cid:96) = f (cid:96) for (cid:96) = − N + 1 , . . . , N, where the expression of (cid:0) L qcl v (cid:1) (cid:96) with v ∈ U is (cid:16) L qcl v (cid:17) (cid:96) := (cid:0) φ (cid:48)(cid:48) F + 4 φ (cid:48)(cid:48) F (cid:1) ( − v (cid:96) +1 + 2 v (cid:96) − v (cid:96) − ) (cid:15) . The blended QCF (B-QCF) operator is obtained through smooth blending of the atomistic andlocal QC models. Let β : R → R be a “smooth” and 2-periodic blending function, then we define F bqcf (cid:96) ( y ) := β (cid:96) F a (cid:96) ( y ) + (1 − β (cid:96) ) F qcl (cid:96) ( y ) , where β (cid:96) := β ( (cid:15)(cid:96) ). Linearization about y F yields the linearized B-QCF operator( L bqcf v ) (cid:96) := β (cid:96) ( L a v ) (cid:96) + (1 − β (cid:96) )( L qcl v ) (cid:96) . Next, we define the blending region I of width K : I : = (cid:8) (cid:96) ∈ {− N + 1 , . . . , N } : 0 < β (cid:96) + j < j ∈ { , ± , ± } (cid:9) , and K : = the cardinality of the set I , (2.4)so that D ( j ) β (cid:96) = 0 for all (cid:96) ∈ {− N + 1 , . . . , N } \ I and j ∈ { , , } . Thus K is the size of thecompact support of D ( j ) β (cid:96) . It is obvious that K < N . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 5
Positive-Definiteness of the B-QCF Operator.
We proved in [21] that the blendingfunction β can be chosen as a quintic polynomial such that(i) The j th derivatives of β satisfy (cid:107) D ( j ) β (cid:107) (cid:96) ∞ ≤ C β ( Kε ) − j , for j = 1 , , . (2.5)(ii) This estimate is sharp in sense that, if β (cid:96) attains both the values 0 and 1, then (cid:107) D ( j ) β (cid:107) (cid:96) ∞ ≥ ( Kε ) − j , for j = 1 , , . (2.6)A linearized operator L w with w ∈ { a , c , bqcf } , is said to be positive definite in the H norm or coercive if there exists a constant γ > (cid:104) L w u , u (cid:105) ≥ γ (cid:107) D u (cid:107) (cid:96) ε ∀ u ∈ U . (2.7)We have proved an asymptotically optimal stability condition on the blending region size of the 1DB-QCF operator in [21]. Theorem 2.1.
Let I and K be defined as in (2.4) , and suppose that β is chosen to satisfy theupper bound (2.5) . Then there exists a constant C = C ( C β ) , such that (cid:104) L bqcf u , u (cid:105) ≥ (cid:0) c − C | φ (cid:48)(cid:48) F | (cid:2) K − / N / (cid:3)(cid:1) (cid:107) D u (cid:107) (cid:96) (cid:15) ∀ u ∈ U , (2.8) where c = min( φ (cid:48)(cid:48) F , φ (cid:48)(cid:48) F + 4 φ (cid:48)(cid:48) F ) is the atomistic stability constant.Moreover, if β (cid:96) takes both the values and , then there exist constants C , C > , independentof I , N , φ (cid:48)(cid:48) F and φ (cid:48)(cid:48) F , such that (cid:104) L bqcf u , u (cid:105) ≤ (cid:16) c + (cid:110) C − C (cid:104) K − / N / (cid:105)(cid:111) | φ (cid:48)(cid:48) F | (cid:17) (cid:107) D u (cid:107) (cid:96) (cid:15) for some u ∈ U \ { } . (2.9)From the conclusion of Theorem 2.1, we can immediately get the following necessary and sufficientconditions on the blending width K for the operator L bqcf to be coercive. Corollary 2.1.
Suppose that L a is positive-definite and that the blending function is sufficientlysmooth. If the blending size K satisfies K (cid:29) N / , then the B-QCF operator L bqcf is positive-definite and this estimate is asymptotically optimal. D Uniform Expansion Experiments.
We conduct numerical experiments in order toverify our theoretical findings. More precisely, we compare the decay rates of the error in criticalstrain as computed by B-QCF with the theoretically predicted rates as we increase the blendingwidth K .We use two kinds of blending functions: a cubic splineˆ B ( x ) = x < , − x + 3 x ≤ x ≤ , x > , (2.10)and a quintic spline ¯ B ( x ) = x < , x − x + 10 x ≤ x ≤ , x > . (2.11)We scale ˆ B ( x ) and ¯ B ( x ) and define the blending functions for the atomistic chains asˆ β (cid:96) := ˆ B (cid:18) (cid:96)K (cid:19) and ¯ β (cid:96) := ¯ B (cid:18) (cid:96)K (cid:19) for (cid:96) = − N + 1 , . . . , N. HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 6
Therefore, atoms with indices from − N + 1 to 0 belong to the continuum region, from 1 to K − K to N belong to the atomistic region. We note that ¯ B ( x )has three bounded derivatives and hence it satisfies (2.5), whereas for ˆ B ( x ) the second derivativehas a jump, hence the third derivative does not exist. Therefore, we expect that only ¯ β will yieldthe asymptotically optimal stability estimates for the B-QCF method (see [21]).For our interaction potential, we use the Morse potential φ ( r ) = [1 − exp( − α ( r − , (2.12)and we cut-off the interactions beyond the second nearest neighbor interactions.We apply a uniform expansion to the atomistic chain: y F := F (cid:15) Z with Dirichlet boundarycondition: u − N +1 = u N = 0 . (2.13)We then compute the critical strains of the atomistic and B-QCF models with different blendingsize K and fixed N . The critical strains are defined as γ w := max { F > L w ( y G ) is positive definite for all G ∈ [1 , F ) } , (2.14)where w ∈ { a , c , bqcf } denotes the respective model. Remark . The stability bounds in Theorem 2.1 hold also for displacements u satisfying a ho-mogeneous Dirichlet boundary condition. To establish this, we note (1) that the bounds hold forconstant displacements as well, and (2) that any function satisfying (2.13) can be extended to aperiodic function (possibly with a nonzero mean). Hence, Corollary 2.1 also holds for displacements u with homogeneous Dirichlet boundary conditions (2.13).The computational results are shown in Figure 1. In Figure 1(a) we plot the dependence of theerrors of quintic blending on K for different values of α . We see that the graph of the error for thequintic blending is very close to the lower bound K − / as given by (2.8) in Theorem 2.1. Also, theerror is lower for larger α , which is also in accordance with the theoretical results. Indeed, when α is large, the strength of the next-nearest neighbor interaction, φ (cid:48)(cid:48) F , is small relative to the nearestneighbor interaction φ (cid:48)(cid:48) F , which contributes to a better stability of B-QCF according to (2.8).Figure 1(b) shows the results of comparison of the cubic and the quintic blending. We see thatthe cubic blending produces the error that seems to decay slower, like K − . On the other hand,the quantitative difference between cubic and quintic is not large on the example considered. Toobserve a significantly higher accuracy of the quintic spline, the computational domain size N hasto be much larger. In addition, for larger α , N has to be even larger for the quintic blending tohave advantage over the cubic blending.3. The B-QCF Operator in D. The Triangular Lattice.
For some integer N ∈ N and (cid:15) := 1 /N , we define the scaled 2Dtriangular lattice L to be L := A Z , where A := [ a , a ] := (cid:15) (cid:20) / √ / (cid:21) , where a i , i = 1 , A ( − N/ , N/ and L := L ∩ Ω . We furthermore set a = ( − / (cid:15), √ / (cid:15) ) T , then the set of nearest-neighbor directions is given by N := {± a , ± a , ± a } . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 7 −10 −8 −6 −4 −2 K | γ bq c f − γ a | Abs strain error for 1D uniform expansion: N=40000 α =3, γ a =0.19721 α =4, γ a =0.16405 α =5, γ a =0.13592N K −5/2 *16 (a) Quintic spline blending −8 −7 −6 −5 −4 −3 −2 −1 K | γ bq c f − γ a | Abs strain error for 1D uniform expansion: N=40000, α =3 cubic splinequintic spline10*K −5/2 (b) Quintic v.s. Cubic Figure 1. (a) The absolute critical strain errors for a 1D uniform expansion. Weset N = 40 , , ∆ γ = 1 /N where ∆ γ is the strain increment used for testingstability, and γ a and γ bqcf are the critical strains for the atomistic and B-QCFmodels, respectively. The dashed line corresponds to the theoretical asymptote.(b) The absolute critical strain errors of quintic and cubic blending functions with N = 40 ,
000 and α = 3. The solid line corresponds to the theoretical asymptote.The set of next nearest-neighbor directions is given by N := {± b , ± b , ± b } , where b := a + a , b := a + a , and b = a − a . We use the notation N := N ∪ N to denote the directions of the neighboring bonds in theinteraction range of each atom (see Figure 2).We identify all lattice functions v : L → R with their continuous, piecewise affine interpolantswith respect to the canonical triangulation T of R with nodes L .3.2. The Atomistic, Continuum, and Blending Regions.
Let
Hex ( R ) denote the closedhexagon centered at the origin, with sides aligned with the lattice directions a , a , a , and di-ameter 2 R . For R a < R b < N ∈ N , we define the atomistic, blending, and continuum regions,respectively, asΩ a := Hex ( (cid:15)R a ) , Ω b := Hex ( (cid:15)R b ) \ Ω a , and Ω c := clos (Ω \ (Ω a ∪ Ω b )) . We denote the blending width by K := R b − R a . Moreover, we define the corresponding latticesites L a := L ∩ Ω a , L b := L ∩ Ω b , and L c := L ∩ Ω c . For simplicity, we will again use L as the finite element nodes, that is, every atom is a repatom.For a map v : L → R and bond directions r, s ∈ N , we define the finite difference operators D r v ( x ) := v ( x + r ) − v ( x ) (cid:15) and D r D s v ( x ) := D s v ( x + r ) − D s v ( x ) (cid:15) . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 8 a a a a a a b b b b b b (a) Neighbor set (b) Domain decomposition Figure 2. (a) The 12 neighboring bonds of each atom. (b) The periodic referencecell L := L ∩ Ω , the atomistic region Ω a := Hex ( (cid:15)R a ) , and the blending regionΩ b := Hex ( (cid:15)R b ) \ Ω a . Here, N = 32 , R a = 3, R b = 7, and K = 4.We define the space of all admissible displacements, U , as all discrete functions L → R whichare Ω-periodic and satisfy the mean zero condition on the computational domain: U := (cid:110) u : L → R : u ( x ) is Ω-periodic and (cid:80) x ∈L u ( x ) = 0 (cid:111) . For a given matrix B ∈ R × , det( B ) >
0, we admit deformations y from the space Y B := (cid:8) y : L → R : y ( x ) = Bx + u ( x ) ∀ x ∈ L , for some u ∈ U (cid:9) . For a displacement u ∈ U and its discrete directional derivatives, we employ the weighted discrete (cid:96) (cid:15) and (cid:96) ∞ (cid:15) norms given by (cid:107) u (cid:107) (cid:96) (cid:15) := (cid:32) (cid:15) (cid:88) x ∈L | u ( x ) | (cid:33) / , (cid:107) u (cid:107) (cid:96) ∞ (cid:15) := max x ∈L | u ( x ) | , and (cid:107) D u (cid:107) (cid:96) (cid:15) := (cid:32) (cid:15) (cid:88) x ∈L (cid:88) i =1 | D a i u ( x ) | (cid:33) / . The inner product associated with (cid:96) (cid:15) is (cid:104) u , w (cid:105) := (cid:15) (cid:88) x ∈L u ( x ) · w ( x ) . The B-QCF operator.
The total scaled atomistic energy for a periodic computational cellΩ is E a ( y ) = (cid:15) (cid:88) x ∈L (cid:88) r ∈N φ ( D r y ( x )) = (cid:15) (cid:88) x ∈L (cid:88) i =1 (cid:2) φ ( D a i y ( x )) + φ ( D b i y ( x )) (cid:3) , (3.1)where φ ∈ C ( R ), for the sake of simplicity. Typically, one assumes φ ( r ) = ϕ ( | r | ); the more generalform we use gives rise to a simplified notation; see also [33]. We define φ (cid:48) ( r ) ∈ R and φ (cid:48)(cid:48) ( r ) ∈ R × to be, respectively, the gradient and hessian of φ . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 9
The equilibrium equations are given by the force balance at each atom, F a ( x ; y ) + f ( x ; y ) = 0 , for x ∈ L , (3.2)where f ( x ; y ) are the external forces and F a ( x ; y ) are the atomistic forces (per unit area (cid:15) ) F a ( x ; y ) := − (cid:15) ∂ E a ( y ) ∂y ( x )= − (cid:15) (cid:88) i =1 (cid:104) φ (cid:48) ( D a i y ( x )) + φ (cid:48) ( D − a i y ( x )) (cid:105) − (cid:15) (cid:88) i =1 (cid:104) φ (cid:48) ( D b i y ( x )) + φ (cid:48) ( D − b i y ( x )) (cid:105) . Again, since u = y − y B , where y B ( x ) = Bx , is assumed to be small, we linearize the atomisticequilibrium equation (3.2) about y B :( L a u a ) ( x ) = f ( x ) , for x ∈ L , where ( L a u ) ( x ), for a displacement u , is given by( L a u ) ( x ) = − (cid:88) i =1 φ (cid:48)(cid:48) ( Ba i ) D a i D a i u ( x − a i ) − (cid:88) i =1 φ (cid:48)(cid:48) ( Bb i ) D b i D b i u ( x − b i ) , for x ∈ L . We use the Cauchy-Born extrapolation rule to approximate the nonlocal atomistic model by alocal continuum Cauchy-Born model [29, 39, 41]. Using the bond density lemma [33, Lemma 3.2](see also [37]), we can write the total QCL energy (the discretized Cauchy-Born energy) as a sumof the bond density integrals E c ( y ) = 1Ω (cid:90) Ω (cid:88) r ∈N φ ( ∂ r y ) dx = (cid:88) x ∈L (cid:88) r ∈N (cid:90) φ (cid:0) ∂ r y ( x + tr ) (cid:1) dt, (3.3)where the factor Ω := √ / L and ∂ r y ( x ) := ddt y ( x + tr ) | t =0 denotes the directional derivative. We compute the continuum force F c ( x ; y ) = − (cid:15) ∂ E c ∂y ( x ) , and linearize the force equation about the uniform deformation y B to obtain( L c u c ) ( x ) = f ( x ) , for x ∈ L . To formulate the B-QCF method, we let the blending function β ( s ) : R → [0 ,
1] be a “smooth”,Ω-periodic function. Then, the (nonlinear) B-QCF forces are given through a convex combinationof F a ( x ; y ) and F c ( x ; y ): F bqcf ( x ; y ) := β ( x ) F a ( x ; y ) + (1 − β ( x )) F c ( x ; y ) , and linearizing the equilibrium equation F bqcf + f = 0 about y B yields( L bqcf u bqcf )( x ) = f ( x ) , for x ∈ L , where ( L bqcf u )( x ) = β ( x )( L a u )( x ) + (1 − β ( x ))( L c u )( x ) . (3.4)The 2D blending function in our computational experiments will be defined radially using cubicand quintic splines:ˆ β ( x ) := ˆ B (cid:18) (cid:15)R b − | x | (cid:15)R b − (cid:15)R a (cid:19) and ¯ β ( x ) := ¯ B (cid:18) (cid:15)R b − | x | (cid:15)R b − (cid:15)R a (cid:19) , HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 10 where ˆ B ( x ) is given by (2.10) and ¯ B ( x ) is given by (2.11). The function ¯ β ( x ) has the smoothnessand satisfies 2D versions of the scaling bounds (2.5) needed for Theorem 3.1 below, whereas ˆ β ( x )does not have a bounded third derivative. We therefore can expect that ˆ β ( x ) will give a largererror asymptotically as compared to ¯ β ( x ).3.4. Positivity of the B-QCF operator in 2D.
Necessary and sufficient conditions for L bqcf tobe positive-definite are given in [21]. To make this paper more concise, we only state the conclusionswithout proof. First, we state a lower bound for (cid:104) L bqcf u , u (cid:105) : Theorem 3.1.
Suppose that β ∈ C and satisfies the scaling bounds (2.5) ; then, (cid:104) L bqcf u , u (cid:105) ≥ γ bqcf (cid:107) D u (cid:107) (cid:96) (cid:15) , where γ bqcf := ˜ γ − C (cid:2) K − / R / b | log( R b /N ) | / (cid:3) , (3.5) where C is a generic constant independent of N , and ˜ γ is the coercivity constant for the operator ˜ L : (cid:104) ˜ L u , u (cid:105) := (cid:104) L c u , u (cid:105) − (cid:15) (cid:88) i =1 (cid:88) x ∈L β ( x − a ) (cid:12)(cid:12) D a i D a i +1 u ( x − a − a ) (cid:12)(cid:12) b i ≥ ˜ γ (cid:107) D u (cid:107) (cid:96) (cid:15) ∀ u ∈ U . One can see very clearly that, whenever N is polynomial in R b and K (cid:29) R / b , then L bqcf can beexpected to be coercive. Both are natural and easy to achieve. We can thus deduce the followingresult for the coercivity of L bqcf : Corollary 3.1.
Suppose that ˜ L is positive-definite and that the blending function β ∈ C andsatisfies the scaling bounds (2.5) . Let the number of atoms R a along the radius be of order N α with ≤ α ≤ . If the blending width K satisfies K (cid:29) | log N | / , α = 0 , | log N | / N α/ , < α < ,N / , α = 1 , then the B-QCF operator L bqcf is positive-definite.Remark . (a) The stability result of Theorem 3.1, and hence of Corollary 3.1, is based on the conjecturethat the operator ˜ L is stable. In [21] we show that ˜ L is indeed stable whenever nearest-neighbor interactions dominate.Moreover, based on the analysis and numerical experiments in [33] for a similar linearizedoperator, we expect that the region of stability for ˜ L is the same as for L a as N, R a , R b → ∞ .We therefore expect that the result of Theorem 3.1 holds (up to a controllable error) ifcoercivity of ˜ L is replaced by coercivity of L a in the hypothesis.(b) One can apply the argument of Remark 2.1 to conclude that the results of Theorem 3.1and Corollary 3.1 are valid for homogeneous Dirichlet boundary conditions as well.By constructing a radial counterexample similar to our 1D counterexample, we can observe thatour conditions in Corollary 3.1 are essentially necessary. HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 11
Theorem 3.2.
Suppose that L a is positive-definite and that the blending function β ∈ C andsatisfies the scaling bounds (2.5) . The number of atoms R a along the radius is of order N α with < α ≤ . If the blending width K is K (cid:28) N α/ , then the B-QCF operator L bqcf cannot bepositive-definite and we can construct a radial counterexample in this case. We note that there is a gap between the necessary and sufficient conditions for 0 < α <
1. Inaddition, we have no necessary condition for α = 0, which corresponds to a fixed atomistic coreindependent of the reference cell Ω.3.5.
2D numerical experiments for B-QCF operators.
In this subsection, we will continuethe numerical experiments for the 2D B-QCF models to verify the theoretical findings by comparingthe decay rates of the error in critical strain as computed by B-QCF with the theoretically predictedrates as we increase the blending width K .(1) Uniform expansion.
We first consider the simplest 2D deformation: we apply a uniform expansion y ( x ) = Bx with B = γ (cid:18) (cid:19) to the perfect lattice L with Dirichlet boundary condition: u ( x ) = 0 ∀ x ∈ ∂ Ω . (3.6)Then we compute the critical strains γ of the atomistic and B-QCF models with differentblending region width K .We note that the 2D conclusions also depend on the size of the atomistic region. Thereforewe let R a = K / in order to narrow the dependence only to the blending width K . Thenthe asymptotical term in (3.5) for sufficiently large N is approximately K − / R / b | log( R b /N ) | / = K − / ( R a + K ) / | log( R b /N ) | / ≈ K − / R / a = K − / = R − a , which means that the error in γ bqcf is systematically improvable.The choice of scaling R a = K / is motivated by the results in [17] which indicate that,generically, one should expect an O ( R − a ) error in the regions of stability between theinfinite lattice atomistic model and the atomistic model in a domain with radius R a . Inthe computation, we assign integer values for K and use the rounded values for R a , that is R a = (cid:98) K / (cid:99) .The critical strains are defined as γ w := max (cid:8) ¯ γ > L w ( B x ) is positive definite for γ ∈ [0 , ¯ γ ) (cid:9) , (3.7)where w ∈ { a , bqcf } denote the models. Here we use the MATLAB function eigs [28] tocompute the smallest eigenvalue of the symmetric part of L w ( B x ) and thus determine thepositive-definiteness of L w ( B x ).We also define the increment of the strain γ in each step by ∆ γ . The results in [10,14, 17, 33] suggest that the theoretical increments be of order O ( N − ) (at least, for findingthe critical strain of a uniform lattice), and we set ∆ γ = 10 − which is sufficiently smallconsidering N = 200 or 300 in our experiments.We plot the difference of the critical strains with different blending width K in Figure 3.The numerical critical strain errors in the left figure approach the analytical asymptote as K increases. There are larger fluctuations of errors as compared to the 1D case, which is HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 12 −4 −3 −2 −1 K | γ bq c f − γ a | Abs strain error for 2D uniform expansion: N=500, R a =K α =3, γ a =0.1838 α =4, γ a =0.1589 α =5, γ a =0.1355K −5/2 |R b /N*ln(R b /N)| (a) Quintic spline blending −3 −2 −1 K | γ bq c f − γ a | Abs strain error for 2D uniform expansion: N=500, R a =K , α =3 cubic splinequintic splineK −5/2 |R b /N*ln(R b /N)| *1.3 (b) Quintic v.s. Cubic Figure 3. (a) The absolute critical strain errors for the 2D uniform expansion.We set N = 500, and we denote the critical strains for the atomistic and B-QCFmodels by γ a and γ bqcf , respectively. The dashed line corresponds to the theoreticalasymptote. (b) The absolute critical strain errors for the quintic and cubic blendingfunctions with N = 500 and α = 3. The solid line corresponds to the theoreticalasymptote.likely due to round-off errors in calculating R a = K / . Thus, the slopes of the errors withquintic blending agree with the theoretical prediction in Theorem 3.1. Also, similarly tothe 1D results, the error is smaller when the nearest neighbor interaction dominates (thatis, when α is large).Although the slope of the errors with cubic blending seems to be one half order less thanthat with quintic blending (see Figure 3(b)), the computed errors for cubic blending areslightly smaller for the relatively small N considered. We expect that for a sufficiently largesystem, the quintic blending would be more accurate. In addition, the 2D errors for uniformexpansion are similar to the 1D results. This is reasonable since the 2D uniform expansionis similar to the 1D deformation.(2) Uniform shear deformation.
We now investigate stability of B-QCF under shear deformation. We apply a y-directionalshear deformation to the hexagonal lattice Ω with Dirichlet boundary conditions (3.6). They-directional shear is y ( x ) = ˜ Bx with˜ B = (cid:18) γ (cid:19) . The critical strain errors between the B-QCF and atomistic models with the quintic blendingare plotted in Figure 4.In Figure 4 we plot the critical strain errors in the following three regimes: (1) N increases, R a = const, K = const, (2) all three parameters increase, and (3) N = const, R a and K increases. The choice of constant parameters, K = 2 and R a = 14, does not follow HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 13
40 80 120 160 20010 −2.8 −2.7 −2.6 −2.5 −2.4 −2.3 −2.2 Rel error for y−shear with quintic spline: α =4N | γ bq c f − γ a | / γ a K=2, R a =14K=N , R a =N −1/2 N −1 (a) Error vs N −2.8 −2.7 −2.6 R a | γ bq c f − γ a | / γ a Rel error for y−shear with quintic spline: α =4 N=200, K=R a3/5 −5 N −1 (b) Error vs R a Figure 4.
The relative critical strain error for the y-directional shear deformation. γ a and γ bqcf are the critical strains for the atomistic and B-QCF models respectively.For N = 200, γ a ≈ . . The dashed line corresponds to the theoretical asymptote.The fluctuations in the plotted error for N = const seems to be due to round-offerrors in calculating R a and K . The method parameters were rounded as follows:in (a) K = (cid:98) N / (cid:99) and R a = (cid:98) N / (cid:99) , and in (b) K = (cid:98) R / a (cid:99) .the scaling K ≈ R / a , and was made to show that such nonoptimal parameters do notsignificantly affect the results in this case. The results indicate that the error in this casedepends on N , but does not depend on R a or K . This means that, for shear deformations,the local continuum approximation and its finite element coarse-graining contributes mostof the error.We explain such a qualitative difference between the uniform expansion and the sheardeformation in the following way. For the uniform expansion the onset of instability isdue to competition of interaction of the nearest neighbors (NNs), contributing to stability,and the second nearest neighbors (NNNs), contributing to instability. On the other hand,for the shear deformation the onset of instability is primarily due to competition betweenelongated and compressed NN bonds. Therefore, for the uniform expansion it is importantto reduce the interface error which distorts the NNN interaction, whereas in shear defor-mation the NNN interactions do not contribute significantly to stability errors. Since, forNN interaction, the atomistic, Cauchy-Born and B-QCF models are identical, the stabilityerror only depends on the domain size.(3) Regions of stability.
We now combine the uniform expansion and shear deformation together and study thestability region of L bqcf for a general class of homogeneous deformations. We consider thefollowing family of deformations which involve shear, expansion, and compression. B = (cid:18) s .
10 1 + r (cid:19) . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 14
Applying these specific homogeneous deformations to the hexagonal lattice in the referencecell and again using the Dirichlet boundary condition, we plot the stability regions (regionswhere the operators are positive definite) in Figure ?? . − − − − − − − − − − − − r s N=200,R a =50, (cid:95) =2 Fully ContinuumFully AtomisticB − QCF:K=20QCFContinuum ContinuumQCF QCF Atomistic andB − QCF overlap (a) α = 2 − − − − − − − − − − − − − r s N=200,R a =68, (cid:95) =4 Fully ContinuumFully AtomisticB − QCF: K=10QCFContinuum, Atomistic and B − QCF overlapQCF QCF (b) α = 4 Figure 5.
The stability regions of the different models. These closed curves are theboundaries of the stability regions for the atomistic, B-QCF, and the local continuummodels, respectively. The curves with indicators are for QCF.We observe that the stability regions of the B-QCF model with different blending sizesare all proper subsets of the atomistic model. In addition, the fully atomistic and continuummodels are very close to each other, which agrees with the stability analysis of the perfectlattice [17]. Also, when α increases, which means the next-nearest neighbor interactionsbecome less important, the difference becomes smaller.There is a visible difference in the stability regions between the QCF model and the exactatomistic model, whereas the difference between the B-QCF model and the atomistic modelis almost not seen. This implies that using a blending region can significantly improve thestability properties of the approximation models.(4) Stability of micro-cracks.
The experiments that we have reported up to this point were based on perfect lattices.Now we apply the B-QCF model to lattices with local defects.The atomistic system is as follows. There is a micro-crack in the center of the domain Ωwith length 5, i.e., 5 atoms are removed from the lattice (see Figure 6). Hence, we redefineaccordingly the positions of atoms in the reference configuration x , the interaction energy E a , etc. We impose a vertical stretching B = (cid:20) γ (cid:21) on the lattice and compute thecritical strains γ c > γ c in the following way. Given γ >
0, we use Newton’siteration method to solve the following force equations for y γ with the initial guess y F = B x : F w ( x ; y γ ) = 0 for x ∈ Ω \ ∂ Ω . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 15 −8 −6 −4 −2 0 2 4 6 8−2−10123 Equilibrium of crack−5 with γ =0.001 Figure 6.
The stable equilibrium configuration of the micro-crack with cracklength= 5 and γ = 0 . (cid:96) ∞ (cid:15) norm of the force residual is of order O (10 − ) . We set the tolerance for the (cid:96) ∞ (cid:15) norm of the force residual of the Newton’s iteration tobe 10 − . To prevent the configuration from “jumping out” of the local energy well corre-sponding to the defect under consideration, we require at each step the (cid:96) ∞ (cid:15) norm of forceresidual to be less than 100 and the positive-definiteness of L w ( y ), where y is the currentconfiguration. If any of the two requirements is not met, then the current γ is regarded asan unstable strain. When the force residual is smaller than the tolerance, the configuration y ∗ is thought to be in its equilibrium of the local energy well. Then we check the positivedefiniteness of corresponding operator L w ( y ∗ ) with the equilibrium configuration y ∗ . Thenonlinear critical strain is thus defined as γ c := max { ¯ γ > L w ( y ∗ ) is positive definite for γ ∈ [0 , ¯ γ ) } . The plot of critical strain for the B-QCF models are shown in Figure 7. Even though wechoose the blending width K ≈ R / slightly smaller then our previous choice ( K ≈ R / ),we observe the nonlinear error decays much faster than the theoretical predicted rates andit can reach the strain increment ∆ γ = 10 − . This phenomenon has been observed in [33]and is likely related to superconvergence of local quantities of interest. The indicator of thesuperconvergence is the concentration of the critical eigenmode corresponding to γ c nearthe defect, which is illustrated in Figure 8.We also study the relative errors of the critical strains for two different choices of theblending width, K = 2 and K ≈ R / a + 2. Motivated by the analysis in [33], the size ofthe atomistic core is chosen to be R a = √ N . According to Figure 9, the relative errors for K ≈ R / a + 2 are approximately 10 times smaller than those for K = 2. But both graphsdecay rapidly as N increases. The rate of decay appears to be quadratic. Remark . The numerical computations in this section are conducted without coarsening. Themain reason for this was that coarsening introduces more approximation parameters and potentiallymore fluctuations in the results. Typically, coarsening does not reduce the stability, and we thereforeexpect that our stability results will remain valid for any coarsening. Another reason to discardcoarsening was to better compare the numerics with the theory. However, the purpose of a/ccoupling is to reduce the number of degrees of freedom of an atomistic computation, thereforecoarsening is required when comparing efficiency (i.e., accuracy against the number of degrees offreedom) of different methods.
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 16 −8 −7 −6 −5 −4 −3 K | γ bq c f − γ a | α =4, γ a =0.0370810 −4.5 K −5/2 |R b *log(R b /N)| Figure 7.
The nonlinear critical strain error for vertical stretching. We set N =200, crack length= 5, and R a = max { K , } . γ a , γ bqcf are the critical strains forthe atomistic and B-QCF models, respectively. The dashed line corresponds to thetheoretical asymptote. −15 −10 −5 0 5 10 15 −6 −4−202468 Figure 8.
The zoomed-in critical eigenvector of critical strain of vertically stretch-ing a micro-crack. We set N = 200, crack length= 5, strain increment ∆ γ = 10 − and α = 4 . The Accuracy of B-QCF
In the previous sections, we investigated the positivity of the B-QCF operator. One motivationfor this study is that these experiments fill a gap in our error analysis of the B-QCF method [19]. We
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 17
20 40 80 120 160 20010 −5 −4 −3 −2 N | γ bq c f − γ a | / γ a Rel critical strain error for uniform expansion: R a =N , α =4 K=2K=N +2N −2 /5 Figure 9.
The relative errors of the critical strains of vertically stretching a micro-crack in log scale plot. We set crack length= 5 and α = 4. γ a , γ bqcf are the criticalstrains for the atomistic and B-QCF models, respectively. The dashed line is thetheoretical asymptote.now briefly review these results and then include some numerical experiments demonstrating thesuperior accuracy of B-QCF over other a/c coupling schemes that we have investigated previouslyin [27].4.1. Implementation of the B-QCF method.
Let
V ⊂ L be a set of vacancy sites and L V := L \ V the corresponding lattice with defects. Let B ∈ R × be the applied far-field strain. Weconsider the atomistic problem y a ∈ arg min (cid:8) E a ( y ) : y : L V → R , y ( x ) ∼ Bx as | ξ | → ∞ (cid:9) . (4.1)We remark that one must carefully renormalize E a in order to rigorously make sense of this problem;see e.g. [19] for the details. The vacancy sites are accounted for in the definition of E a by simplyremoving the relevant pair interactions.We wish to approximate this problem with a practical variant (i.e., with coarsening) of the B-QCF method. To that end, we choose R a , R b = R a + K, N ∈ N in such a way that all vacancy sitesare contained in the atomistic region Ω a , which is a hexagon with side length R a . The blendingregion is defined analogously. The full computational domain is given by Ω, which is a hexagonwith side length N . We triangulate Ω in such a way that it matches the canonical triangulation ofthe triangular lattice in Ω.Let T h denote the set of triangles, let N h denote the nodes of the triangulation, and let N free h := N h \ ( V ∪ ∂ Ω) denote the free nodes . HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 18
Let P h denote the space of all functions v h : Ω → R , that are continuous and piecewise affinewith respect to the triangulation T h . The space of admissible trial functions is then given byY h := (cid:8) y h ∈ P h : y h ( x ) = Bx for x ∈ ∂ Ω (cid:9) . Each deformation y h ∈ Y h is understood to be extended by Bx outside of Ω and thereby gives riseto an admissible atomistic configuration.We define the discretized Cauchy-Born energy functional as E c ( y h ) := (cid:88) T ∈T h vol( T ) W (cid:0) ∇ y h | T (cid:1) , where vol( T ) in 2D is the area of the triangle T . We can define the discretized B-QCF operator,for a given blending function β , as follows: F bqcf ( x ; y h ) := (1 − β ( x )) ∂ E a ( y ) ∂y ( x ) (cid:12)(cid:12)(cid:12) y = y h + β ( x ) ∂ E c ( z h ) ∂z h ( x ) (cid:12)(cid:12)(cid:12) z h = y h for x ∈ N free h . In the B-QCF method, we aim to find a solution y bqcf h ∈ Y h satisfying F bqcf ( x ; y bqcf h ) = 0 ∀ x ∈ N free h . (4.2)We remark that this method has essentially five approximation parameters that must be chosencarefully: the atomistic region size R a , the blending width K , the computational domain size N ,the blending function β, and the finite element mesh T h .4.1.1. Practical considerations.
To implement (4.2) in practice, we need to specify further detailsof the method:(1) In our choice of blending function, we deviate from the optimal choice of a C , -blendingfunction and instead choose only a C , blending function, which is more easily constructed.This is justified, firstly, by our foregoing numerical experiments which suggest that littleadditional accuracy in the stability regions can be gained in the pre-asymptotic regime byusing quintic splines (i.e., C , -blending), and secondly, because the consistency error doesnot depend on the regularity of the blending function.We choose the blending function proposed in [27], which minimizes (cid:107)∇ β (cid:107) L , or a discretevariant thereof, in a precomputation step (see [27] for the details).(2) In addition to the blending region Ω b we ensure that two additional “layers” of atoms outsideof it belong to N h . This makes the implementation of the atomistic force contribution in(4.2) straightforward.Moreover, we ensure that the vacancy sites do not affect the forces on atoms x where β ( x ) (cid:54) = 0. This ensures that all the Cauchy-Born force contributions in (4.2) are the correctCauchy-Born forces.(3) To obtain an appropriate initial guess for the B-QCF solutions, we first solve the corre-sponding energy-based blended QCE method (B-QCE) [27] with the same approximationparameters, using a preconditioned line search method. The details are described in [27].The B-QCE solution is then taken as a starting guess for the B-QCF Newton iteration tosolve (4.2). If no B-QCE code is readily available, then a natural alternative would be toimplement a damped Newton method for B-QCF.We remark, that the Jacobian matrix of the B-QCF operator is straightforward to as-semble from the Hessians of the atomistic and Cauchy-Born energy. Nevertheless, for large3D simulations, more sophisticated solution methods may be required. HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 19 (4) We are now only left to choose the remaining approximation parameters R a , K, N and themesh T h .4.2. Error versus computational cost.
We briefly review the main ideas of our analysis in [19]without technical details. A first key result is that if the atomistic solution is stable ( δ E a ( y a ) ispositive definite) and the linearized B-QCF operator δF bqcf ( · ; Bx ) is positive definite, then choosing R a and K sufficiently large implies that δF bqcf ( · ; y a ) is also positive definite, that is, the B-QCFmethod is stable under these conditions. To achieve this in practice, we need to choose K (cid:29) R a (recall from section 4.1.1 that we have chosen a sub-optimal β ).From this stability result, we can deduce the existence of a B-QCF solution in a neighborhoodof the atomistic solution, and an error estimate in terms of the best approximation error (thebest approximation of y a from the finite element space Y h ). and of the modeling error (the forcediscrepancy of the B-QCF and atomistic models). We estimate the error in the strain ∇ y a − ∇ y bqcf h in terms of the “smoothness” of y a , which is measured in terms of bounds on the derivatives ∇ j y a .The derivatives of the discrete functions y a are understood as derivatives of a smooth interpolant.(See [19] for the details.)Dropping an unimportant term for the sake of readability, our error estimate reads (cid:107)∇ y a − ∇ y bqcf h (cid:107) L ( R ) (cid:46) C stab (cid:16) C β (cid:107)∇ y a (cid:107) L ( R \ ω a ) + (cid:107) h ∇ y a (cid:107) L (Ω \ ω a ) + (cid:107)∇ y a (cid:107) L ( R \ ω ) (cid:17) , (4.3)where (cid:107)∇ y a (cid:107) L ( R \ ω a ) measures the modeling error, (cid:107) h ∇ y a (cid:107) L (Ω \ ω a ) the finite element discretiza-tion error and (cid:107)∇ y a (cid:107) L ( R \ ω ) the error in the far-field due to the artificial boundary condition (thetwo latter errors comprise the best approximation error). The domains ω a , ω are slightly smallerhexagonal subsets of, respectively, Ω a and Ω, with comparable side lengths.In addition, C stab is a stability constant that is uniformly bounded for R a (cid:28) K , and C β := K − / R / a log (cid:12)(cid:12) R a /N (cid:12)(cid:12) is a β -dependent prefactor, which arises from a crucial inequality, (cid:107)∇ ( βv ) (cid:107) L ≤ C β (cid:107)∇ v (cid:107) L , in theconsistency analysis of B-QCF.We choose K ≈ R a and N a polynomial of R a (we will see momentarily why this is natural),then C β is uniformly bounded and in addition, we choose R a (cid:28) K , which we require for stability.With this choice, it is easy to see that C β (cid:107)∇ y a (cid:107) L ( R \ ω a ) (cid:46) (cid:107) h ∇ y a (cid:107) L (Ω \ ω a ) (recall that we areworking in units where atomic spacing is 1), and hence we can simply ignore the modeling errorterm from now on.We recall from [27] that the atomistic method (ATM) is given by the B-QCF method with β ≡ . We also recall the corresponding error estimates (dropping less important terms) for the atomistic(ATM) and the B-QCE methods [19, 27] (cid:107)∇ y a − ∇ y atm (cid:107) L ( R ) (cid:46) (cid:107)∇ y a (cid:107) L ( R \ ω ) , (4.4) (cid:107)∇ y a − ∇ y bqce h (cid:107) L ( R ) (cid:46) (cid:107)∇ β (cid:107) L ( R \ ω a ) + (cid:107) h ∇ y a (cid:107) L (Ω \ ω a ) + (cid:107)∇ y a (cid:107) L ( R \ ω ) . (4.5)To better understand the best approximation error, we need to understand the regularity of y a .Since the problems only involve defects with zero Burgers vector, it is reasonable to assume basedon linear elasticity, that |∇ j y a ( x ) | ∼ | x | − j − . (We stress that this estimate only applies in the far-field. In the preasymptotic regime differentrates of decay might be observed, e.g., |∇ j y a ( x ) | ∼ | x | / − j for the micro-crack case discussedin § HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 20
Having this explicit knowledge about the elastic field, we can optimize our choice of finite elementtriangulation. Using the construction in [33] and also used successfully in our B-QCE experimentsin [27], we obtain a triangulation T h (as a function of R b and N ), for which the following estimateholds: (cid:107) h ∇ y a (cid:107) L (Ω \ ω a ) + (cid:107)∇ y a (cid:107) L ( R \ ω ) (cid:46) R − a + N − . Thus, we choose N ≈ R a to balance these two error contributions.Finally, we note that, with this construction, the number of degrees of freedom in Y h , DoF :=dimY h = 2 N free h is approximately equal to DoF ≈ R a . (In particular, the number of degrees offreedom in the atomistic, blending and continuum regions are comparable.)In summary, choosing K ≈ R a , N ≈ R a , the blending function β according to the constructionproposed in [27], and the finite element mesh according to the construction proposed in [33], weobtain from (4.3) the error estimate (cid:107)∇ y a − ∇ y bqcf h (cid:107) L ( R ) (cid:46) DoF − . (4.6)We note that N = R a in the ATM method, and consequently we obtain from (4.5) (cid:107)∇ y a − ∇ y atm (cid:107) L ( R ) (cid:46) DoF − / ; (4.7)thus demonstrating an improved rate of convergence for the B-QCF method in comparison withthe ATM method.We remark that this is optimal for P1-finite element type coarse-graining schemes, as the mod-eling error is in fact dominated by the finite element error. In particular, it is a substantialimprovement over the B-QCE method, for which the corresponding error estimate obtained from(4.4) is (cid:107)∇ y a − ∇ y bqce h (cid:107) L ( R ) (cid:46) DoF − / . We note that the B-QCE method can be shown to have a higher rate of convergence than the ATMmethod for defects with nonzero Burgers vector (such as dislocations) which have a lower rate ofdecay. The finite element coarse-graining of the B-QCE method can more efficiently approximatethe larger region where the strain gradient is significant; [19, 27] for the details.4.3.
Numerical rates.
We test our analytical predictions against the two numerical examples,for which we already tested the B-QCE method in [27]. In both examples, we choose the Morseinteraction potential φ ( r ) = [1 − exp( − α ( r − , with stiffness parameter α = 4.We compare the B-QCF method with a pure atomistic computation on a finite domain, with theQCE and B-QCE methods (cf. [27] for a detailed description of these three methods) and with thepure QCF method, which is simply the B-QCF method with K < β ( x ) ∈ { , } ).Finally, we have also included a highly optimized B-QCE variant where we choose K ≈ R a and N ≈ R a , which is a very unexpected scaling, but yields improved errors in the preasymptoticregime; see [27, Remark 4.3]. We denote this method by B-QCE+ in the error graphs.4.3.1. The di-vacancy example.
We choose the vacancy set V = { , e } and the macroscopic strain B = (cid:18) .
03 0 . . . (cid:19) · B , where B is a minimizer of W (3% uniform stretch and 3% shear from ground state). The setup ofthe B-QCF method for the di-vacancy problem is shown in Figure 10. HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 21 −100 −50 0 50 100−100−50050100 −10 −5 0 5 10−10−50510
Figure 10.
Setup of the B-QCF method for the di-vacancy example, for a specificchoice of approximation parameters, shown in deformed equilibrium. The size/colorof the atoms in the center correspond to decreasing values of (1 − β ( x )). −4 −3 −2 −1 DoF k ∇ y a c − ∇ y a k L Convergence rates for the di−vacancy example
ATMQCEB−QCEB−QCE+QCFB−QCF
DoF −1 DoF −1/2
DoF −1/2
Figure 11.
Plots of computational cost (DoF) versus error in the energy-normfor various a/c coupling methods approximating the di-vacancy problem describedin section 4.3.1.In Figure 11, we plot the degrees of freedom (DoF) against the error in the energy-norm forthe various a/c coupling methods that we consider. As predicted by our analysis, the B-QCFmethod clearly outperforms all other methods, with the exception of the QCF method, which isbarely distinguishable from the B-QCF method in this graph. Unfortunately, we cannot offer asatisfactory theory for the QCF method at present.We also remark that, due to the high consistency error committed in the interface region, theB-QCE does not even outperform a plain atomistic computation in this particular example. (But
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 22 −100 −50 0 50 100−100−50050100
Figure 12.
Setup of the B-QCF method for the micro-crack example, for a specificchoice of approximation parameters, shown in deformed equilibrium. The size/colorof the atoms in the center correspond to decreasing values of (1 − β ( x )).it will clearly outperform the fully atomistic method (ATM) in the micro-crack example, where theelastic field is much more significant.)4.3.2. The micro-crack example.
In the micro-crack (or void) example, we choose the vacancy set V = {− e , . . . , e } and the macroscopic strain B = (cid:18) . . . . (cid:19) · B , where B is a minimizer of W (3% tensile stretch and 3% shear from ground state). The setup ofthe B-QCF method for the micro-crack problem is shown in Figure 12.In Figure 13 we plot the degrees of freedom (DoF) against the error in the energy-norm, for thevarious a/c coupling methods that we consider. In this example the picture is less clear than in thedi-vacancy example due to a more significant preasymptotic regime, which is caused by the moresignificant deformation admitted by the microcrack. In the preasymptotic regime we observe thatthe QCE and B-QCE methods perform much better than expected, but eventually fall back to thepredicted rates. By contrast, the B-QCF and QCF methods display clear systematic convergenceat the predicted rate throughout.We also note that, in this example, the B-QCE+ method performs comparable to the B-QCFand QCF methods, at least in the preasymptotic regime accessible in the experiment.5. Conclusion
We have formulated an atomistic-to-continuum force-based coupling, which we call the blendedforce-based quasicontinuum (B-QCF) method. In this paper, we numerically studied the stabilityas well as accuracy of the B-QCF method. We computed the critical strain errors between theatomistic and B-QCF models with different sizes of the blending region under different types ofdeformations.The main theoretical conclusion in [21] is that the required blending width to ensure coercivityof the linearized B-QCF operator is surprisingly small. For both 1D and 2D uniform expansion,the computational results of the linearized operators perfectly match the analytic predictions. In
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 23 −2 −1 DoF k ∇ y a c − ∇ y a k L Convergence rates for the microcrack example
ATMQCEB−QCEB−QCE+QCFB−QCF
DoF −1 DoF −1/2
DoF −1/2
Figure 13.
Plots of computational cost (DoF) versus error in the energy-normfor various a/c coupling methods approximating the micro-crack problem describedin section 4.3.2.addition, the stability for a general class of homogeneous deformations of the 2D B-QCF operatorbecomes almost the same as that of the atomistic model by using a very small blending region,in contrast to the fact that the stability region of the force-based quasicontinuum (QCF) method,that is, the B-QCF method without blending region, is just a proper subset of the fully atomisticmodel. However, the critical strain error for the B-QCF operator applied to shear deformationseems to only linearly depend on the system size and is thus insensitive to blending width.For the problem of a microcrack in a two-dimensional crystal, we studied the nonlinear stabilityof the B-QCF operators. The critical strain error decays faster than the prediction, and it can beas small as the strain increment. However, we find that the error increases a little bit when theblending size becomes larger, which is possibly due to round-off error.Moreover, we implemented a practical version of the B-QCF method. We briefly reviewed theaccuracy results in terms of computational cost [19]. The numerical experiments, di-vacancy andmicrocrack demonstrate the superior accuracy of B-QCF over other a/c coupling schemes that wehave investigated previously in [27].The BQCF method with a surprisingly small blending region is an appealing choice for numericalsimulations of atomistic multi-scale problems as it is always consistent and can be guaranteed byboth theory and benchmark testing to be positive definite when the fully atomistic operator ispositive definite. 6.
Acknowledgments
We appreciate helpful discussions with Brian Van Koten.
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 24
References [1] S. Badia, P. Bochev, R. Lehoucq, M. L. Parks, J. Fish, M. Nuggehally, and M. Gunzburger. A force-based blend-ing model for atomistic-to-continuum coupling.
International Journal for Multiscale Computational Engineering ,5:387–406, 2007.[2] S. Badia, M. Parks, P. Bochev, M. Gunzburger, and R. Lehoucq. On atomistic-to-continuum coupling by blend-ing.
Multiscale Model. Simul. , 7(1):381–406, 2008.[3] P. T. Bauman, H. B. Dhia, N. Elkhodja, J. T. Oden, and S. Prudhomme. On the application of the Arlequinmethod to the coupling of particle and continuum models.
Comput. Mech. , 42(4):511–530, 2008.[4] T. Belytschko and S. P. Xiao. Coupling methods for continuum model with molecular model.
InternationalJournal for Multiscale Computational Engineering , 1:115–126, 2003.[5] T. Belytschko, S. P. Xiao, G. C. Schatz, and R. S. Ruoff. Atomistic simulations of nanotube fracture.
Phys. RevB , 65, 2002.[6] X. Blanc, C. Le Bris, and F. Legoll. Analysis of a prototypical multiscale method coupling atomistic and con-tinuum mechanics.
M2AN Math. Model. Numer. Anal. , 39(4):797–826, 2005.[7] W. Curtin and R. Miller. Atomistic/continuum coupling in computational materials science.
Modell. Simul.Mater. Sci. Eng. , 11(3):R33–R68, 2003.[8] M. Dobson and M. Luskin. Analysis of a force-based quasicontinuum approximation.
M2AN Math. Model. Numer.Anal. , 42(1):113–139, 2008.[9] M. Dobson and M. Luskin. An analysis of the effect of ghost force oscillation on the quasicontinuum error.
Mathematical Modelling and Numerical Analysis , 43:591–604, 2009.[10] M. Dobson, M. Luskin, and C. Ortner. Accuracy of quasicontinuum approximations near instabilities.
Journalof the Mechanics and Physics of Solids , 58:1741–1757, 2010.[11] M. Dobson, M. Luskin, and C. Ortner. Sharp stability estimates for force-based quasicontinuum methods.
SIAMJ. Multiscale Modeling and Simulation , 8:782–802, 2010.[12] M. Dobson, M. Luskin, and C. Ortner. Stability, instability and error of the force-based quasicontinuum approx-imation.
Archive for Rational Mechanics and Analysis , 197:179–202, 2010.[13] M. Dobson, M. Luskin, and C. Ortner. Iterative methods for the force-based quasicontinuum approximation.
Computer Methods in Applied Mechanics and Engineering , 200:2697–2709, 2011.[14] M. Dobson, C. Ortner, and A. V. Shapeev. The Spectrum of the Force-Based Quasicontinuum Operator for aHomogeneous Periodic Chain.
Multiscale Model. Simul. , 10(3):744–765, 2012.[15] W. E, J. Lu, and J. Yang. Uniform accuracy of the quasicontinuum method.
Phys. Rev. B , 74(21):214115, 2006.[16] J. Fish, M. A. Nuggehally, M. S. Shephard, C. R. Picu, S. Badia, M. L. Parks, and M. Gunzburger. ConcurrentAtC coupling based on a blend of the continuum stress and the atomistic force.
Comput. Methods Appl. Mech.Engrg. , 196(45-48):4548–4560, 2007.[17] T. Hudson and C. Ortner. On the stability of Bravais lattices and their Cauchy–Born approximations.
M2ANMath. Model. Numer. Anal. , 46:81–110, 2012.[18] B. V. Koten and M. Luskin. Analysis of energy-based blended quasicontinuum approximations.
SIAM. J. Numer.Anal. , 49:2182–2209, 2011.[19] H. Li, M. Luskin, C. Ortner, A. V. Shapeev, and B. Van Koten. Blended atomistic/continuum hybrid methods.manuscript.[20] X. H. Li and M. Luskin. A generalized quasi-nonlocal atomistic-to-continuum coupling method with finite rangeinteraction.
IMA Journal of Numerical Analysis , 32:373–393, 2012.[21] X. H. Li, M. Luskin, and C. Ortner. Positive-definiteness of the blended force-based quasicontinuum method.
SIAM J. Multiscale Modeling & Simulation , 10:1023–1045, 2012. arXiv:1112.2528v1.[22] P. Lin. Convergence analysis of a quasi-continuum approximation for a two-dimensional material without defects.
SIAM J. Numer. Anal. , 45(1):313–332 (electronic), 2007.[23] W. K. Liu, H. Park, D. Qian, E. G. Karpov, H. Kadowaki, and G. J. Wagner. Bridging scale methods fornanomechanics and materials.
Comput. Methods Appl. Mech. Engrg. , 195:1407–1421, 2006.[24] J. Lu and P. Ming. Stability of a force-based hybrid method in three dimension with sharp interface.
ArXive-prints , Dec. 2012.[25] J. Lu and P. Ming. Convergence of a force-based hybrid method in three dimensions.
Communications on Pureand Applied Mathematics , 66(1):83–108, 2013.[26] M. Luskin and C. Ortner. Atomistic-to-continuum coupling.
Acta Numerica , 22:397–508, 4 2013.[27] M. Luskin, C. Ortner, and B. Van Koten. Formulation and optimization of the energy-based blended quasicontin-uum method.
Computer Methods in Applied Mechanics and Engineering , 253:160–168, 2013. arXiv: 1112.2377.
HEORY-BASED BENCHMARKING OF THE B-QCF METHOD. 25 [28] MATLAB. version 8 (R2012b) . The MathWorks Inc., Natick, Massachusetts, 2010.[29] R. Miller and E. Tadmor. The quasicontinuum method: overview, applications and current directions.
Journalof Computer-Aided Materials Design , 9:203–239, 2003.[30] R. Miller and E. Tadmor. Benchmarking multiscale methods.
Modelling and Simulation in Materials Scienceand Engineering , 17:053001 (51pp), 2009.[31] P. Ming and J. Z. Yang. Analysis of a one-dimensional nonlocal quasi-continuum method.
Multiscale Model.Simul. , 7(4):1838–1875, 2009.[32] C. Ortner. A priori and a posteriori analysis of the quasinonlocal quasicontinuum method in 1D.
Math. Comp. ,80(275):1265–1285, 2011.[33] C. Ortner and A. V. Shapeev. Analysis of an energy-based atomistic/continuum approximation of a vacancy inthe 2D triangular lattice.
Math. Comp. , 82:2191–2236, 2013. arXiv:1104.0311.[34] C. Ortner and L. Zhang. Construction and sharp consistency estimates for atomistic/continuum coupling methodswith general interfaces: a 2D model problem.
SIAM J. Numer. Anal. , 50, 2012.[35] S. Prudhomme, H. Ben Dhia, P. T. Bauman, N. Elkhodja, and J. T. Oden. Computational analysis of modelingerror for the coupling of particle and continuum models by the Arlequin method.
Comput. Methods Appl. Mech.Engrg. , 197(41-42):3399–3409, 2008.[36] P. Seleson and M. Gunzburger. Bridging methods for atomistic-to-continuum coupling and their implementation.
Communications in Computational Physics , 7:831–876, 2010.[37] A. V. Shapeev. Consistent Energy-Based Atomistic/Continuum Coupling for Two-Body Potentials in One andTwo Dimensions .
SIAM Journal on Multiscale Modeling and Simulation , 9:905–932, 2011.[38] A. V. Shapeev. Consistent energy-based atomistic/continuum coupling for two-body potentials in three dimen-sions.
SIAM Journal on Scientific Computing , 34(3):B335–B360, 2012.[39] V. B. Shenoy, R. Miller, E. B. Tadmor, D. Rodney, R. Phillips, and M. Ortiz. An adaptive finite element approachto atomic-scale mechanics–the quasicontinuum method.
J. Mech. Phys. Solids , 47(3):611–642, 1999.[40] T. Shimokawa, J. Mortensen, J. Schiotz, and K. Jacobsen. Matching conditions in the quasicontinuum method:Removal of the error introduced at the interface between the coarse-grained and fully atomistic region.
Phys.Rev. B , 69(21):214104, 2004.[41] E. B. Tadmor, M. Ortiz, and R. Phillips. Quasicontinuum analysis of defects in solids.
Philosophical MagazineA , 73(6):1529–1563, 1996.[42] S. P. Xiao and T. Belytschko. A bridging domain method for coupling continua with molecular dynamics.