Goal-Oriented Atomistic-Continuum Adaptivity for the Quasicontinuum Approximation

# Goal-oriented Atomistic-Continuum Adaptivity for the Quasicontinuum Approximation

Marcel Arndt  and  Mitchell Luskin
April 23, 2019
###### Abstract.

We give a goal-oriented a posteriori error estimator for the atomistic-continuum modeling error in the quasicontinuum method, and we use this estimator to design an adaptive algorithm to compute a quantity of interest to a given tolerance by using a nearly minimal number of atomistic degrees of freedom. We present computational results that demonstrate the effectiveness of our algorithm for a periodic array of dislocations described by a Frenkel-Kontorova type model.

###### Key words and phrases:
quasicontinuum, atomistic-continuum, error estimation, a posteriori, goal-oriented, Frenkel-Kontorova, dislocation, defect
###### 2000 Mathematics Subject Classification:
65Z05, 70C20, 70G75
This work was supported in part by DMS-0304326 and by the Minnesota Supercomputing Institute. This work is also based on work supported by the Department of Energy under Award Number DE-FG02-05ER25706.

## 1. Introduction

Multiscale methods offer the potential to solve complex problems by utilizing a fine-scale model only in regions that require increased accuracy. However, it is usually not known a priori which region requires the increased accuracy, and an adaptive algorithm is needed to compute a given quantity of interest to a given tolerance with nearly optimal computational efficiency.

The quasicontinuum (QC) method [9, 10, 11, 4] is a multiscale computational method for crystals that retains sufficient accuracy of the atomistic model by utilizing a strain energy density obtained from the atomistic model by the Cauchy-Born rule in regions where the deformation is nearly uniform. The atomistic model is needed to accurately model the material behavior in regions of highly non-uniform deformations around defects such as dislocations and cracks.

The approximation error within the quasicontinuum method can be decomposed into the modeling error which occurs when replacing the atomistic model by a continuum model, and the coarsening error which arises from reducing the number of degrees of freedom within the continuum region. This paper purely focuses on the estimation of the modeling error. The optimal strategy to determine the mesh size in the continuum region will be studied in a forthcoming paper.

The development of goal-oriented error estimators for mesh coarsening in the quasicontinuum method has been given in [7, 6], and goal-oriented error estimators for atomistic-continuum modeling have recently been given in [2]. In all these works, the error is measured in terms of a user-definable quantity of interest instead of a global norm. Goal-oriented error estimation in general is based on duality techniques and has already been used for a variety of applications such as mesh refinement in finite element methods [1, 3] and control of the modeling error in homogenization [8].

In [2], an a posteriori error estimator for the modeling error in the quasicontinuum method has been developed, analyzed, and applied to a one-dimensional Frenkel-Kontorova model of crystallographic defects [5]. In this paper, we summarize this approach and adapt it to a different setting. Instead of clamped boundary conditions, we use periodic boundary conditions here which are physically more realistic and allow for more succinct formulas. Moreover, an asymmetric quantity of interest is used here rather than the symmetric one studied in [2] to give further insight into the behavior of the error estimator.

The one-dimensional periodic Frenkel-Kontorova model chosen here allows for an easy study of the error estimator and keeps the formulas short, but still exhibits enough of the features of multidimensional problems for a realistic study of the error estimator. In addition to the nearest-neighbor harmonic interactions between the atoms in the classical Frenkel-Kontorova model, we add next-nearest-neighbor harmonic interactions to obtain a non-trivial quasicontinuum approximation.

We describe the atomistic model and its quasicontinuum approximation in Section 2, and we formulate the error estimator in Section 3. We then develop in Section 3 an algorithm which employs the error estimator for an adaptive strategy, and we conclude by exhibiting and interpreting the results of our numerical experiments.

## 2. Atomistic Model and Quasicontinuum Model

As an application for the method of error estimation described in this paper, we study a periodic array of dislocations in a single crystal. We employ a Frenkel-Kontorova type model [5] to give a simplified one-dimensional description of these typically two-dimensional or three-dimensional phenomena. To model the elastic interactions, we consider atoms in a periodic chain that interact by Hookean nearest-neighbor and next-nearest-neighbor springs. The dislocation is modeled by the interaction with a substrate which gives rise to a misfit energy, see Figure 1.

The vector describes the positions of atoms which generate the positions of an infinite chain of atoms by the relation

 yi+2M=yi+(2M+1)a0for i=−∞,…,∞, (2.1)

where denotes the lattice constant. The relation (2.1) gives an average strain of due to a periodic array of dislocations that stretch the chain by one lattice constant every atoms.

The total energy for this atomistic system reads

 Ea(y)=Ea,e(y)+Ea,m(y) (2.2)

