A/C coupling with higher-order finite elements

Analysis of patch-test consistent atomistic-to-continuum coupling with higher-order finite elements


We formulate a patch test consistent atomistic-to-continuum coupling (a/c) scheme that employs a second-order (potentially higher-order) finite element method in the material bulk. We prove a sharp error estimate in the energy-norm, which demonstrates that this scheme is (quasi-)optimal amongst energy-based sharp-interface a/c schemes that employ the Cauchy–Born continuum model. Our analysis also shows that employing a higher-order continuum discretisation does not yield qualitative improvements to the rate of convergence.

Key words and phrases:
atomistic models, coarse graining, atomistic-to-continuum coupling, quasicontinuum method, error analysis
HW was supported by MASDOC doctoral training centre, EPSRC grant EP/H023364/1. CO was supported by ERC Starting Grant 335120.

1. Introduction

Atomistic-to-continuum (a/c) coupling is a class of coarse-graining methods for efficient atomistic simulations of systems that couple localised atomistic effects described by molecular mechanics with long-range elastic effects described by continuum models using the finite-element method. We refer to [6], and references therein for an extensive introduction and references.

The presented work explores the feasibility and effectiveness of introducing higher-order finite element methods in the a/c framework, specifically for quasi-nonlocal (QNL) type methods.

The QNL-type coupling, first introduced in [15], is an a/c method that uses a “geometric consistency condition” [3] to construct the coupling between the atomistic and continuum model. The first explicit construction of such a scheme for two-dimensional domains with corners is described in [12] for a neareast-neighbour many-body site potential. We call this construction ”G23” for future reference. This approach satisfies force and energy patch tests (often simply called consistency), which in particular imply absence of ghost forces.

We will supply the G23 scheme with finite element methods of different orders and investigate the rates of convergence for the resulting schemes. Our conclusion will be that second-order finite element schemes are theoretically superior to first-order schemes, while schemes of third and higher order do not improve the rate of convergence. This is due to the fact that the consistency error of the a/c scheme is dominated by the modelling error committed at the a/c interface. We will also explore, for some basic model problems, how well second-order schemes fare in practise against first-order schemes.

1.1. Outline

The theory of high-order finite element methods (FEM) in partial differential equations, and applications in solid mechanics is well established; see [14] and references therein. However, most work on the rigorous error analysis of a/c coupling has been restricted to P1 finite element methods; the only exception we are aware of is [11], which focuses on blending-type methods.

In the present work we estimate the accuracy of a QNL method employing a P2 FEM in the continuum region against an exact solution obtain from a fully atomistic model. Since stability of QNL type couplings is a subtle issue [9] we will primarily analyse the consistency errors, taking account the relative sizes of the fully resolved atomistic region and of the entire computational domain (Sections 5.1-5.4). We will then optimize these relative sizes as well as the mesh grading in the continuum region in order to minimize the total consistency error (Section 5.5). We will observe that, using P1-FEM in the continuum region, the error resulting from FEM approximations is the dominating contributor of the consistency estimates, which implies that increasing the order of the FEM can indeed improve the accuracy of the simulation. We will show that, using Pk-FEM with , the FEM approximation error is dominated by the interface error which comes purely from the G23 construction, and in particular demonstrate that the P2-FEM is sufficient to achieve the optimal convergence rate for the consistency error. Finally, assuming the stability of G23 coupling (see [9] why this must be an assumption), we prove a rigorous error estimate in §6.

Finally, we conduct numerical experiments on a 2D anti-plane model problem to test our analytical predictions. The numerical results display the predicted error convergence rates for the fully atomistic model, P1-FEM G23 model, and P2-FEM G23 model.

2. Preliminaries

Our setup and notation follows [12]. As our model geometry we consider an infinite 2D triangular lattice,

We define the six nearest-neighbour lattice directions by , and , where denotes the rotation through the angle . We supply with an atomistic triangulation, as shown in Figure 1, which will be convenient in both analysis and numerical simulations. We denote this triangulation by and its elements by . We also denote , and , for .

Figure 1. The lattice and its canonical triangulation.

We identify a discrete displacement map , , with its continuous piecewise affine interpolant, with weak derivative , which is also the pointwise derivative on each element . For , the spaces of displacements are defined as

