Iterative Methods for the Force-based QC Approximation

Iterative Methods for the Force-based Quasicontinuum Approximation: Analysis of a 1D Model Problem

M. Dobson M. Luskin  and  C. Ortner
July 5, 2019

Force-based atomistic-continuum hybrid methods are the only known pointwise consistent methods for coupling a general atomistic model to a finite element continuum model. For this reason, and due to their algorithmic simplicity, force-based coupling methods have become a popular class of atomistic-continuum hybrid models as well as other types of multiphysics models. However, the recently discovered unusual stability properties of the linearized force-based quasicontinuum (QCF) approximation, especially its indefiniteness, present a challenge to the development of efficient and reliable iterative methods.

We present analytic and computational results for the generalized minimal residual (GMRES) solution of the linearized QCF equilibrium equations. We show that the GMRES method accurately reproduces the stability of the force-based approximation and conclude that an appropriately preconditioned GMRES method results in a reliable and efficient solution method.

Key words and phrases:
atomistic-to-continuum coupling, quasicontinuum method, iterative methods, stability
2000 Mathematics Subject Classification:
M. Dobson: CERMICS - ENPC, 6 et 8 avenue Blaise Pascal, Cité Descartes - Champs sur Marne, 77455 Marne la Vallée Cedex 2, France,
M. Luskin (Corresponding Author): School of Mathematics, 206 Church St. SE, University of Minnesota, Minneapolis, MN 55455, USA,
Christoph Ortner: Mathematical Institute, St. Giles’ 24–29, Oxford OX1 3LB, UK,
This work was supported in part by DMS-0757355, DMS-0811039, the Department of Energy under Award Numbers DE-FG02-05ER25706 and DE-SC0002085, the University of Minnesota Supercomputing Institute, the University of Minnesota Doctoral Dissertation Fellowship, the NSF Mathematical Sciences Postdoctoral Research Fellowship, and the EPSRC critical mass programme “New Frontier in the Mathematics of Solids.”

1. Introduction

The motivation for coupled atomistic/continuum models of solids is that the accuracy of an atomistic model is often only needed in localized regions of the computational domain, but a coarse-grained continuum model is necessary for the simulation of large enough systems to include long-range effects  [19, 16, 2, 5, 15, 25, 3, 24, 18, 28]. The force-based approach has become very popular because it provides a particularly simple and accurate [13] method for coupling two physics models without the development of an accurate hybrid coupling energy. It operates by creating disjoint subdomains in which the equilibrium equations at each degree of freedom are obtained by assigning forces directly from one of the physics models. In addition to coupling atomistic and continuum models, such an approach has also been found to be attractive, for example, in the coupling of regions modeled by quantum mechanics to regions modeled by molecular mechanics, since accurate hybrid coupling energies require an interfacial region that is too computationally demanding for the quantum mechanics model [4].

The force-based quasicontinuum (QCF) approximation is attractive because of its simple and efficient implementation and because it is the only known pointwise consistent quasicontinuum (QC) approximation for coupling a general atomistic model with a Cauchy-Born continuum model [13]. By consistent we mean that the absence of ghost forces under homogeneous deformations. Its main drawback is that it results in a non-conservative force field [6], that is, the QCF forces are not compatible with any energy functional. Several creative attempts have been made to develop hybrid coupling energies that satisfy the patch test (there are no resultant forces under uniform strain)  [29, 14], which is a weaker compatibility condition than pointwise consistency and leads to reduced accuracy.

In this paper, we consider the force-based quasicontinuum approximation (QCF),


but, for simplicity, we will focus mainly on its linearization about a reference state,

see Section 2 for the precise definitions. Recent analyses of the linearized QCF operator [13, 12] have identified both further advantages as well as disadvantages of the force-based coupling approach. In addition to being non-symmetric, which is related to the fact that is non-conservative, the linearized QCF operator also suffers from a lack of positive-definiteness [13]. In the present paper, we show that this somewhat unusual stability property of the operator presents a challenge for the development of efficient and stable iterative solution methods that is overcome by the GMRES methods we propose.

1.1. Framework for iterative solution methods

We consider three related approaches to the development of iterative methods for the QCF equilibrium equations (1). A popular approach [21] to solve the force-based equations (1) modifies a nonlinear conjugate gradient algorithm by replacing the univariate optimization of an energy, used for step size selection [23], with the computation of a step size such that the residual is (approximately) orthogonal to the current search direction. We will show in Section 4.2 that, due to the indefiniteness of , this method is not numerically stable for our QCF model problem.