with elastic energy given by (recall (2.1))

 Ea,e(y)=12k1M∑i=−M+1(yi+1−yi−a0)2+12k2M+1∑i=−M+2(yi+1−yi−1−2a0)2=12k1[M−1∑i=−M+1(yi+1−yi−a0)2+(y−M+1−yM+2Ma0)2]+12k2[M−1∑i=−M+2(yi+1−yi−1−2a0)2+(y−M+1−yM−1+(2M−1)a0)2+(y−M+2−yM+(2M−1)a0)2] (2.3)

and misfit energy given by

 Ea,m(y) =12k0M∑i=−M+1(yi−a0⌊yia0+12⌋)2. (2.4)

Here , and denote the elastic constants. In the misfit energy, denotes the largest integer smaller than or equal to We can obtain the following more symmetric form of the elastic energy (2.3) by realizing that the forces of constraint corresponding to the strain (2.1) move the equilibrium spacing of the chain to

 Ea,e(y)=12k1[M∑i=−M+1(yi+1−yi−2M+12M(((i+1) mod 2M)−i)a0)2+12Ma20]+12k2[M∑i=−M+1(yi+2−yi−2M+12M(((i+2) mod 2M)−i)a0)2+2Ma20] (2.5)

where all atom indices are understood modulo , and denotes modulo within the interval . In the following, we neglect the constant terms since they do not have any effect when finding energy minimizers later.