We equip with the -seminorm, . From [8] we know that is dense in in the sense that, if , then there exist such that strongly in .

A homogeneous displacement is a map , where .

For a map , we define the finite difference operator


Note that .

2.1. 2D many-body nearest neighbour interactions

We assume that the atomistic interaction is described by a nearest-neighbour many-body site energy potential ,, with . Furthermore, we assume that satisfies the point symmetry

The energy of a displacement , given by

is well-defined since the infinite sum becomes finite. To formulate a variational problem in the energy space , we need the following lemma to extend to .


is continuous and has a unique continuous extension to , which we still denote by . Furthermore, the extended is -times continuously Fréchet differentiable.


See Lemma 2.1 in [4]. ∎

We model a point defect by adding an external potential with for all , where is some given radius (the defect core radius), and for all constants . For example, we can think of modelling a substitutional impurity. See also [5, 7] for similar approaches.

We then seek the solution to


For we define the first and second variations of by

We use analogous definitions for all energy functionals introduced in later sections.

2.2. The Cauchy–Born Approximation

The Cauchy–Born strain energy function, corresponding to the interatomic potential is

where is the volume of a unit cell of the lattice . Thus is the energy per volume of the homogeneous lattice .

2.3. The G23 coupling method

Let denote the set of all lattices sites for which we want to maintain full atomistic accuracy. We denote the set of interface lattice sites by

and we denote the remaining lattice sites by . Let be the Voronoi cell associated with site . We define the continuum region ; see Figure 2.

Figure 2. The domain decomposition with a layer of interface atoms.

A general form for the GRAC-type a/c coupling energy [3, 12] is

where . The parameters are to be determined in order for the coupling scheme to satisfy the “patch tests”:

is locally energy consistent if, for all ,


is force consistent if, for all ,


is patch test consistent if it satisfies both (2.3) and (2.4).

For the sake of brevity of notation we will often write . Following [12] we make the following standing assumption (see Figure 3 for examples).

(A0) Each vertex has exactly two neighbours in , and at least one neighbour in .

Figure 3. The first two configurations are allowed. The third configuration is not allowed as the interface atom at the corner has no nearest neighbour in the continuum region, and should instead be taken as an atomistic site.

Under this assumption, the geometry reconstruction operator is then defined by

see Figure 4. The resulting coupling method is called G23 and the corresponding energy functional . This choice of coefficients (and only this choice) leads to patch test consistency (2.3) and (2.4).

Figure 4. The geometry reconstruction coefficents at the interface sites.

For future reference we decompose the canonical triangulation as follows:


2.4. Notation for a P finite element scheme

In the atomistic and interface regions, the interactions are represented by discrete displacement maps, which are identified with their linear interpolant. Here, we identify the displacement map with its P1 interpolant. No approximation error is committed.

On the other hand, in the continuum region where the interactions are approximated by the Cauchy–Born energy, we could increase the accuracy by using P-FEM with . In later sections we will review that the Cauchy–Born approximation yields a 2nd-order error, whereas employing the P1-FEM in the continuum region would reduce the accuracy to first order. In fact, we will show in that, with optimized mesh grading, P2-FEM is sufficient to obtain a convergence rate that cannot be improved by other choices of continuum discretisations. High-order P-FEM with will increase the computational costs but yield the same error convergence rate (see § 3.5).

Let denote the inner radius of the atomistic region,

where denotes the ball of radius centred at . In order for the defect to be contained in the atomistic region we assume throughout that .

Let denote the entire computational domain and denote the inner radius of , i.e.,

Let be a finite element triangulation of which satisfies, for ,

In other words, and coincide in the atomistic and interface regions, whereas in the continuum region the mesh size may increase towards the domain boundary. The optimal rate at which the mesh size increases will be determined in later sections.

We note that the concrete construction of will be based on the choice of the domain parameters and ; hence, when emphasizing this dependence, we will write . We assume throughout that the family is uniformly shape-regular, i.e., there exists such that,


This assumption eliminates the possibility of extreme angles on elements, which would deteriorate the constants in finite element interpolation error estimates. For the most part we will again drop the parameters from the notation by writing but implicitly will always keep the dependence.

