Goaloriented AtomisticContinuum Adaptivity for the Quasicontinuum Approximation
Abstract.
We give a goaloriented a posteriori error estimator for the atomisticcontinuum 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 FrenkelKontorova type model.
Key words and phrases:
quasicontinuum, atomisticcontinuum, error estimation, a posteriori, goaloriented, FrenkelKontorova, dislocation, defect2000 Mathematics Subject Classification:
65Z05, 70C20, 70G751. Introduction
Multiscale methods offer the potential to solve complex problems by utilizing a finescale 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 CauchyBorn rule in regions where the deformation is nearly uniform. The atomistic model is needed to accurately model the material behavior in regions of highly nonuniform 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 goaloriented error estimators for mesh coarsening in the quasicontinuum method has been given in [7, 6], and goaloriented error estimators for atomisticcontinuum modeling have recently been given in [2]. In all these works, the error is measured in terms of a userdefinable quantity of interest instead of a global norm. Goaloriented 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 onedimensional FrenkelKontorova 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 onedimensional periodic FrenkelKontorova 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 nearestneighbor harmonic interactions between the atoms in the classical FrenkelKontorova model, we add nextnearestneighbor harmonic interactions to obtain a nontrivial 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 FrenkelKontorova type model [5] to give a simplified onedimensional description of these typically twodimensional or threedimensional phenomena. To model the elastic interactions, we consider atoms in a periodic chain that interact by Hookean nearestneighbor and nextnearestneighbor 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
(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
(2.2) 
with elastic energy given by (recall (2.1))
(2.3) 
and misfit energy given by
(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
(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
(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 atomwise contributions:
(2.7) 
with
(2.8) 
The corresponding continuum energy can be derived following [2] to be
(2.9) 
If
(2.10) 
then
(2.11) 
denotes the mixed atomisticcontinuum 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 CauchyBorn hypothesis. The lengthy expression of the resulting quasicontinuum energy
(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
(2.13) 
where the matrices are given by
(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
(2.16) 
where
(2.17)  
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 atomisticcontinuum 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
(3.1) 
for some vector . We then have the splitting
(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
(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
(3.4)  
Then we have the basic identity for the error of the quantity of interest
(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
(3.6) 
for all , where the norm of some vector is defined by . We define the computable bounds
(3.7) 
by
(3.8)  
where
(3.9)  
Optimization of the bounds with respect to and leads in [2] to the following choice of the parameters:
(3.10) 
Theorem 3.1.
We have that
(3.11) 
where the computable error estimator is defined as
(3.12) 
We also developed the following weaker estimator in [2] using the CauchySchwarz 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 atomisticcontinuum modeling decisions.
Theorem 3.2.
We have that
(3.13) 
where the computable global error estimator and the computable local error estimators and , associated with atoms and elements, respectively, are defined as

(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 FrenkelKontorova model.
The error estimator gives a better estimate of the error because involves additional inequalities. However, allows for an atomwise decomposition, whereas does not. This is due to the fact that in the definition of is equal to the square of a sum of atomwise 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 atomwise error tolerance . While decreasing , the algorithm adaptively tags larger and larger regions atomistic until the estimate for the goaloriented error finally reaches . The complete algorithm reads:

Choose . Model all atoms as a continuum. Set .

Compute error estimator from (3.12).

If then stop.

Compute local error estimators and from (3.14).

Set .

Make all atoms atomistic ( and ) for which
(4.1) 
Go to (2).
The factor describes the rate at which the atomwise 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 FrenkelKontorova 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
(4.2) 
The global tolerance is chosen to be .
iteration  atomistic region  

1  none  1.000000e10  6.860546e03 
2  1.000000e11  1.238016e07  
3  1.000000e12  2.600112e08  
4  1.000000e13  3.922946e09  
5  1.000000e14  4.104868e10  
6  1.000000e15  4.105166e11 
Table 1 shows how the successive adaptive determination of the atomisticcontinuum modeling proceeds. After six iterations, the atomwise 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 atomisticcontinuum 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 .
atomistic  

region  
none  1.416421e03  6.860545e03  4.843577  1.231314e02  8.693133  
4  10  1.863104e03  6.107510e03  3.278136  1.049800e02  5.634680  
9  20  1.000572e05  3.358722e04  33.56803  6.621488e04  66.17705  
14  30  1.430363e04  3.187552e04  2.228492  5.140285e04  3.593694  
19  40  1.675490e05  2.626711e05  1.567727  3.691344e05  2.203142  
24  50  7.361419e07  1.190138e06  1.616723  1.693910e06  2.301065  
29  60  3.139276e08  5.157753e08  1.642975  7.388556e08  2.353586  
34  70  1.146997e09  2.001550e09  1.745035  2.934377e09  2.558312 
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 atomisticcontinuum 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 atomisticcontinuum adaptivity for the quasicontinuum approximation of a FrenkelKontorova 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 forcebased 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 goaloriented 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.