Formulation and optimization of the energy-based blended quasicontinuum method
FFORMULATION AND OPTIMIZATION OF THEENERGY-BASED BLENDED QUASICONTINUUM METHOD
M. LUSKIN, C. ORTNER, AND B. VAN KOTEN
Abstract.
We formulate an energy-based atomistic-to-continuum couplingmethod based on blending the quasicontinuum method for the simulation ofcrystal defects. We utilize theoretical results from [24, 32] to derive optimalchoices of approximation parameters (blending function and finite elementgrid) for microcrack and di-vacancy test problems and confirm our analyticalpredictions in numerical tests. Introduction
A major goal of materials science is to predict the macroscopic properties of ma-terials from their microscopic structure. For this purpose, it is necessary to under-stand the behavior of defects in these materials. We propose a computational tool,the energy-based blended quasicontinuum method (BQCE), for simulating defectssuch as cracks, dislocations, vacancies, and interstitials in crystalline materials.Accurate modeling of the region near a defect requires the use of computationallyexpensive atomistic models. Such models are practical only for small problems.However, a defect may interact with a large region of the material through long-range elastic fields. Thus, accurate simulation of defects requires the use of a largecomputational domain; typically, the size required rules out the use of atomisticmodels for the entire region of interest.Fortunately, the long-range elastic fields generated by a defect are well describedby continuum models which can be efficiently computed using the finite-elementmethod. Thus, defects can be accurately and efficiently simulated by coupled mod-els which use an atomistic model near the defect and a continuum model elsewhere.We call any such model an atomistic-to-continuum coupling .Many atomistic-to-continuum couplings have been proposed in recent years [4–6, 13, 16, 21, 28, 30]; see [19, 31] for a survey of atomistic-to-continuum couplingsand computational benchmark tests. These couplings fall into two major classes:energy-based and force-based. Energy-based couplings provide an approximationto the atomistic energy of a configuration of atoms. Force-based couplings provide anon-conservative force-field which approximates the forces on each atom under theatomistic model. Our BQCE method is an energy-based coupling. Both types of
Date : November 10, 2018.2000
Mathematics Subject Classification.
Key words and phrases. atomistic models, coarse graining, atomistic-to-continuum coupling,quasicontinuum method.ML and BVK were supported in part by the NSF PIRE Grant OISE-0967140, DOE Award de-sc0002085, AFOSR Award FA9550-12-1-0187, and the University of Minnesota SupercomputingInstitute. 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 M. LUSKIN, C. ORTNER, AND B. VAN KOTEN couplings have intrinsic advantages; the development of energy-based couplings isespecially important for finite-temperature applications since equilibrium statisticalproperties and transition rates can be directly approximated [10, 18].Here, and throughout, we only consider concurrent coupling methods. An al-ternative approach are upscaling methods such as [26], which achieve increasedresolution near defects by adding finer scales to a coarse scale description ratherthan by decomposing space into fine (atomistic) and coarse (continuum) scale mod-els.The primary source of error for most (concurrent) energy-based couplings is the ghost force . We say that a coupling suffers from ghost forces if it predicts non-zero forces on the atoms in a perfect lattice. Although many attempts have beenmade to develop an energy-based coupling free from ghost forces such couplings arecurrently known only for a limited range of problems [11, 25, 28, 30].Shapeev’s method [28] applies to one and two-dimensional simple crystals withan atomistic energy based on a pair interaction model, and can be extended to 3Dif a modified “continuum model” is used [29]. GR-AC was proposed in [11,30], andhas recently been implemented for a two dimensional crystal with nearest neighbormany-body interactions in [25]. No ghost force free methods are currently knownfor three-dimensional crystals, for multi-lattice crystals (except in 1D [22]), or foratomistic models with general many-body interactions. The field-based coupling ofIyer and Gavini [16] is another interesting approach; however, it is unclear at presentwhether it is competitive in terms of computational complexity and generality.In our BQCE method, the ghost forces cannot be eliminated but can be con-trolled in terms of an additional approximation parameter (the blending width) [24,32]. BQCE applies to a wide range of problems for which no ghost force free methodsare known; these problems include three-dimensional crystals with general many-body interactions as well as multi-lattices. This makes it an attractive method forsuch challenging and physically important problems.The key feature of BQCE is a blending region where the atomistic and continuumcontributions to the total energy are smoothly mixed. The ghost forces of theBQCE method can be made arbitrarily small by increasing the size of this blendingregion [24,32]. BQCE shares the idea of a blending region with the bridging domainmethod [6], the AtC coupling [4], and the Arlequin method [5]. By contrast, theenergy-based quasicontinuum method (QCE) [21], Shapeev’s method [28], and GR-AC [11, 25, 30] exhibit an abrupt transition between the atomistic and continuummodels. We call any method with a blending region a blended method , and we callthe weights which mix the atomistic and continuum contributions to the energy a blending function .Both the bridging domain method and AtC coupling are very general formula-tions, each of them incorporating BQCE and QCE as special cases. Our BQCE for-mulation provides a set of specific instructions for the successful implementation ofa blended method. We identify two important practical differences between BQCEand the bridging domain and AtC coupling methods. First, the BQCE methodspecifies a strong coupling between the atomistic and continuum regions, whereasweak couplings based on Lagrange multipliers or the penalty method have been usedin most work involving the bridging domain and AtC coupling methods. Second, in The acronym “GR-AC” was introduced in [25]. It stands for “geometry reconstruction-basedatomistic-to-continuum coupling method.”
ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 3
BQCE, we blend the atomistic site-energy with a continuum site-energy based onthe continuum site-energy defined in some formulations of the QCE method (seeSection 3). This guarantees that BQCE correctly predicts the total energy of aperfect lattice subjected to uniform strain.Our approach to blending is supported by rigorous analysis. In [32], we showedthat the ghost force error of BQCE in 1D does indeed decrease with the size ofthe blending region, and we found that the error is minimized when the blendingfunction is a Hermite cubic-spline. We also found that the error of the BQCEmethod in predicting lattice instabilities can be reduced by increasing the size ofthe blending region.In this paper, we have extended the atom-based formulation (4.2) given for1D problems in [32] to a volume-based multiD formulation (4.3) that allows finiteelement coarse-graining in the continuum region.The implementation of BQCE requires the choice of two approximation param-eters: a blending function β and a finite-element mesh T which is used to computethe continuum contribution to the energy. In Section 4.4, we give optimal choicesof β and T to minimize global error norms for the problem of a point defect in a2D crystal, based on theoretical results in [24].We demonstrate the validity of these results in computational test problems inwhich we simulate a microcrack (see Figure 1) and a di-vacancy. Outline.
In Sections 2 and 3, we introduce an atomistic model for simple latticeswith general many-body interactions and its corresponding continuum Cauchy–Born approximation. In Section 4, we give a precise formulation of BQCE in thissetting. In Section 4.4, we offer guiding principles on choosing the blending function β and the mesh T based on our analysis in [24]. The results of this analysis aresummarized in Table 1. In Section 5, we describe the details of our numericalexperiment. We use an atomistic energy based on the Embedded Atom Method.We did not choose our atomistic energy to model any specific physical material;instead, the atomistic energy is a toy model chosen for its simplicity. The observedrates of convergence are in agreement with the rates predicted in Section 4.4. Inparticular, the error of BQCE in the W , semi-norm decreases as DoF − whereDoF is the number of degrees of freedom. We summarize our results and discussopen questions concerning local quantities of interest in the concluding Section 6.2. The Atomistic Energy
Let Λ be a d -dimensional Bravais lattice. We call Λ the reference lattice . In thepresent work, we consider only monatomic crystals . That is, we assume that eachsite of the lattice Λ is the reference position of a single atom, and that all atomsare identical. Let Y := { y : Λ → R d ; y ( ξ ) (cid:54) = y ( η ) for all ξ (cid:54) = η } be the set of deformations of Λ. Let Ω a ⊂ Λ be a finite subset of Λ. We call Ω a the atomistic computational domain .In the Embedded Atom Method (EAM) [8], the energy of Ω a subjected to de-formation y takes the form(2.1) E a ( y ) := (cid:88) ξ ∈ Ω a (cid:26) (cid:88) η ∈ Λ η (cid:54) = ξ φ (cid:0) | y ( η ) − y ( ξ ) | (cid:1) + G (cid:18) (cid:88) η ∈ Λ η (cid:54) = ξ ρ (cid:0) | y ( η ) − y ( ξ ) | (cid:1)(cid:19)(cid:27) , M. LUSKIN, C. ORTNER, AND B. VAN KOTEN where φ is a pair potential, ρ is an electron density function, and G is an embeddingfunction. We call the inner sum(2.2) E a ξ ( y ) := (cid:88) η ∈ Λ η (cid:54) = ξ φ (cid:0) | y ( η ) − y ( ξ ) | (cid:1) + G (cid:18) (cid:88) η ∈ Λ η (cid:54) = ξ ρ (cid:0) | y ( η ) − y ( ξ ) | (cid:1)(cid:19) the atomistic site-energy of atom ξ . We observe that the sum defining the energy ofan atom is finite in practice even though summation ranges over the infinite latticeΛ. This is because the pair potential φ and electron density function ρ are takento have a cut-off radius r c such that φ ( r ) = ρ ( r ) = 0 for all r ≥ r c . For defining the atomistic model and the BQCE method, we can in fact considera more general class of potentials than EAM. We only require that the total energycan be decomposed into a sum of localized site-energies associated with each atom.By localized , we mean that the energy associated with an atom ξ does not dependon the positions of atoms beyond a certain cut-off distance. Such an assumptionmay be violated for certain energies arising from quantum mechanics, but doeshold for most empirical potentials including EAM, the bond-angle potentials, andso forth. In addition, we require that the site-energies are homogeneous ; that is,the energies of atoms which have the same local environment are the same.To make these assumptions precise, we let E a ξ ( y ) denote the site-energy associatedwith atom ξ under deformation y and require that it is of the form(2.3) E a ξ ( y ) := V (cid:0)(cid:8) y ( η ) − y ( ξ ) : η ∈ Λ \ { ξ } , | y ( η ) − y ( ξ ) | < r c (cid:9)(cid:1) , where V is the site-energy potential. We assume that the resulting site-energies aretwice continuously differentiable, that is, E a ξ ∈ C ( Y ). The restricted dependenceof E a ξ on atoms within the cut-off radius may also be expressed in the form ∂ E a ξ ∂y ( η ) ( y ) = 0 for all η with | y ( η ) − y ( ξ ) | ≥ r c . This quantifies the requirement that the site-energy is localized.Given E a ξ we define the energy of Ω a subjected to a deformation y by(2.4) E a ( y ) := (cid:88) ξ ∈ Ω a E a ξ ( y ) . We call an energy of the form (2.4) where E a ξ satisfies (2.3) homogeneous . When wedefine the Cauchy–Born strain energy density corresponding to (2.4) in Section 3,we use that E a is homogeneous: if E a is not homogeneous, the energy per unitvolume in a perfect lattice subjected to uniform strain may be more difficult toobtain (see [1] and references therein for examples). Remark 2.1.
The locality assumption can, in principle, be replaced by an assump-tion that the interaction strength decays sufficiently rapidly with increasing distancebetween atoms.
Remark 2.2.
Homogeneity of the site-energy is our main assumption that is vio-lated for multi-lattices. We show in [24] how to generalize our formulation for thatscenario.
ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 5 The Cauchy–Born Site Energy
BQCE is a coupling of an atomistic energy based on (2.4) with a coarse-grainedcontinuum elastic energy based on the Cauchy–Born strain energy density (3.1).Let vor( ξ ) denote the Voronoi cell of a site ξ ∈ Λ. Then the Cauchy–Born strainenergy density W : R d × d → R ∪ {±∞} corresponding to V is defined by(3.1) W ( F ) := 1 | vor(0) | E a0 ( y F ) , where y F ∈ Y is the homogeneous deformation y F ( ξ ) = F ξ . W ( F ) may be inter-preted as the energy per unit volume in Λ subjected to strain F . Note that theassumption of homogeneity (2.3) ensures that E a ξ ( y F ) = E a0 ( y F ) for all ξ ∈ Λ.We will use the Cauchy–Born strain energy density (3.1) to derive a coarse-grained continuum energy suitable for coupling with the atomistic energy (2.4).First, we define a space of coarse-grained deformations. Let Λ rep ⊂ Λ denote aset of representative atoms (repatoms) , let T be a triangulation with vertices Λ rep ,and let P ( T ) denote the set of all functions y : R d → R d that are continousand piecewise affine with respect to T . We call P ( T ) the set of coarse-graineddeformations . The Cauchy–Born energy of a deformation y ∈ P ( T ) in a domainΩ is then given by E c ( y ) := (cid:90) Ω W ( ∇ y ( x )) dx. The Cauchy–Born approximation is analyzed, for example, in [7, 12, 15].The definition of the QCE method [21] and our construction of the BQCEmethod in the next section use a Cauchy–Born site-energy E c ξ , which is analogousto the atomistic site-energy E a ξ . For y ∈ P ( T ) and ξ ∈ Λ, we define E c ξ by E c ξ ( y ) : = (cid:90) vor ( ξ ) W ( ∇ y ( x )) dx = (cid:88) T ∈T | vor( ξ ) ∩ T | W ( ∇ y | T ) . (3.2)In formula (3.2), | vor( ξ ) ∩ T | denotes the volume of the intersection of vor( ξ ) withthe element T . We observe that the sum on the right hand side of equation (3.2)is finite because only finitely many elements T ∈ T can intersect vor( ξ ).4. The Blended Quasicontinuum Energy
Formulation of the BQCE method.
The BQCE method is an atomistic-to-continuum coupling based on the QCE method of Tadmor et al. [21]. In QCE,the reference domain Ω a is partitioned into an atomistic region A and a continuumregion C , and the QC energy E qc : P ( T ) → R is defined by(4.1) E qc ( y ) := (cid:88) ξ ∈A E a ξ ( y ) + (cid:88) ξ ∈C E c ξ ( y ) . In BQCE, the atomistic and continuum energies per atom are weighted averages.Given a blending function β : Ω a → [0 ,
1] the BQCE energy E β : P ( T ) → R isdefined by(4.2) E β ( y ) := (cid:88) ξ ∈ Ω a β ( ξ ) E c ξ ( y ) + (1 − β ( ξ )) E a ξ ( y ) . M. LUSKIN, C. ORTNER, AND B. VAN KOTEN
We observe that the QCE energy with continuum region C is the same as the BQCEenergy with β chosen as the characteristic function of C . Our formulation of BQCEis similar in spirit to the bridging domain method [6], the AtC coupling [4], andthe Arlequin method [5]; it differs from [4–6] in that we specify a strong couplingbetween the atomistic and continuum models, and accordingly our method is basedon the minimization of an energy. Moreover, our formulation in terms of atomisticand continuum energies per atom guarantees that the BQCE energy is exact underhomogenous deformations.The BQCE energy can be rewritten in the form E β ( y ) = (cid:88) ξ ∈ Ω a β ( ξ ) (cid:90) vor ( ξ ) W ( ∇ y ) dx + (1 − β ( ξ )) E a ξ ( y )= (cid:88) ξ ∈ Ω a (cid:88) T ∈T β ( ξ ) | vor( ξ ) ∩ T | W ( ∇ y | T ) + (cid:88) ξ ∈ Ω a (1 − β ( ξ )) E a ξ ( y )= (cid:88) T ∈T v βT W ( ∇ y | T ) + (cid:88) ξ ∈ Ω a (1 − β ( ξ )) E a ξ ( y ) , (4.3)where the BQCE-effective volume of the element T is defined by(4.4) v βT := (cid:88) ξ ∈ Ω a β ( ξ ) | vor( ξ ) ∩ T | . Remark 4.1.
The triangulation T need not cover the entire domain Ω a . For T atriangulation which covers only part of R d , define Ω( T ) : = ∪ T ∈T T, and P ( T ) : = { y : Ω( T ) ∪ Λ → R d : y piecewise affine w.r.t. T on Ω( T ) } . We observe that E β ( y ) is defined for y ∈ P ( T ) if for every ξ ∈ Ω a such that β ( ξ ) > we have vor( ξ ) ⊂ Ω( T ) . In particular, it is not necessary to assume thatthe triangulation T is refined to atomistic scale anywhere in the domain Ω a . It ispossible that the use of a mesh which is not refined to atomistic scale may make theimplementation of BQCE easier and more efficient in some cases. Remark 4.2 (Multi-lattices) . If one interprets a multi-lattice as a simple latticewith a basis (see, e.g., [9,24,33]), then the formulation of the atomistic energy (2.4) and of the BQCE energy (4.2) require no changes. The only difference is that inthe general multi-lattice case, the site energy E a ξ ( y ) has fewer symmetries (however,see [33] for examples with high symmetry). Far-field boundary conditions.
A typical application of the BQCE methodis the simulation of a defect or defect region in an infinite crystal. To that end,we require far-field boundary conditions at the domain boundary. We propose twochoices.4.2.1.
Dirichlet boundary conditions.
Let
S ⊂ T be a finite subset of T , and letΩ( S ) := ∪ S ∈S S be a polygonal domain. When Dirichlet boundary conditions areimposed, the deformation of the boundary ∂ Ω( S ) of Ω( S ) is fixed to agree withsome y ∈ P ( T ). Precisely, we let(4.5) Adm := (cid:8) y ∈ P ( T ) : y ( x ) = y ( x ) for all x ∈ ∂ Ω( S ) (cid:9) ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 7 denote the space of admissible deformations. We then solve the problem(4.6) Find y ∈ argmin z ∈ Adm E β ( z ) , where we interpret argmin z ∈ Adm E β ( z ) as the set of local minimizers of E β .4.2.2. Periodic boundary conditions.
A popular method to construct artificial far-field boundary conditions is to formulate the problem in a periodic cell. To that end,suppose that Ω( T ) = R d , T is periodic, and let S ⊂ T be the finite element meshon one periodic cell. That is, suppose that for some nonsingular matrix A ∈ R d × d whose columns are elements of Λ, we have T = (cid:91) n ∈ Z d (cid:8) An + S (cid:9) , and that the union is disjoint.Let F ∈ R d × d be nonsingular, and call F a macroscopic strain . Given F , wedefine the admissible set as(4.7) Adm := (cid:8) y ∈ P ( T ) : y ( x + An ) = y ( x ) + F An for all x ∈ R d (cid:9) . The associated variational problem can again be stated as (4.6).The columns of the matrix A should be understood as the directions in which thedisplacement u ( x ) = y ( x ) − F x is periodic. Equivalently, the set (cid:8) Ax : x ∈ [0 , d (cid:9) is a periodic cell. The macroscopic strain F should be understood as the averagestrain imposed on the crystal. If the size of the periodic cell is large, periodicboundary conditions with macroscopic strain F approximate far-field boundaryconditions with uniform strain F imposed at infinity.Note that, while we only discuss a static formulation, we see no obstacle thatprevents both our formulation and analysis to apply also to quasi-static problems, inwhich the energy is minimized at each loading step. However, we stress that furthercareful analysis is required to understand the accuracy, in particular stability, ofthe BQCE method as the deformation approaches bifurcation points (e.g., as in theformation or motion of crystal defects or the propagation of a crack).4.3. Degrees of freedom (DoF).
In Section 4.4, we present error estimates forthe BQCE method in terms of the computational complexity , by which we meanthe computational cost to compute the BQCE energy and its gradient. This isproportional to the number of degrees of freedom of the space of admissible functionsAdm. Similarly, in Section 5, we plot errors against the number of degrees offreedom in the computed solution.The number of degrees of freedom (DoF) is simply the size of the nonlinearsystem that we solve to minimize the BQCE energy. In terms of the notationintroduced in the previous section, this is the number of finite element nodes in thetriangulation S plus the number of lattice sites in Λ \ Ω( S ) for which β ( ξ ) < Computational complexity and optimal parameters.
In [24], we con-jecture an error estimate for BQCE in 2D and 3D. The conjecture is based on ourerror analysis of BQCE in 1D [32] and on the consistency estimates for BQCE in2D and 3D in [24]. Following [23, Sec. 7.1], we use our conjectured error estimateto derive complexity estimates , which are bounds on the error of the BQCE method
M. LUSKIN, C. ORTNER, AND B. VAN KOTEN in terms of the number of degrees of freedom. We use our complexity estimates toguide the choice of optimal approximation parameters T and β .Let y a be a local minimizer of the atomistic energy, and let y β be a local minimizerof BQCE with the same boundary conditions as y a . Let h be the mesh size functionof the triangulation T ( h ( x ) = diam( T ) for a.e. x ∈ T ), and let β be the blendingfunction. Then we conjecture an estimate of the form:Err := (cid:107)∇ y a − ∇ y β (cid:107) L p (cid:46) (cid:107) h ∇ y a (cid:107) L p ( C ) + (cid:107)∇ β (cid:107) L p =: CG + GF , (4.8)for given p ∈ [1 , ∞ ]. For p = 2, we have proven this result in [24]; for p (cid:54) = 2, itwould be difficult to establish this result rigorously; and for p ∈ { , ∞} , it is unclearin what generality the result holds. In (4.8), ∇ y a and ∇ β should be interpretedas the second derivatives of smooth interpolants of y a and β , and C := supp( β ).The first term, CG = (cid:107) h ∇ y a (cid:107) L p ( C ) , is the finite element coarsening error , whilethe second term, GF = (cid:107)∇ β (cid:107) L p , measures the effect of the ghost forces.We now describe the problem of a point defect in a 2D crystal. To quantify thenotion of a point defect, we assume that for some α > |∇ y a ( x ) | (cid:39) r − α , where r = | x | . It has been observed in numerical experiments that α = 2 for a dislocation, andthat α = 3 for a vacancy [14, 23].We assume that the reference domain Ω is a roughly circular region of radius N atomic spacings centered at the origin. We let K > K > β of the form(4.10) β ( x ) := r < K ,β (cid:16) | x |− K K (cid:17) if K ≤ r < K + K , K + K ≤ r, where β : [0 , → [0 ,
1] is a twice continuously differentiable function with β (0) = β (cid:48) (0) = β (cid:48) (1) = 0 and β (1) = 1.We summarize our complexity estimates for BQCE in Table 1. These estimatesare proved in [24]. We distinguish three cases based on the value of γ := αpp +2 . In thesecond column, we give the optimal rates of convergence for BQCE. The optimalrates are attained when K is given in terms of K by the formula appearing in thethird column and when the mesh size function h is given by(4.11) h ( x ) = (cid:18) | x | K (cid:19) γ . Remark 4.3.
When γ > and α > , another choice of parameters results inconvergence at the optimal rate but with a better constant: we use K = K µ for µ = α − p − p in our numerical experiments below. This choice guarantees that the finite elementerror decreases at the same rate as the ghost force error. See [24] for additionalsuggestions related to choosing parameters. ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 9
Table 1.
Complexity estimates and optimal approximation pa-rameters for BQCE [24]. The estimates above are for the problemof a point defect in a two-dimensional crystal. The variables K , K , N , and α are defined in Section 4.4, and γ := αpp +2 . In all cases,the optimal rate of convergence is attained when h = (cid:16) | x | K (cid:17) γ andwhen K is given in terms of K by the formula in the third columnabove.Case Complexity Estimate Optimal Parameters γ > (cid:107)∇ y a − ∇ y β (cid:107) L p (cid:46) DoF max { p − , p − α } K = K γ = 1 (cid:107)∇ y a − ∇ y β (cid:107) L p (cid:46) DoF max { p − , − } K = K ln (cid:16) NK (cid:17) γ < (cid:107)∇ y a − ∇ y β (cid:107) L p (cid:46) DoF max { p − , − } K = K γ N − γ Remark 4.4.
To make the explicit calculation of optimal parameters possible, wehave made simplifying assumptions on the forms of β and h . In (4.10) , we assumethat β is radial. We explain how to construct β given general atomistic, continuum,and blending regions in Section 5.2 below. Similarly, we assume that h is continuouseven though a nonconstant continuous function cannot be the mesh size functionof a triangulation. In practice, we implement BQCE using a triangulation whosemesh size function approximates the continuous h defined in (4.11) ; we explain howto do this in Section 5.2. We expect that the complexity estimates in Table 1 willhold for practical choices of β and h , not only for the simplified β and h used inthe derivation of the estimates. Remark 4.5.
The rates of convergence depend on the geometry of the problemand on the norm in which the error is measured. We observe that the error ofBQCE does not decrease with
DoF when measured in the W , -seminorm; see thecase p = 1 in Table 1. This is because the W − , -norm of the ghost force does notdecrease as the size of the blending region increases. We refer the reader to [24] forfurther discussion of the convergence of BQCE for other geometries and in othernorms. Remark 4.6.
We have proved that linear blending functions have a reduced rateof convergence compared to smooth blending functions of the form (4.10) (see [32]for 1D and [24] for 2D). In Table 2, we give comparative convergence rates forthe problem of a point defect in a two-dimensional crystal with decay α = 3 . Weemphasize, however, that this is only an asymptotic estimate. Indeed in our micro-crack experiment in Section 5.3, we observe a significant preasymptotic regime wherethe two blending functions yield comparable errors.
Remark 4.7.
Although our subsequent numerical experiments are restricted to de-fects with zero Burgers vector and consequently γ > , we consider more generalsituations in the above analysis, in order to emphasize the generality of our ap-proach. The implementation of a/c couplings for defects with non-zero Burgersvector poses additional challenges (what should one use for the reference domain inthe neighborhood of a single dislocation?) that we will address in future work. As linear BQCE smooth BQCE (cid:107)∇ y a − ∇ y β (cid:107) L p (cid:107)∇ y a − ∇ y β (cid:107) L p p = 2 (cid:46) DoF − (cid:46) DoF − p = ∞ (cid:46) DoF − (cid:46) DoF − Table 2.
Complexity estimates comparing the optimized linearand smooth BQCE method for the problem of a point defect in atwo-dimensional crystal with decay α = 3 . a matter of fact, we are not aware of numerical results from a/c coupling methods,simulating a single dislocation core in an otherwise defect-free lattice (as opposedto the simulation of a dislocation pair with total zero Burgers vector). Comparison of BQCE with other methods.
We use the estimates aboveto compare the complexity of BQCE with direct atomistic simulation and withpatch test consistent couplings, such as Shapeev’s method.4.5.1.
Shapeev’s method and other patch test consistent methods.
When α is small(slow decay of the elastic field), BQCE converges at the same rate as patch testconsistent methods such as Shapeev’s method. In fact, if α ≤ W , -error of Shapeev’s method decreases as DoF − with the number ofdegrees of freedom [23]. We predict the same rate of convergence for BQCE when α ≤
2. On the other hand, when α is larger, patch test consistent couplings mayconverge faster than BQCE. When α > W , -error of Shapeev’s method decreases as DoF − α , whereas the W , -error of BQCE decreases as DoF − . Roughly speaking, this is due to thefact that, for small α , the coarse-graining error dominates, whereas for large α , theghost force error dominates.According to the estimates above, Shapeev’s method always converges at least asfast as BQCE. However, we remind the reader that patch test consistent methodsare currently known only for a limited range of problems [11, 25, 28, 30]. As faras we are aware, BQCE is the only convergent energy-based method for 2D or3D crystals with general many-body interactions. Moreover, even when patch testconsistent methods are known, BQCE may be an attractive method for slowlydecaying defects, since it converges at the same rate as patch test consistent methodsfor these problems. See [24] for a more detailed analytic comparison of BQCE withShapeev’s method.4.5.2. Reduced atomistic computations (ATM).
In our numerical experiments be-low, we compare the performance of BQCE with the reduced atomistic computation(ATM) , in which far-field boundary conditions are simulated by imposing Dirichletboundary conditions directly on the atomistic model. Precisely, let Ω ⊂ Λ be an atomistic computational domain , and let F be a macroscopic strain as in Section 4.2.We then define the admissible spaceAdm := (cid:8) y : Λ → R d : y ( ξ ) = F ξ for all ξ ∈ Λ \ Ω } , ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 11
BQCE ATM 2 < α
ATM rate vs. (cid:107)∇ y a − ∇ y β (cid:107) L p (cid:107)∇ y a − ∇ y Ωa (cid:107) L p BQCE rate α > p = 2 (cid:46) DoF − (cid:46) DoF − + (3 − α ) α = 3 same α < α > p = ∞ (cid:46) DoF − (cid:46) DoF − (3 − α ) α = 3 same α < Table 3.
Complexity estimates for the optimized BQCE andATM for the problem of a point defect in a two-dimensional crystalwith decay α > α is defined in Section 4.4.and we solve the variational problem(4.12) Find y ∈ argmin z ∈ Adm E a ( z ) . We regard ATM as the simplest means of simulating far-field boundary condi-tions, and we feel that an atomistic-to-continuum coupling is successful only whenit is more accurate than ATM. In [24], we show that (cid:13)(cid:13) ∇ y Ωa − ∇ y a (cid:13)(cid:13) L p (cid:46) N p +1 − α (cid:39) DoF p + − α , where N is the diameter of Ω, y Ωa solves (4.12), p ∈ [1 , ∞ ], and α >
2. Comparingthe estimate above with Table 1, we observe that BQCE converges more quicklythan ATM when α <
3, at the same rate when α = 3, and more slowly when α > , which we summarize in Table 3. Our numerical experiments concern thecase α = 3, and we do observe that BQCE and ATM converge at the same rate.However, for optimally chosen parameters, BQCE is still much more accurate thanATM by a significant constant factor.5. Numerical Example
Setup of the atomistic model.
In the 2D triangular lattice defined byΛ := (cid:18) / √ / (cid:19) · Z we choose a hexagonal domain Ω a with embedded microcrack as described in Fig-ure 1. The sidelength of the domain is N = 500 atomic spacings, and a de-fect is introduced by removing a segment of atoms at the center of Ω a , Λ def11 := {− e , . . . , e } in the microcrack example in Section 5.3, and Λ def2 := { , e } inthe di-vacancy example in Section 5.4, where the subscripts refer to the number ofatoms removed in the simulation.We supply the domain with Dirichlet boundary conditions by setting y a ( ξ ) = F ξ for ξ ∈ Λ \ Ω a , where F is a macroscopic strain that we will specify below. − −
50 0 50 100 − − − − Figure 1.
Deformed configuration in atomic units of a BQCEsolution for a microcrack in a computational domain with approx-imately 3 N atoms ( N = 100 in this figure, but N = 500 in thebenchmarks described in Sections 5.3 and 5.4). The color and sizeof the atom positions indicate the value of the blending function.The site energy is given by the EAM toy-model (2.2), with φ ( r ) = (cid:2) e − a ( r − − e − a ( r − (cid:3) , ρ ( r ) = e − br , and G (¯ ρ ) = c (cid:2) (¯ ρ − ¯ ρ ) + (¯ ρ − ¯ ρ ) (cid:3) , where a, b, c, ¯ ρ ∈ R are parameters of the model. In all our computational experi-ment we choose a = 4 . , b = 3 , c = 5 , ¯ ρ = 6 e − b , r cut1 = 1 . , r cut2 = 2 . . The model mimics qualitative properties of practical EAM models such as thosedescribed in [17, 34].To simplify the implementation of the model, we have introduced a cut-off in thereference configuration instead of the deformed configuration. In (2.2), we computeonly sums over first, second and third neighbors. Since there are no rearrangementsof atom connectivity in our experiments, we do not expect that this affects theoutcome of the results.5.2.
Setup of the BQCE method.
Our implementation of the BQCE method isbased on (4.3), using standard finite element assembly techniques. The constructionof the blending function is governed by two approximation parameters: • K ∈ N denotes the number of atomic layers surrounding the microcrackwhere β = 1; • K ∈ N denotes the number of atomic layers in the blending region.For the microcrack problem we expect α = 3 in the context of Section 4.4. Hence,according to Table 1, the choice K := K (so that the number of atoms in theatomistic region and blending region are comparable) results in the optimal rate ofconvergence for both p = 2 and p = ∞ . ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 13
Let d hop ( ξ, η ) denote the hopping distance in the triangular lattice (with naturalextension to sets), then we defineΛ a := (cid:8) ξ ∈ Ω a : d hop ( ξ, Λ crack ) ≤ K (cid:9) , andΛ b := (cid:8) ξ ∈ Ω a : K < d hop ( ξ, Λ crack ) ≤ K + K (cid:9) . We consider three choices of the blending function, which are all easily defined forgeneral interface geometries: • QCE: choosing β to be the characteristic function of the atomistic regionΛ a , and K = 0, yields the QCE method defined in (4.1). • Linear Blending:
Let d ( ξ ) denote the hopping distance from the atomisticregion, then we choose β lin ( ξ ) := max(1 , d ( ξ ) /K ) . • Smooth Blending:
The three nearest-neighbour lattice vectors are given by a := (cid:18) (cid:19) , a := (cid:18) / √ / (cid:19) , and a := (cid:18) − / √ / (cid:19) . Define ∆ i β ( ξ ) := β ( ξ + a i ) − β ( ξ ) + β ( ξ − a i ), andΦ( β ) := (cid:88) ξ ∈ Ω a (cid:88) i =1 | ∆ i β ( ξ ) | . Then we define β smooth := argmin (cid:8) Φ( β ) : β ( ξ ) = 0 in Λ a and β ( ξ ) = 1 in Ω a \ Λ a ∪ Λ b (cid:9) . Roughly speaking, this choice generalizes the cubic spline blending sug-gested in [32]. It provides a practical variant for interface geometries thatare not circular.The third approximation parameter is the finite element mesh in the continuumregion. We coarsen the finite element mesh away from the boundary of the blendingregion according to the rule suggested by the complexity estimates in Table 1. Asa matter of fact it turns out that the mesh size growth is too rapid to create shape-regular meshes, hence we also impose the restriction that neighboring element canat most grow by a prescribed factor; this introduces an additional logarithmic factorin the complexity estimates [23, Sec. 7.1].The resulting energy functional is minimized using a preconditioned Pol´ak–Ribi`ere conjugate gradient algorithm described in [31]. We removed the terminationcriterion in this algorithm and allowed it to converge to its maximal precision, thatis, until the numerically computed descent direction ceases to be an actual descentdirection for the energy; this occurs at an accuracy of around 10 − in atomic units.We then perform three Newton steps, after which the residual is in the range ofmachine precision, to ensure that the results are not affected by the accuracy of thenonlinear solver.5.3. Microcrack.
In the microcrack experiment, we remove a long segment ofatoms, Λ def11 = {− e , . . . , e } from the computational domain. The body is thenloaded in mixed mode I & II, by setting F := (cid:18) γ II γ I (cid:19) · F , − − DoF ! ∇ y a c − ∇ y a ! L ATMQCEBQCE: smoothBQCE: linear
DoF − DoF
DoF − DoF − Figure 2.
Convergence rates in the energy-norm (the H -seminorm) for the microcrack benchmark problem described inSection 5.3. This corresponds to the case p = 2 and α = 3 inTable 1. In the above, y a denotes the minimizer of the atomisticenergy, and y ac denotes the corresponding minimizer of the ATM,QCE or BQCE problems.where F ∝ I minimizes the Cauchy–Born stored energy function W , γ I = γ II = 0 . K , and compute the error relativeto the exact atomistic solution.The relative errors in the W , -seminorm are displayed in Figure 2; the relativeerrors in the W , ∞ -seminorm are displayed in Figure 3; the errors in the energy aredisplayed in Figure 4. We observe a significant pre-asymptotic regime where therate of convergence is faster than predicted in our theory. For larger computations,the rate of convergence approaches closely the predicted rate.In particular, it is worth noting that the advantage of a smooth blending functiononly becomes significant in the asymptotic regime and is less pronounced than ourtheory might suggest. Since the precomputation of β smooth is computationallycheap and straightforward, we nevertheless propose the choice β smooth for mostsimulations.Finally, we remark that the initial dip in the energy error in Figure 4 is mostlikely due to a change of sign in the energy error.5.4. Di-vacancy.
We perform a second numerical experiment, in order to demon-strate that the effects we observed in the microcrack example, which are not coveredby our theory, are indeed caused by a pre-asymptotic regime. To this end we con-sider a di-vacancy defect, where only two neighboring sites Λ def2 := { , e } are ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 15 − − − − DoF ! ∇ y a c − ∇ y a ! L ∞ ATMQCEBQCE: smoothBQCE: linear
DoF − DoF − DoF Figure 3.
Convergence rates in the W , ∞ -seminorm for themicrocrack benchmark problem described in Section 5.3. This cor-responds to the case p = ∞ and α = 3 in Table 1. In the above, y a denotes the minimizer of the atomistic energy, and y ac denotesthe corresponding minimizer of the ATM, QCE or BQCE prob-lems. − − − − DoF | E a c − E a | ATMQCEBQCE: smoothBQCE: linear
DoF − DoF − DoF
Figure 4.
Convergence rates in the energy, for the microcrackbenchmark problem described in Section 5.3. In the above, E a denotes the minimum of the atomistic energy, and E ac denotes theminimum of either the ATM, QCE or BQCE energies. − − DoF ! ∇ y a c − ∇ y a ! L ATMQCEBQCE: smoothBQCE: linear
DoF − DoF − DoF
Figure 5.
Convergence rates in the energy-norm (the H -seminorm) for the di-vacancy benchmark problem described in Sec-tion 5.4. This corresponds to the case p = 2 and α = 3 in Table 1.In the above, y a denotes the minimizer of the atomistic energy,and y ac denotes the corresponding minimizer of either the ATM,BQCE or QCE problems.removed from the lattice. We apply 3% isotropic stretch and 3% shear loading, bysetting F := (cid:18) s γ II s (cid:19) · F , where F minimizes W , s = γ II = 0 . W , -seminorm are displayed in Figure 5; the relativeerrors in the W , ∞ -seminorm are displayed in Figure 6; the errors in the energy aredisplayed in Figure 7. For the W , and W , ∞ -seminorms we now observe preciselythe rates of convergence DoF − / and DoF − predicted by our theory. However, weobserve a convergence rate DoF − / for the energy, rather than the square of therate for the W , -seminorm which is DoF − . We can, at present offer no rigorousexplanation for this improved convergence rate, but speculate that it is caused bya cancellation effect that is not captured by our theory.It is important to point out that, in this second numerical experiment, BQCEis not substantially more efficient than ATM. In fact, in the W , ∞ -seminorm theerrors are almost identical. This example clearly demonstrates that an atomisticcoarse-graining scheme need not always lead to improved computational results, butthat the improvement over a judiciously performed atomistic simulation is highlyproblem dependent. We also remark that for similar examples (e.g., a single vacancywith F = F ), we observed that ATM has a higher accuracy than BQCE in boththe W , and W , ∞ -seminorms. ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 17 − − − DoF ! ∇ y a c − ∇ y a ! L ∞ ATMQCEBQCE: smoothBQCE: linear
DoF − DoF − DoF Figure 6.
Convergence rates in the W , ∞ -seminorm for thedi-vacancy benchmark problem described in Section 5.4. This cor-responds to the case p = ∞ and α = 3 in Table 1. In the above, y a denotes the minimizer of the atomistic energy, and y ac denotesthe corresponding minimizer of the ATM, BQCE or QCE prob-lems. − − − − DoF | E a c − E a | ATMQCEBQCE: smoothBQCE: linear
DoF − DoF
DoF − Figure 7.
Convergence rates in the energy, for the di-vacancybenchmark problem described in Section 5.4. In the above, E a denotes the minimum of the atomistic energy, and E ac denotes thecorresponding minimum of either the ATM, QCE or BQCE ener-gies. Conclusion
We have formulated an atomistic-to-continuum coupling, which we call theblended energy-based quasicontinuum method (BQCE). Our formulation requiresonly two approximation parameters for the implementation of BQCE: the blendingfunction β and the finite element mesh T . For the problems of a microcrack anda di-vacancy in a two-dimensional crystal, we utilized theoretical results from [24]to obtain optimal choices of approximation parameters (blending function and fi-nite element grid) to minimize a global error norm and confirmed our analyticalpredictions in numerical tests.An interesting open question is how well different atomistic-to-continuum cou-plings (in particular, QCE, BQCE and the patch test consistent methods) performif the quantity of interest is localized in the atomistic region. In [27], this questionwas investigated through a series of numerical experiments in which the Arlequinmethod was used to simulate a one-dimensional chain with a defect. The authorsconcluded that the local error of the Arlequin method could be controlled by movingthe blending region away from the defect and increasing its size. For the quasicon-tinuum method, local quantities of interest were investigated in [2, 3, 20]. We hopeto develop rigorous error estimates and optimal approximation parameters for lo-calized quantities of interest in a future work.At present BQCE is a competitive choice for numerical simulations of atom-istic multi-scale problems as it applies to a much broader class than any knownpatch test consistent method. For example, Shapeev’s method applies only to two-dimensional, simple lattice crystals with pair interactions. Our BQCE methodcan be applied to challenging and physically important problems featuring three-dimensional, multi-lattice crystals and arbitrary many-body interaction potentials. References [1] A. Abdulle, P. Lin, and A. V. Shapeev. Numerical methods for multilattices. arXiv:1107.3462.[2] Marcel Arndt and Mitchell Luskin. Error estimation and atomistic-continuum adaptivityfor the quasicontinuum approximation of a Frenkel-Kontorova model.
SIAM J. MultiscaleModeling & Simulation , 7:147–170, 2008.[3] Marcel Arndt and Mitchell Luskin. Goal-oriented adaptive mesh refinement for the qua-sicontinuum approximation of a Frenkel-Kontorova model.
Computer Methods in AppliedMechanics and Engineering , 197:4298–4306, 2008.[4] S. Badia, M. Parks, P. Bochev, M. Gunzburger, and R. Lehoucq. On atomistic-to-continuumcoupling by blending.
Multiscale Model. Simul. , 7(1):381–406, 2008.[5] P. T. Bauman, H. B. Dhia, N. Elkhodja, J. T. Oden, and S. Prudhomme. On the applicationof the Arlequin method to the coupling of particle and continuum models.
Comput. Mech. ,42(4):511–530, 2008.[6] T. Belytschko and S. P. Xiao. Coupling methods for continuum model with molecular model.
International Journal for Multiscale Computational Engineering , 1:115–126, 2003.[7] X. Blanc, C. Le Bris, and P.-L. Lions. From molecular models to continuum mechanics.
Arch.Ration. Mech. Anal. , 164(4):341–381, 2002.[8] M. S. Daw and M. S. Baskes. Embedded-Atom Method: Derivation and Application toImpurities, Surfaces, and other Defects in Metals.
Phys. Rev. B
20, 1984.[9] M. Dobson, R. Elliot, M. Luskin, and E. Tadmor. A multilattice quasicontinuum for phasetransforming materials: Cascading cauchy born kinematics.
Journal of Computer-Aided Ma-terials Design , 14:219–237, 2007.[10] L. M. Dupuy, E. B. Tadmor, F. Legoll, R. E. Miller, and W. K. Kim. Finite-temperaturequasicontinuum. manuscript, 2011.
ORMULATION AND OPTIMIZATION OF THE BQCE METHOD 19 [11] W. E, J. Lu, and J. Yang. Uniform accuracy of the quasicontinuum method.
Phys. Rev. B ,74(21):214115, 2006.[12] W. E and P. Ming. Cauchy-Born rule and the stability of crystalline solids: static problems.
Arch. Ration. Mech. Anal. , 183(2):241–297, 2007.[13] H. Fischmeister, H. Exner, M.-H. Poech, S. Kohlhoff, P. Gumbsch, S. Schmauder, L. S.Sigi, and R. Spiegler. Modelling fracture processes in metals and composite materials.
Z.Metallkde. , 80:839–846, 1989.[14] F. C. Frank and J. H. van der Merwe. One-dimensional dislocations. I. static theory.
Proc.R. Soc. London , A198:205–216, 1949.[15] T. Hudson and C. Ortner. On the stability of bravais lattices and their cauchy–born approx-imations.
M2AN Math. Model. Numer. Anal. , 46, 2012.[16] M. Iyer and V. Gavini. A field theoretical approach to the quasi-continuum method.
Journalof the Mechanics and Physics of Solids , 59(8):1506 – 1535, 2011.[17] R. A. Johnson and D. J. Oh. Analytic embedded atom method model for bcc metals.
Journalof Materials Research
4, 1195–1201, 1989.[18] W. K. Kim, E. B. Tadmor, M. Luskin, D. Perez, and A. Voter. Hyper-qc: An acceleratedfinite-temperature quasicontinuum method using hyperdynamics. manuscript, 2011.[19] R. Miller and E. Tadmor. A unified framework and performance benchmark of fourteenmultiscale atomistic/continuum coupling methods.
Modelling Simul. Mater. Sci. Eng. , 17,2009.[20]
J. T. Oden, S. Prudhomme, and P. Bauman , Error control for molecular statics problems ,Int. J. Multiscale Comput. Eng., 4 (2006), pp. 647–662.[21] M. Ortiz, R. Phillips, and E. B. Tadmor. Quasicontinuum analysis of defects in solids.
Philo-sophical Magazine A , 73(6):1529–1563, 1996.[22] C. Ortner and A. Shapeev. work in progress.[23] C. Ortner and A. V. Shapeev. Analysis of an Energy-based Atomistic/Continuum CouplingApproximation of a Vacancy in the 2D Triangular Lattice.
ArXiv e-prints , 1104.0311, Apr.2010.[24] C. Ortner and B. Van Koten. manuscript.[25] C. Ortner and L. Zhang. Construction and sharp consistency estimates for atom-istic/continuum coupling methods with general interfaces: a 2d model problem.arXiv:1110.0168.[26] H. Park, E. Karpov, P. Klein, and W.K.Liu. Three-dimensional bridging scale analysis ofdynamic fracture.
Journal of Computational Physics , 207:588–609, 2005.[27] S. Prudhomme, H. Ben Dhia, P. T. Bauman, N. Elkhodja, and J. T. Oden. Computationalanalysis of modeling error for the coupling of particle and continuum models by the Arlequinmethod.
Comput. Methods Appl. Mech. Engrg. , 197(41-42):3399–3409, 2008.[28] A. Shapeev. Consistent energy-based atomistic/continuum coupling for two-body potentialsin one and two dimensions.
Multiscale Model. Simul. , 9(3):905–932, 2011.[29] A. V. Shapeev. Consistent energy-based atomistic/continuum coupling for two-body poten-tials in three dimensions. arXiv:1108.2991.[30] T. Shimokawa, J. Mortensen, J. Schiotz, and K. Jacobsen. Matching conditions in the qua-sicontinuum method: Removal of the error introduced at the interface between the coarse-grained and fully atomistic region.
Phys. Rev. B , 69(21):214104, 2004.[31] B. Van Koten, X. H. Li, M. Luskin, and C. Ortner. A computational and theoretical inves-tigation of the accuracy of quasicontinuum methods. In I. Graham, T. Hou, O. Lakkis, andR. Scheichl, editors,
Numerical Analysis of Multiscale Problems , volume 83 of
Lect. NotesComput. Sci. Eng.
Springer, 2012. arXiv:1012.6031.[32] B. Van Koten and M. Luskin. Analysis of energy-based blended quasi-continuum approxima-tions.
SIAM Journal on Numerical Analysis , 49(5):2182–2209, 2011.[33] B. Van Koten and C. Ortner. Symmetries of 2-lattices and second order accuracy of theCauchy–Born model. arXiv:1203.5854.[34] R.R. Zope and Y. Mishin. Interatomic potentials for atomistic simulations of the Ti-Al system.
Phys. Rev. B
68, 024102, 2003.
M. Luskin, 127 Vincent Hall, 206 Church St. SE, Minneapolis, MN 55455, USA
E-mail address : [email protected] C. Ortner, Mathematics Institute, Zeeman Building, University of Warwick, Coven-try CV4 7AL, UK
E-mail address : [email protected] B. Van Koten, 127 Vincent Hall, 206 Church St. SE, Minneapolis, MN 55455, USA
E-mail address ::