Similar to (2.5), we define the atomistic, interface and continuum elements as and , respectively. Note that and . We also let denote the number of degrees of freedom of .

We define the finite element space of admissible displacements as


In defining we have made two approximations to the class of admissible displacements: (1) truncation to a finite computational domain and (2) finite element coarse-graining.

The computational scheme is to find


is embedded in via point evaluation. Through this identification, is well-defined for all .

We will make this identification only when we evaluate . The reason for this is a conflict when interpreting elements as lattice functions is that we identify lattice functions with their continuous interpolants with respect to the canonical triangulation , which would be different from the function itself. However, for the evaluation of this issue does not arise.

3. Summary of results

3.1. Regularity of

The approximation error analysis in later sections requires estimates on the decay of the elastic fields away from the defect core. These results follow from a natural stability assumption:

(A1) The atomistic solution is strongly stable, that is, there exists ,


where is a solution to (2.2).


Suppose that (A1) is satisfied, then there exists a constant such that, for ,


See Theorem 2.3 in [4]. ∎

3.2. Stability

In [9] it is shown that there is a “universal” instability in 2D interfaces for QNL-type a/c couplings: it is impossible to prove in full generality that is a positive definite operator, even if we assume (3.1). Indeed, this potential instability is universal to a wide class of generalized geometric reconstruction methods. However, it is rarely observed in practice. To circumvent this difficulty, we make the following standing assumption:

(A2) The homogeneous lattice is strongly stable under the G23 approximation, that is, there exists which is independent of such that, for sufficiently large,


Since (3.2) does not depend on the solution it can be tested numerically. But a precise understanding under which conditions (3.2) is satisfied is still missing. In [9] a method of stabilizing 2D QNL-type schemes with flat interfaces is introduced, which could replace this assumption, however we are not yet able to extend this stabilizing method for interfaces with corners, such as the configurations discussed in this paper.

3.3. Main results

To state the main results it is convenient to employ a smooth interpolant to measure the regularity of lattice functions. In Lemma 6.1, we define such an interpolant for , for which there exists a universal constant such that, for all , ,

where .

Consistency error estimate

In (5.6) we define a quasi-best approximation operator , which truncates an atomistic displacement to enforce the homogeneous Dirichlet boundary condition, and then interpolates it onto the finite element mesh.

Our main result is the following consistency error estimate.


If is a solution to (2.2) then we have, for all ,


where corresponds to the continuum region of , and with .

Optimizing the approximation parameters

Before we estimate the error , we optimize the approximation parameters in the computational scheme. This means that the radius of the atomistic region, the radius of the entire computational domain and the mesh size should satisfy certain balancing relations. We only outline the result of this optimisation and refer to § 6.7 for the details.

Due to the decay estimates on the dominating terms in (3.3) turn out to be


(We will see momentarily that the mesh size plays a minor role.) These two terms result from the nature of the coupling scheme and the far-field truncation error. In particular, both of these cannot be improved by the choice of discretisation of the Cauchy–Born model, e.g., order of the FEM. We also note that, if we had employed a P1-FEM, then the limiting factor would have been .

We can balance the two terms in (3.4) by choosing . It then remains to determine a mesh-size so that the finite element error contribution,

remains small in comparison. We show that the scaling is a suitable choice, with , under which both terms become of order .

Thus, we have determined the approximation parameters in terms of a single parameter . The quasi-optimal relations for P2-FEM discretisatino of the Cauchy–Born model are summarised in Table 1.

consistency error
Table 1. Quasi-optimal relations between approximation parameters for P2-GR23 and, for comparision, for P1-GR23.

Suppose that satisfy the relations of Table 1, the consistency error estimate (3.3) in terms of the number of degrees of freedom can be written as


Error estimate

To complete our summary of results, we now use the Inverse Function Theorem to obtain error estimates for the strains and the energy.


Suppose that (A0), (A1) and (A2) are satisfied and that the quasi-optimal scaling of from Table 1 is satisfied. Then, for sufficiently large atomistic region size , a solution to (2.8) exists which satisfies the error estimates


where is the number of degrees of freedom.