The second approach we consider is the nonlinear splitting

to construct the nonlinear iteration equation


The iterative solution of the nonlinear splitting method (2) can then be obtained from the minimization of the sum of and the potential energy of the dead load where

that is,

For this approach to be accurate under conditions near the formation or motion of defects, care must be taken to ensure that the energy accurately reproduces the stability of the approximated atomistic system. We will see in Section 4.1 that using the original quasicontinuum energy defined in (18), which results in the ghost force correction (GFC) scheme, does not reliably reproduce the stability of the atomistic system [10] and can give a reduced critical strain for a lattice instability.

To develop the final approach, we recall the Newton method


where is the residual

The GMRES methods proposed and analyzed in this paper apply to the solution of the linear Newton equations (3) or their approximations. Since the QCF equilibrium equations are generally solved along a quasi-static process [7], a good initial guess is usually available and a small number of iterations of the outer iteration (3) is sufficient to maintain stability and accuracy.

1.2. Outline

We begin in Sections 2 and 3 by introducing the most important quasicontinuum approximations and outlining their stability properties, which are mostly straightforward generalizations of results from [13, 12, 10]. We also present careful numerical studies of the spectral properties of which are particularly useful for the analysis of Krylov subspace methods in Section 5.

In Section 4, we revisit the ghost force correction (GFC) scheme [27] which, as was pointed out in [6], can be understood as a linear stationary iterative method (2) for solving the QCF equilibrium equations. We show that, even though the QCF method itself is stable up to a critical strain , the GFC scheme becomes unstable at a significantly reduced strain for our model problem. This leads us to conclude (though the simple examples we analyze here can only be first indicators) that the GFC method is not universally reliable near instabilities. We note, however, that the GFC method can be expected on the basis of both theoretical [10] and computational results [10, 21] to be more accurate near instabilities than the use of the uncorrected QCE energy as explained in Section 4.1. Numerical results have also shown that the GFC method can give an accurate approximation of critical loads if the atomistic-to-continuum interface is sufficiently far from the defect [21, Figure 16], at a cost of a larger atomistic region than likely required by the accuracy of the QCF approximation.

The quasi-nonlocal energy of [29] given by (20) is a more reliable and accurate energy to use in the splitting iteration (2). It has been shown to reproduce the atomistic stability of one-dimensional atomistic systems with next-nearest neighbor interactions [10], and the error for multi-dimensional atomistic systems is likely to be acceptable if the longer-range interactions decay sufficiently fast. The splitting iteration (2) can then be used as part of a continuation algorithm for a quasi-static process [7] that provides the reliable detection of the stability of the atomistic system [10] as well as the improved accuracy for the deformation given by the force-based approximation [13].

We conclude Section 4 by proving the numerical instability of the modified conjugate algorithm [21] for our QCF model problem. We present these two examples to demonstrate the subtleties in designing an iterative algorithm for the solution of the QCF system and to underscore the need for thorough numerical analysis in the development of stable and efficient iterative methods for the QCF system.

We conclude by considering in Section 5 the generalized minimal residual method (GMRES) for the solution of the indefinite and non-symmetric QCF system. We provide an analysis of basic as well as preconditioned GMRES methods. We find in this section that a non-standard preconditioned GMRES method, based on the discrete -inner product, appears to have excellent stability properties up to the critical strain and a more reliable termination criterion.

2. Quasicontinuum Approximations and Their Stability

In this section, we give a condensed description of the prototype QC approximations and their stability properties. We refer the reader to [12, 10] for more details. Many details of this section can be skipped on a first reading and only referred back to when required.

2.1. Notation

Before we introduce the atomistic model and its QC approximations, we define the notation that will be used throughout the paper.

We consider a one-dimensional atomistic chain whose atoms have the reference positions for We will constrain the displacement of boundary atoms which gives rise to the displacement space

We will equip the space with various norms which are discrete variants of the usual Sobolev norms that arise naturally in the analysis of elliptic PDEs. For displacements and we define the norms,

and we let denote the space equipped with the norm. The inner product associated with the norm is

In fact, we use and to denote the -norm and -inner product for arbitrary vectors which need not belong to . In particular, we further define the norm


where , , and we let denote the space equipped with the norm. Similarly, we define the space and its associated norm, based on the centered second difference for (We remark that, for , we have that and .)

For a linear mapping where are vector spaces equipped with the norms we denote the operator norm of

If , then we use the more concise notation

If is invertible, then we can define the condition number by

When is symmetric and positive definite, we have that

where the eigenvalues of are