We consider a single vacancy between the atoms and . If we assume that the leftmost atoms for are situated in the interval and that the rightmost atoms for are situated in the interval , then the misfit energy can be expressed more simply as

 Ea,m(y) =12k0M∑i=−M+1(yi−bi)2withbi={(i−1)a0for i≤0,ia0for i≥1. (2.6)

To reduce the work in computing (2.2), we employ the quasicontinuum method [9, 10, 11] which has been successfully used for many applications. The quasicontinuum method consists of two steps: the passage to a continuum energy within the continuum region of the chain, and a subsequent coarsening in the continuum region to reduce the number of degrees of freedom.

In the first step, we replace the atomistic energy of all atoms from the continuum region by a continuum energy. To this end, we split the total energy into atom-wise contributions:

 Ea(y)=M∑i=−M+1Ea,i(y) (2.7)

with

 (2.8)

The corresponding continuum energy can be derived following [2] to be

 (2.9)

If

 δai={1if atom i is modeled % atomistically,0if atom i is modeled as continuum,andδci=1−δai, (2.10)

then

 Eac(y)=M∑i=−M+1[δaiEa,i(y)+δciEc,i(y)] (2.11)

denotes the mixed atomistic-continuum energy.

In the second step of the quasicontinuum approximation, the chain is coarsened in the continuum region by choosing representative atoms, more briefly called repatoms. The chain is then fully modeled in terms of the repatoms. The missing atoms are implicitly reconstructed by linear interpolation according to the Cauchy-Born hypothesis. The lengthy expression of the resulting quasicontinuum energy

 Eqc(y) (2.12)

is not needed in this paper since we focus on the estimation of the modeling error. Hence we refer to [2] for the formula and its derivation. The error arising from coarsening will be studied in a forthcoming paper.

For the subsequent argumentation, it is useful to rewrite the energies in matrix notation. We have

where the matrices are given by

 Di,i=−1,Di,i+1=1,(Ea)i,i=k1+2k2,(Ea)i,i+1=(Ea)i+1,i=k2,\omit\span\omit(Eac)i,i=12k1(δai+δai+1)+12k2(δai−1+δai+δai+1+δai+2)+(12k1+2k2)(δci+δci+1),\omit\span\omit(Eac)i,i+1=(Eac)i+1,i=12k2(δai+δai+2),(Ka)i,i=k0, (2.14)

with and all indices to be understood modulo as before. The vectors and are defined as

 (2.15)

We require that the elastic moduli satisfy to ensure that is positive definite and that the misfit modulus to ensure that is positive definite.

We are interested in finding energy minimizing configurations , , and of , , and , respectively. The minimizers and satisfy the linear equations

 Maya=fa,Macyac=fac, (2.16)

where

We refer to (2.16) as the primal problems. Note that the minimizers are uniquely determined due to the convexity of the energy.

## 3. Error Estimation

In the preceding section, we described how the quasicontinuum method gives an approximation of the atomistic solution by passing from the fully atomistic model to a mixed atomistic-continuum formulation, and then we briefly mentioned how a further approximation, , can be obtained by coarsening in the continuum region.

Instead of measuring the error in some global norm, we measure the error of a quantity of interest denoted by for some function . We assume that is linear and thus has a representation

 Q(y)=qTy (3.1)

for some vector . We then have the splitting

 |Q(ya)−Q(yqc)|=|Q(ya−yac)+Q(yac−yqc)|≤|Q(ya−yac)|+|Q(yac−yqc)| (3.2)

of the total error into the modeling error, , and the coarsening error, , everything in terms of the quantity of interest. In this paper, we restrict ourselves to the estimation of the modeling error. The coarsening error will be analyzed in a forthcoming paper.

An important tool for the error estimation in terms of a quantity of interest are the dual problems for the influence or generalized Green’s functions and given by

 Maga=q,Macgac=q. (3.3)

The matrices and are symmetric since they stem from an energy, and we thus do not need to use their transpose for the dual problems.

We denote the errors and the residuals, both for the deformation and the influence function by

 e :=ya−yac, R(y) :=Ma(ya−y)=fa−May, (3.4) ^e :=ga−gac, ^R(g) :=Ma(ga−g)=q−Mag.

Then we have the basic identity for the error of the quantity of interest

 Q(ya)−Q(yac)=qTe=gaTMae=(gac+^e)TMae=gacTR(yac)+^eTMae. (3.5)

The quantities and are considered to be computable since the continuum degrees of freedom give local interactions, whereas and are not considered to be computable since they require a full atomistic computation. Thus the first term is easily computable, and the challenge is to estimate . Let us note that in applications to mesh refinement for linear finite elements, the residual term vanishes due to Galerkin orthogonality, whereas in other applications it can be dominant over the second term.

We utilize two error estimators derived in [2] and briefly summarized here. Our first error estimator is based on the generalized parallelogram identity

 ^eTMae=14∥σe+σ−1^e∥2Ma−14∥σe−σ−1^e∥2Ma (3.6)

for all , where the -norm of some vector is defined by . We define the computable bounds

 η±low≤∥σe±σ−1^e∥Ma≤η±upp (3.7)

by

 η±upp :=∥∥PD[σ(yac−aa)±σ−1gac]∥∥Ea, (3.8) η±low :=∣∣(yac+θ±gac)Tr±∣∣∥yac+θ±gac∥Ma

where

 P :=I−(Ea)−1Eac, (3.9) r± :=σR(yac)±σ−1^R(gac).

Optimization of the bounds with respect to and leads in [2] to the following choice of the parameters:

 σ:=√∥PDgac∥Ea∥PD(yac−aa)∥Ea,θ±:=r±TyacgacTMayac−r±Tgac∥yac∥2Mar±TgacgacTMayac−r±Tyac∥gac∥2Ma. (3.10)
###### Theorem 3.1.

We have that

 |Q(ya)−Q(yac)|≤η1, (3.11)

where the computable error estimator is defined as

 η1:=max(∣∣gacTR(yac)+14(η+low)2−14(η−upp)2∣∣,∣∣gacTR(yac)+14(η+% upp)2−14(η−low)2∣∣). (3.12)

We also developed the following weaker estimator in [2] using the Cauchy-Schwarz inequality in place of the parallelogram identity in (3.6). We note that this estimator can be decomposed among the degrees of freedom and can thus be utilized in adaptive atomistic-continuum modeling decisions.

###### Theorem 3.2.

We have that

 |Q(ya)−Q(yac)|≤η2≤M∑i=−M+1ηat2,i+M∑i=−M+1ηel2,i (3.13)

where the computable global error estimator and the computable local error estimators and , associated with atoms and elements, respectively, are defined as

η2 ηat2,i ηel2,i :=∣∣gacTR(yac)∣∣+∥PD(yac−aa)∥Ea∥PDgac∥Ea, :=∣∣gaciR(yac)i∣∣, i=−M+1,…,M, +12∣∣(PDgac)i((Ea−Eac)Dgac)i∣∣, i=−M+1,…,M.
(3.14)

## 4. Numerics

Now we use the two a posteriori error estimators given in Section 3 to formulate an algorithm which adaptively decides between atomistic and continuum modeling. Finally, we present and discuss the numerical results for the periodic array of dislocations described by the Frenkel-Kontorova model.

The error estimator gives a better estimate of the error because involves additional inequalities. However, allows for an atom-wise decomposition, whereas does not. This is due to the fact that in the definition of is equal to the square of a sum of atom-wise components and not the sum of the square of these components. We can thus let decide whether a given global tolerance is already achieved or not and use the decomposition of to decide where the atomistic model is needed for a better approximation. In this way, we combine the better efficiency of with the error localization of .

We start with a fully continuum model. We then switch to the atomistic model wherever the local error exceeds an atom-wise error tolerance . While decreasing , the algorithm adaptively tags larger and larger regions atomistic until the estimate for the goal-oriented error finally reaches . The complete algorithm reads:

1. Choose . Model all atoms as a continuum. Set .

2. Solve primal problem (2.16) for and dual problem (3.3) for .

3. Compute error estimator from (3.12).

4. If then stop.

5. Compute local error estimators and from (3.14).

6. Set .

7. Make all atoms atomistic ( and ) for which

 ηtot2,i:=ηat2,i+12(ηel2,i−1+ηel2,i)≥τat. (4.1)
8. Go to (2).

The factor describes the rate at which the atom-wise tolerance is decreased during the adaptive process. We found that gives an efficient method for this problem.

Now we come to the results for the Frenkel-Kontorova dislocation model for a periodic chain of 1000 atoms, that is . The elastic constants are set to be and . To define the quantity of interest, we choose the average displacement of atoms . This leads to

 q=(qi)i=−M+1,…,M,qi=1for% 11≤i≤30,qi=0otherwise. (4.2)

The global tolerance is chosen to be .

Table 1 shows how the successive adaptive determination of the atomistic-continuum modeling proceeds. After six iterations, the atom-wise tolerance is small enough so that , that is the desired accuracy has been achieved.

Figure 2 shows the decomposition of the error estimator for the fully continuum model in iteration 1 of the adaption process. One can clearly read off from the graph that the error is large near the dislocation between atoms 0 and 1 and near atoms 11 and 30, and that it decays exponentially away from these points. We note the slight nonsymmetry of the atomistic-continuum modeling due to using a goal function which averages over atoms to the right of the dislocation, but not to its left. The graph on the right shows the decomposition of in the final iteration 6 with an atomistic region given by indices . It exhibits the same nonsymmetry, but the error is considerably smaller with peaks at the boundary between the atomistic region and the continuum region. In both diagrams, the fluctuations come from the limited relative machine precision of about .

Finally, Table 2 shows the efficiency of the error estimators and for different atomistic regions. gives the actual error which can be computed for this relatively small problem. In real applications, it is of course not available. One can clearly see that gives a better estimate than , which numerically confirms our conjecture that is a better estimator than . An unusually high value for the efficiency occurs when the atomistic-continuum boundary sweeps through the region where the quantity of interest is measured. After this, the efficiencies converge to decent values around 1.7 and 2.5 for and , respectively. We note that for clamped boundary conditions and a symmetric quantity of interest, better efficiencies of 1.4 and 2, respectively, have been obtained [2].

## References

• [1] M. Ainsworth and J. T. Oden, A Posteriori Error Estimation in Finite Element Analysis, Wiley, New York, 2000.
• [2] M. Arndt and M. Luskin, Error estimation and atomistic-continuum adaptivity for the quasicontinuum approximation of a Frenkel-Kontorova model, 2007, arXiv:0704.1924.
• [3] W. Bangerth and R. Rannacher, Adaptive Finite Element Methods for Differential Equations, Lectures in Mathematics, ETH Zürich, Birkhäuser, Basel, 2003.
• [4] M. Dobson and M. Luskin, Analysis of a force-based quasicontinuum approximation, 2006, arXiv:math.NA/0611543.
• [5] M. Marder, Condensed Matter Physics, John Wiley & Sons, 2000.
• [6] J. T. Oden, S. Prudhomme, and P. Bauman, Error control for molecular statics problems, Int. J. Multiscale Comput. Eng., 4 (2006), pp. 647–662.
• [7] J. T. Oden, S. Prudhomme, A. Romkes, and P. Bauman, Multiscale modeling of physical phenomena: Adaptive control of models, SIAM J. Sci. Comput., 28 (2006), pp. 2359–2389.
• [8] J. T. Oden and K. S. Vemaganti, Estimation of local modeling error and goal-oriented adaptive modeling of heterogeneous materials: Part I: Error estimates and adaptive algorithms, J. Comput. Phys., 164 (2000), pp. 22–47.
• [9] E. B. Tadmor, R. Miller, R. Phillips, and M. Ortiz, Nanoindentation and incipient plasticity, J. Mater. Res., 14 (1999), pp. 2233–2250.
• [10] E. B. Tadmor, M. Ortiz, and R. Phillips, Quasicontinuum analysis of defects in solids, Philos. Mag. A, 73 (1996), pp. 1529–1563.
• [11] E. B. Tadmor, R. Phillips, and M. Ortiz, Mixed atomistic and continuum models of deformation in solids, Langmuir, 12 (1996), pp. 4529–4534.
You are adding the first comment!
How to quickly get a good reply:
• Give credit where it’s due by listing out the positive aspects of a paper before getting into which changes should be made.
• Be specific in your critique, and provide supporting evidence with appropriate references to substantiate general statements.
• Your comment should inspire ideas to flow and help the author improves the paper.

The better we are at sharing our knowledge with each other, the faster we move forward.
The feedback must be of minimum 40 characters and the title a minimum of 5 characters