3.4. Setup of the numerical tests

For our numerical tests, we consider an anti-plane displacement . We choose a hexagonal atomistic region with side length and one layer of atomistic sites outside as the interface. To construct the finite element mesh, we add hexagonal layers of elements such that, for each layer , , with ; see Figure 5. The procedure is terminated once the radius of the domain exceeds . This construction guarantees the quasi-optimal approximation parameter balance to optimise the P2-FEM error. The derivation is given in Section 6.7.

In our tests we compare the P2-G23 method against

  • a pure atomistic model with clamped boundary condition: the construction of the domain is as in the P2-G23 method, but without continuum region;

  • a P1-G23 method: the construction is again identical to that of the P2-G23 method, but the P2-FEM in the definition of is replaced by a P1-FEM. The same mesh scaling as for P2 is used (see also [3] where this is shown to be quasi-optimal).

The site potential is given by a nearest-neighbour embedded atom toy model,

with and . This is the anti-plane toy model as the one used in [4].

The external potential is defined by , which can be thought of as an elastic di-pole. A steepest descent method, preconditioned with a finite element Laplacian and fixed (manually tuned) step-size, is used to find a minimizer of , using as the starting guess.

In order to compare the errors, we use a comparison solution with atomistic region size and other computational parameters scaled as above.

The numerical results, with brief discussions, are shown in Figures 69. The two most important observations are the following:

  • the numerical tests confirm the analytical predictions for the energy-norm error, but the experimental rates for the energy error are better than the analytical rates. Similar observations were also made in [4].

  • With our specific setup, the improvement of the P2-GR23 over P1-GR23 is clearly observed when plotting the error against , but when plotted against the improvement is only seen in the asymptotic regime. This indicates that further work is required, such as a posteriori adaption, to optimise the P2-GR23 in the pre-asymptotic regime as well.

Figure 5. An example of the computatioanl mesh. The the vertices marked by ”•” are the atomistic sites; the vertices marked by ”” are the interface sites.
Figure 6. Error in energy norm plotted against . We clearly observe the predicted rate of convergence.
Figure 7. Error in energy norm plotted against the number of degrees of freedom. The improvement of P2-FEM is now only seen asymptotically.
Figure 8. The energy error plotted against . The observed rate of convergence is better than the rate predicted in Theorem 3.3.3.
Figure 9. The energy error plotted against the number of degrees of freedom. The improvement of P2-FEM over P1-FEM can again only be observed asymptotically.

3.5. Extension to high-order FEM

If we apply higher-order FEM in the continuum region, then to extend our error analysis we would need a smooth interpolant of with higher regularity than . A suitable extension given in [5] is, for arbitrary , a piecewise polynomial of degree with properties analogous to those stated in Lemma 6.1. The resulting higher-order decay rate indicates that the use of high-order FEM could be beneficial.

However, as we have pointed out in § 3.3.2, if we employ the mesh grading with in the continuum region, the total approximation error cannot be improved by using P-FEM with , since the dominating term is the interface error for , which results from the construction of G23 coupling and is not affected by the choice of FEM.

If we consider a coarser mesh for high-order FEM in hopes of reducing the number of degrees of freedom, i.e., choosing , then applying analogous calculations to those in §6.7 gives us the following result:

Employing P-FEM with , in order to match the convergence rate of the Cauchy–Born error term , the highest mesh coarsening rate is

This means that the optimal mesh grading that P-FEM can achieve without compromising accuracy is no greater than . However, in that case, the number of degrees of freedom is always .

In summary, despite the possibility of (slightly) reducing the number of degrees of freedom, considering its computational complexity, we conclude that higher-order FEM is not worthwhile to pursue. However, we emphasize that this conclusion would need to be revisited if a coupling method with higher-order interface error as well as continuum model error could be devised.

4. Conclusion

We obtained a sharp energy-norm error estimate for the G23 coupling method with P2-FEM discretisation of the continuum model. Furthermore, we demonstrated that, with P1-FEM discretisation the FEM coarsening error is the dominating term in the consistency error estimate, whereas for P2-FEM discretisation the interface error becomes the dominating term. In particular, a P2-FEM discretisation yields a more rapid decay of the error. Crucially though, since for P-FEM with the interface contribution dominates the total error the P2-FEM is already optimal. That is, increasing to will not improve the rate of convergence, but increase the computational cost and algorithmic complexity.