If a linear mapping is symmetric and positive definite, then we can define the -inner product and -norm by

We define the discrete Laplacian by


A definition of the inner product and norm that is equivalent to (4) can now be given by


Since is symmetric and positive definite, we can also define the inner product and “negative” norm by


2.2. The atomistic model

We consider a one-dimensional atomistic chain whose atoms have the reference positions for and interact only with their nearest and next-nearest neighbors. (For an explanation why we require instead of atoms as previously stated, we note that the atoms with indices will later be removed from the model, and refer to Remark 2.2.3 for further details.) We denote the deformed positions by , and we constrain the boundary atoms and their next-nearest neighbors to match the uniformly deformed state, where is a macroscopic strain, that is,


The total energy of a deformation is given by



Here, is a scaled two-body interatomic potential (for example, the normalized Lennard-Jones potential, ), and , are external forces. We do not apply a force at the atoms , which will later be removed from the model. The equilibrium equations are given by the force balance conditions at the unconstrained atoms,


where the atomistic force (per lattice spacing ) is given by


2.2.1. Linearization of .

To linearize (11) we let , , be a displacement from the homogeneous state that is, we define

We then linearize the atomistic equilibrium equations (10) about the homogeneous state and obtain a linear system for the displacement ,

where is given by

Here and throughout we define

where is the interatomic potential in (9). We will always assume that and which holds for typical pair potentials such as the Lennard-Jones potential under physically realistic strains . For example, if is the Lennard–Jones potential, and if then . This shows that the force to compress a chain to achieve a strain for which is several orders of magnitude larger than the force to fracture the chain.

2.2.2. Stability of .

The stability properties of can be best understood by using a representation derived in [10],


where is the continuum elastic modulus

Following the argument in [10, Prop. 1], we prove the following equality in  [11] which describes the stability of the uniformly stretched chain.

Proposition 1. If , then

where for some universal constant .

2.2.3. The critical strain

The previous result shows, in particular, that is positive definite, uniformly as , if and only if . For realistic interaction potentials, is positive definite in a ground state . For simplicity, we assume that , and we ask how far the system can be “stretched” by applying increasing macroscopic strains until it loses its stability. In the limit as , this happens at the critical strain which solves the equation


Remark 1. We introduced the two additional atoms with indices so that the uniform deformation is an equilibrium of the atomistic model. As a matter of fact, our choice of boundary condition here is very close in spirit to the idea of “artificial boundary conditions” (see [13, Section 2.1] or [17]), which are normally used to approximate the effect of a far field. In the quasicontinuum approximations that we present next, these additional boundary atoms are not required. ∎

2.3. The Local QC approximation (QCL)

The local quasicontinuum (QCL) approximation uses the Cauchy-Born approximation to approximate the nonlocal atomistic model by a local continuum model [6, 20, 24]. In our context, the Cauchy-Born approximation reads

and results in the QCL energy, for satisfying the boundary conditions (8),


Imposing the artificial boundary conditions of zero displacement from the uniformly deformed state, we obtain the QCL equilibrium equations



In particular, we see from (15) that the QCL equilibrium equations are well-defined with only a single constraint at each boundary (see also Remark 2.2.3), and we can restrict our consideration to with the boundary conditions and .

Linearizing the QCL equilibrium equations (15) about the uniformly deformed state results in the system

where , for a displacement is given by

The increased efficiency of the local QC approximation is obtained when its equilibrium equations (15) are coarsened by reducing the degrees of freedom, using piecewise linear interpolation between a subset of the atoms [6, 20]. For the sake of simplicity of exposition, we do not treat coarsening in this paper.

We note that

where is the discrete Laplacian (5). Since the QCL operator is simply a scaled discrete Laplace operator, its stability analysis is straightforward:

In particular, it follows that is stable if and only if , that is, if and only if where is the critical strain defined in (13).

2.4. The force-based QC approximation (QCF)

In order to combine the accuracy of the atomistic model with the efficiency of the QCL approximation, the force-based quasicontinuum (QCF) method decomposes the computational reference lattice into an atomistic region and a continuum region , and assigns forces to atoms according to the region they are located in. Since the local QC energy (14) approximates in (9) by it is clear that the atomistic model should be retained wherever the strains are varying rapidly. The QCF operator is given by [6, 7]


and the QCF equilibrium equations by

We recall that is a non-conservative force field and cannot be derived from an energy [6].

For simplicity, we specify the atomistic and continuum regions as follows. We fix , , and define

Linearization of  (16) about reads


where the linearized force-based operator is given explicitly by