Numerically, we observe that the improvement of P2-GR23 over P1-GR23 is only modest at low , hence a P2-GR23 scheme would be primarily of interest if high accuracy of the solution is required. However, our numerical results indicate that there is scope for further optimisation, using a posteriori type techniques.

While our estimates for the error in energy-norm are sharp, our numerical results show the estimates for the energy errors are suboptimal. We hightlight the leading term in the error analysis which overestimate the error in Section 7.2. We are unable, at present, to obtain an optimal energy error estimate. This appears to be an open problem throughout the literature on hybrid atomistic multi-scale schemes; see e.g. [4].

In summary we conclude that using P2-FEM is a promising improvement to the efficiency of a/c coupling methods, but that some further work, both theoretical and for its implementation may be need to exploit its full potential.

5. Reduction to consistency

Assuming the existence of an atomistic solution , we seek to prove the existence of satisfying


and to estimate the in a suitable norm.

The error analysis consists of consistency and stability estimates. Once these are established we apply the following theorem to obtain the existence of a solution and the error estimate. The proof of this theorem is standard and can be found in various references, e.g. [13, Lemma 2.2].


[The inverse function theorem] Let be a subspace of , equipped with , and let with Lipschitz-continuous derivative :

where denotes the -operator norm.

Let satisfy


such that satisfy the relation

Then there exists a (locally unique) such that ,

To ensure Dirichlet boundary conditions, we adapt the approximation map defined in [4]. Let be a cut-off function such that

For , define


Let be the vertices of and be the mid-point of an edge . Then, the set of all active P2 finite element nodes is given by

This includes all P1 nodes as well as the P2 nodes (edge midpoints) associatd with edges entirely in the P2 region.

Furthermore, let be the interpolation operator such that, for , for , , and


We also introduce ghost nodes on the edges shared by interface and continuum elements:


Then, for , , where and are the vertices of the edge on which lies. Hence, the and interpolants coincide on .

We can now define the projection map (quasi-best approximation operator) as


5.1. Stability

To put Theorem 5 (Inverse Function Theorem) into our context, let

To make (5.2) and (5.3) concrete we will show that there exist such that, for all ,

Ignoring some technical requirements, the inverse function theorem implies that, if is sufficiently small, then there exists such that

Finally adding the best approximation error gives the error estimate

The Lipschitz and consistency estimates require assumptions on the boundedness of partial derivatives of . For , define the first and second partial derivatives, for , by

and similarly for the third derivatives . We assume that the second and third derivatives are bounded


With the above bounds it is easy to show that


From the bounds above we can obtain the following Lipschitz continuity and stability results.


There exists such that


The result follows directly from the global bounds of derivatives of and the fact that and that is compactly supported hence is also Lipschitz. ∎


Under the assumptions (A1) and (A2), if , then there exits such that, when is sufficiently large,


The proof of this result is a straightforward adaption of the proof of [5, Lemma 4.9], which is an analogous result for blending-type a/c coupling. ∎

6. Consistency estimate with a P2-FEM

6.1. Outline of the consistency estimate

We begin by decomposing the consistency error into


where is given and we can choose arbitrarily.

For , for . But the test function in is a piecewise linear lattice function. While we postpone the construction of , we will ensure that it is defined in such a way that for all , where is an extra layer of atomistic sites outside . With this assumption in place, we can further decompose into the following parts,


where is the smooth interpolant of defined in Lemma 6.1 below. By in we mean the gradient of the canonical linear interpolant of . To estimate we require an approximation error estimate for . To estimate we will exploit the fact that the atomistic triangulation is uniform to prove a super-convergence estimate. Finally, for the modelling error, , we employ the techniques developed in [12].

To define the smooth interpolant , we use the construction from [5], namely a -conforming multi-quintic interpolant. Although the interpolant defined in [5] is for lattice functions on , we can use the linear transformation from to to obtain a modified interpolant.


(a) For each , there exists a unique such that, for all ,