We note that, since atoms near the artificial boundary belong to , only one boundary condition is required at each end.

We know from [12] that the stability analysis of the QCF operator is highly non-trivial. We will therefore treat it separately and postpone it to Section 3.

2.5. The original energy-based QC approximation (QCE)

In the original energy-based quasicontinuum (QCE) method [24], an energy functional is defined by assigning atomistic energy contributions in the atomistic region and continuum energy contributions in the continuum region. In the context of our model problem, it can be written as



The QCE method does not satisfy the patch test [22, 9, 8, 27], which be seen from the existence of “ghost forces” at the interface, that is, . Consequently, the linearization of the QCE equilibrium equations about takes the form (see [9, Section 2.4] and [8, Section 2.4] for more detail)


where, for we have

and where the vector of “ghost forces,” , is defined by

For space reasons, we only list the entries for The equations for follow from symmetry.

We prove in [11] the following new sharp stability estimate for the QCE operator which implies that the operator gives an O() approximation for the critical strain,

Lemma 2. If , , and , then

where . Asymptotically, as , we have

This result will be used in Section 4.1 where we analyze the ghost-force correction iteration, interpreted as a linear stationary iterative method for with preconditioner .

2.6. The quasi-nonlocal QC approximation (QNL)

The QCF method is the simplest idea to circumvent the patch test failure of the QCE method. An alternative approach was suggested in [29, 14], which is based on a modification of the energy at the interface. In this model, a next-nearest neighbor interaction term is left unchanged if at least one of the atoms belong to the atomistic region or an interface region (which is implicitly defined by (20)), and is otherwise replaced, preserving symmetry, by a Cauchy-Born approximation,

This idea leads to the energy functional


where we set . The QNL approximation satisfies the patch test; that is, is an equilibrium of the QNL energy functional.

The linearization of the QNL equilibrium equations about the uniform deformation is



We observe from (21) that is not pointwise consistent at and

Repeating our stability analysis for the periodic QNL operator in [10, Sec. 3.3] verbatim, we obtain the following result.

Proposition 3. If , and , then

Remark 2. Since the linearized operators and depend only on , and . ∎

3. Stability and Spectrum of the QCF operator

In this section, we collect various properties of the linearized QCF operator, which are, for the most part, variants of our results in [13, 12]. We begin by stating a result for the lack of positive-definiteness of which lies at the heart of many of the difficulties one encounters in analyzing the QCF method.

Theorem 4 (Lack of Positive-Definiteness of QCF, Theorem 1, [13]). If and then, for sufficiently large the operator is not positive-definite. More precisely, there exist and such that, for all and ,

As a consequence of Theorem 3, we analyzed the stability of in alternative norms. Following the proof of [12, Theorem 3] verbatim (see also [12, Remark 3]) gives the following sharp stability result.

Proposition 5. If and , then is invertible with

If then is singular.

This result shows that is operator stable up to the critical strain at which the atomistic model loses its stability as well (cf. Section 2.2). In the remainder of this section, we will investigate, in numerical experiments, the spectral properties of the operator for strains such that and .

3.1. Spectral properties of in

The spectral properties of the operator are crucial for analyzing the performance of iterative methods in Hilbert spaces. The basis of our analysis of in the Hilbert space is the remarkable observation that, even though is non-normal, it is nevertheless diagonalizable and its spectrum is identical to that of . We first observed this in [12, Section 4.4] for the case of periodic boundary conditions. Repeating the same numerical experiments for Dirichlet boundary conditions, we obtain similar results. Table 1, where we display the error between the spectrum of and gives rise to the following conjecture.

Table 1. The difference between the spectra of and The table displays the norm of errors in the ordered vectors of eigenvalues for various choices of with , for increasing , . All entries are zero to the precision of the eigenvalue solver.

Conjecture 6. For all the operator is diagonalizable and its spectrum is identical to the spectrum of .

We denote the eigenvalues of (and ) by

The following lemma provides a lower bound for an upper bound for and consequently an upper bound for . Assuming the validity of Conjecture 3.1, this translates directly to a result on the spectrum of .

Lemma 7. If and , then


It follows from Proposition 2.6 and (22) that

since the infimum of the Rayleigh quotient is attained for where  [31, Exercise 13.9] and has the value


The estimate for the maximal eigenvalue follows similarly from

and the representation (21). ∎

For the analysis of iterative methods, particularly the GMRES method, we are also interested in the condition number of the basis of eigenvectors of as tends to infinity. Assuming the validity of Conjecture 3.1, we can write where is diagonal. In Figure 1, we plot the conditio