An atomistic-based boundary element method for the reduction of molecular statics models

Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

An atomistic-based boundary element method for the reduction of molecular statics models Xiantao Li 1 Department of Mathematics, Pennsylvania State University, University Park, PA 16802, USA

a r t i c l e

i n f o

Article history: Received 24 February 2011 Received in revised form 13 January 2012 Accepted 8 March 2012 Available online 24 March 2012 Keywords: Molecular statics Boundary element Crystalline solids

a b s t r a c t We propose a new reduced computational model, called atomistic-based boundary element model (ABEM), derived from a full atomistic model for a crystalline solid system. The procedure is based on a domain decomposition method, which allows the separation of the atoms near crystal defects from the surrounding region, in which the displacement of the atoms is smooth. A reduction method, which is similar to the boundary integral method for continuum models, is developed to eliminate the atoms in the surrounding region. The reduction procedure gives rise to a system of equations only involving the atoms at the remote boundary and at the boundaries of the atomistic regions containing local defects. In this paper, we will discuss the derivation of the model, the implementation, and further reduction methods. We also present applications to some test problems. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction Atomistic models, which take into account detailed atomic interactions, have been important tools in studying mechanical properties of crystalline defects. One major difficulty in implementing such models is the dimensionality: usually the model involves a huge number of atoms, which makes the computation intractable. However, typical observations (e.g. see [45,17,13]) are that only a small fraction of the atoms are directly involved in the formation and migration of local defects, and for the majority of the atoms, the displacement is quite smooth, indicating that the models are reducible. The main purpose of this paper is to develop a reduction procedure to remove the excessive atomic degrees of freedom. The reduction is achieved in two steps: 1. We adopt a domain decomposition method, dividing the entire domain into several subdomains containing local defects, and one surrounding region that contains the remaining atoms, as illustrated in Fig. 1. The full atomistic description will be retained within the subdomains around defects, and these regions will be referred to as atomistic regions. 2. We reduce the dimension of the surrounding region by removing the atoms in its interior. The removal procedure has been motivated by the boundary integral formulation for continuum models.

1

E-mail address: [email protected] Research supported by NSF Grant DMS-1016582.

0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.03.006

In the second step, we will use the lattice Green’s function [38], and the summation by parts technique to eliminate the atoms in the interior of the surrounding region. As a result, we obtain a system of equations coupling the boundaries, including the physical boundary and the boundaries of the atomistic regions. Hence the number of degrees of freedom is greatly reduced. Due to its similarity to the traditional boundary element method, the current method will be called the atomistic-based boundary element method (ABEM). The reduced equations, expressed in terms of the displacement of the atoms at the boundaries, will be referred to as boundary atom equations (BAE). The first advantage of ABEM is the reduction of the dimension. A distinct advantage of traditional boundary element methods is the ability to deal with exterior problems, where the problem is posed in an infinite domain. For the present problem, the exterior domain is the one away from the local defects but still within the remote physical boundary. Although it is not infinite, in most practical problems the physical boundary is of a much larger scale, and the size of the exterior is much greater than that of the atomistic regions. Therefore, the reduction is considerable. The second advantage is that it is meshless: no mesh or grid points are needed, which can otherwise be difficult to generate since the atomistic regions and the physical boundary are of quite different scales. This advantage is more pronounced when the defects migrate under the external loading and the atomistic regions have to be re-defined as the computation continues. It is worthwhile to point out another crucial difference between the present approach and the traditional boundary element approach. For the traditional formulation, the process consists of deriving the boundary integral equation, and discretizing the

2

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

boundary. The current formulation follows a different path: It starts with a discrete model, and then the reduced equation is derived using summation by parts. At this point, no discretization of the boundary is needed. Under the current formulation, it is easy to identify the material interactions of multiple scales, a characteristic typical in many material problems. The Green’s function Gjjn for the discrete model depends on the relative distance between the jth and nth atoms. It reduces to the Green’s function for the continuum elasto-statics model for a large distance. The convergence is typically observed in about 6–10 atomic spacings. The boundaries, where the boundary atom equations are applied, will include the remote physical boundary, and the boundaries of the atomistic regions. Therefore, the stiffness matrix in ABEM can be written in a block form and each block can be assembled from the lattice Green’s functions. At this level, we see that each off-diagonal block represents the coupling between one atomistic region and another atomistic region, or the coupling between the atomistic regions and the remote boundary. These coupling will be approximately the same as that in a continuum model since the lattice Green’s function in this regime is well approximated by the continuum Green’s function. As a result, the coupling is of a continuum nature. This connection will be drawn later in Section 3.4. Meanwhile, within each diagonal block, the entries will represent the coupling between atoms within the same atomistic region, and since the distance between those atoms is typically small, the Green’s function, hence the coupling, is different from that at the continuum level. We have employed the lattice Green’s functions to derive the boundary atom equations. Such Green’s functions have been used directly as boundary conditions (BC) for problems where the surrounding region is infinite [31,34,43]. The main idea is to use the Green’s function to adjust the position of the atoms at the boundary based on the displacement of the atoms inside the domain. A more efficient formulation has been proposed in [19]. For finite exterior regions, the multiscale BC obtained by Karpov et al. [14,26] is essentially a special case of the BAE presented here. However, their derivation of the multiscale BC relies on Fourier transform, and as a result, the formulation only applies to periodic domains. Our formulation, however, is much more general, and it can be applied to domains of arbitrary geometry. Another class of methods for such multiscale problems have been based on coupling the molecular statics model with a continuum elasticity model. For instance, Mullins and Dokainish [29] proposed such a hybrid method in which a finite element linear elasticity model is used to provide boundary conditions for the atomistic regions. The influence of atomistic regions on the finite element region is taken into account by calculating the forces on some embedded atoms inside the elements near the interface. Later, Kohlhoff et al. [18] proposed a similar coupling method, called the finite element method combined with atomistic modeling (FEAt). In FEAt, the mesh at the interface is gradually refined to the atomic scale so that the finite element nodes coincide with the atoms. In addition, transition regions are created, where the nodal values of the finite elements and the displacement of the atoms are made consistent. Several other examples of coupling methods can be found in [2,5,16,37,30,15]. See [27] for a review of these methods and some comparative studies. The main departure of the current approach from these existing methods lies in the fact that it does not introduce the continuum model to begin with. The BAE provides a seamless transmission condition between subdomains. This avoids the problem of having to interface atoms with finite elements, which has been a rather complicated modeling and numerical issue [10,22,21,23]. The rest of the paper is organized as follows. In Section 2, we introduce the molecular statics model, and the partition of the domain. We then discuss in Section 3 the reduction procedure in

the surrounding region. Section 4 describes two domain decomposition methods. Then we present numerical examples in Section 5. 2. The molecular statics model This section describes the model to be considered in the paper. We consider a sample of a solid material with N atoms. We start with a full atomistic model involving all the N atoms in the system. Let ui ; 1 6 i 6 N, be the displacement of atom i from the reference position, Ri . The current position is given by, ri ¼ Ri þ ui . We consider the molecular statics model that seeks the mechanical equilibrium of the system by minimizing the energy,

min Vðu1 ; u2 ; . . . ; uN Þ þ

P

bi  ui ;

ð1Þ

i

where V is a potential energy is a possible body force. Alternatively, one can solve the force balance equation,

f i ¼ bi ;

fi ¼ 

@ V: @ui

ð2Þ

It is important to point out that the reference state may not have to be the one defined by the perfect lattice structure. Indeed, when dislocations or cracks are present, the coordinate of the neighbors of an atom may be different from that in the underlying lattice due to broken bond or slip. In this case, it is more suitable to choose a current configuration as the reference coordinate [44]. For example the reference coordinate can be chosen based on an analytical solution of the linear elasticity model for dislocations or cracks [40]. It is also important to emphasize that these reference coordinates are chosen only to identify the correct neighbors for atoms in the surrounding region. Typically, an empirical potential can be written as,

V¼

N P

V i:

i¼1

Namely the energy can be defined at each atom. In addition, most potential models have a cut-off radius rcut , in that the jth atom exerts a force on the ith atom only when their distance is smaller than rcut . This implies that,

@2V ¼ 0; @ui @uj

if kri  rj k P r cut :

ð3Þ

The surrounding region Atomistic regions

Physical boundary Fig. 1. A schematic picture of the domain decomposition. The shaded area is the surrounding region.

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

3

the observation that in the region B, which is away from local defects, the displacement there is typically smooth, in the sense that kui  uj k  kRi  Rj k. Within this approximation, we solve the linear equation,



P

Dk ujþk ¼ bj ;

ð6Þ

k2N

for all j 2 B, while the atoms in D and @ Xþ are left alone. Here Dk is the force constant, which can be computed from the potential energy around the reference state. Namely,

Dij ¼

Fig. 2. The partition of the domain. The physical boundary, and the boundary of the atomistic regions are shown in solid lines. Each boundary is defined by a curve, and the nearby atoms outside and inside those boundaries are indicated in the figure. The reduction procedure will remove the atoms in the interior of the surrounding region (hidden) and yield an equation involving those atoms close to the boundaries (shown). Further, the sets @Bþ and @B are shown as open and filled circles, respectively. The atoms in the interior of the atomistic regions (shown as squares) do not affect the reduction procedure because of the short range interaction.

This allows us to define the atoms close to the boundary, which will be described below. We first partition the entire sample into subdomains. Let X be the entire domain in the reference coordinate. We also define subdomains, Xa , to be called atomistic regions, in which the atoms will be kept in the computation. As a result, the actual computational domain based on the molecular statics model, called D, consists of several subdomains of X, i.e. D ¼ [a Xa . In addition, we let B be the surrounding region, B may or may not intersect Xa depending on the domain decomposition method (see Section 4). The partition of the domain is depicted in Fig. 2. Our reduction procedure will remove the atoms in B and replace that part of the molecular statics model by a more efficient one. Throughout this paper, we will use the notation i 2 S for any set S to indicate that the reference position of the ith atom is in the set S, i.e. Ri 2 S. We define the outer and inner boundary as follows,

@ Xþ ¼ fRi R X; jRi  Rj j 6 r cut ; for some j 2 Xg;

ð4Þ

and

@ X ¼ fRi 2 X; jRi  Rj j 6 rcut ; for some j R Xg:

ð5Þ

@2V : @ui @uj

ð7Þ

From (3), we see that the force constant matrices are only nonzero when the two atoms are within the cut-off radius. Since in the region B the crystal structure is intact, and the atomic interaction is local, we will use the force constant Dk computed around a perfect lattice structure. The translational invariance of the force constant matrices come from the translational symmetry of the potential V. We assume a single crystal structure so that Dij is uniform throughout the region B. In Eq. (6), bj is the external force on atom j, similar to the body force in continuum mechanics. We have used the convention that the index j þ k is referring to the atom with reference position at Rj þ Rk . Therefore, if the zeroth atom is the atom with reference position at the origin, then N is a set containing the neighbors of the zeroth atom that are within the cut-off distance. With the above notation, we can write j þ N to include the neighbors of the jth atom with reference position Rj . These neighbors may be outside B. In the second step of the reduction procedure, we introduce the lattice Green’s function, Gjjn , the solution of the lattice statics model for the infinite lattice with an applied point-force at the point n. Namely,

P

Dk Gjþkjn ¼ Idjn :

ð8Þ

k2N

Here the index j indicates the point where the function is being evaluated, and n is the center of the Green’s function. The existence of the lattice Green’s functions has been studied in [24]. For computational purposes, we will use the Fourier integral representation (A.1) provided in the appendix. In this paper, we consider simple lattices, for which Dij satisfies the following properties [1],

Dk ¼ DTk ;

ð9aÞ

Dk ¼ Dk ; P Dk ¼ 0;

ð9bÞ ð9cÞ

Similarly, we can define @ X a . This is also shown in Fig. 2.

k2N

3. Reduction of the linear molecular statics model

As a result of (9a), Gjjn represented by the Fourier integral (A.1) is symmetric, i.e., Gjjn ¼ GTjjn . Using the property (9c), Eq. (6) can be written as,

Dij ¼ 0;

3.1. Derivation In this section, we provide the details of the reduction procedure for the molecular statics model in B. The molecular statics model in the atomistic regions D will be treated later in the next section. We will assume that the lattice points in B can be mapped to a subset of the infinite Bravais lattice. Secondly, we introduce the harmonic approximation [1] in this domain. This is based on

Table 1 Comparison of the number of iterations. Overlapping size (4.032 Å)

34

46

68

Number of iterations

57

39

26



P

if jRi  Rj j > r cut :

  Dk ujþk  uj ¼ bj :

ð9dÞ

ð10Þ

k2N

Now for our domain B, which is a connected (possibly multiconnected) domain, let @Bþ and @B be respectively the outer and inner parts of the boundary, as discussed in the previous section, and set @B ¼ @Bþ [ @B . For instance, for the domain  shown in Fig. 2, we would have @Bþ ¼ @ Xþ [ X 1 [ X2 , and  þ þ  @B ¼ @ X [ X1 [ X2 . We first formulate a boundary element model to compute the solution of the linear molecular statics model in B subject to certain boundary conditions along @Bþ . Remember that the purpose of our reduction procedure is to reduce the number of degrees of freedom associated with the atoms in B. We begin the derivation by showing a summation by

4

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

parts formula. Without loss of generality, we consider two scalar discrete functions, fj and g j , defined on B and its boundary. Then the following formula can be directly verified,

P



 P  fjþk  fj g j ¼  fj g j  g jk þ

j2B

j2B

þ

P

P

djn uj ¼

j2B

PP

j2B;jþkRB

fj g jk :

ð11Þ

Gjkjn Dk uj :

ð12Þ

j2Bk2N

In the second step, we took the transpose of Eq. (8) and used the fact that the Green’s function is symmetric. Using the symmetry of the force constant (9b) and (9c), we can write this equation as,

PP ðGjþkjn  2Gjjn þ Gjkjn ÞDk uj ;

2un ¼

k2N j2B

P P

¼

   Gjþkjn  Gjjn  Gjjn  Gjkjn Dk uj :

k2N j2B

We now proceed with our derivation using the summation by parts formula (11). We have,

PP 2un ¼  ðGjþkjn  Gjjn ÞDk ðujþk  uj Þ k j2B



P P

ðGjjn  Gjkjn ÞDk uj

k j2B;jkRB

þ

P P

ðGjþkjn  Gjjn ÞDk ujþk ;

k j2B;jþkRB

¼

P P

Gjjn Dk ðuj  ujk Þ  ðGjjn  Gjkjn ÞDk uj

k j2B;jkRB

þ

P P

j2B

At the last step, we have substituted (10) into the equation above. Collecting terms, we have,

P i2@B ;j2@B

þ

i2@Bþ ;j2@B

ðGijn  Gjjn ÞDij ðui þ uj Þ  2

The modeling error is given by the difference between the solution of the full molecular statics model (2) and the linearized model (6). It is clear that the error, denoted by ei , satisfies the equation,

P

Dk eiþk ¼ gi ðuÞ;

ð15Þ

k2N

where u on the right hand side stands for the exact solution, the solution of the full molecular statics model (2). This is a system of finite difference equations. There are two necessary steps in estimating ei : (1) the estimate of gi ; (2) the stability of the difference operator on the left hand side. We have analyzed the anharmonic part of the force for pair potential and the EAM potential models. With proper scaling, it can be shown that the linearization error is bounded by e times divided differences of the displacement up to third order, a discrete analog of derivatives. Here e is the lattice spacing. Meanwhile, for the stability of the difference operator, one can introduce a discrete coercivity condition, which is similar to the analysis of elliptic partial differential equations with variable coefficients. Such coercivity condition is similar to the phonon stability condition [11]. This analysis can also be used to analyze approximations of the boundary atom equations, including the approximation of remote boundary conditions, and the approximations of the lattice Green’s functions, which will be discussed in later sections. However, the details of such analysis is outside the scope of this paper, and it will appear in a separate paper. In this paper, we will mainly focus our attention on the formulation of the methods and practical implementations.

P

un ¼

j2@B



J n;j uj 

P

In;i ui ;

ð16Þ

i2@Bþ

where,

ðGijn þ Gjjn ÞDij ðui  uj Þ 

P

ð14Þ

We first assume that the displacement of the atoms in @Bþ is given. In this case, we write Eq. (13) as: For each n 2 @B ,

P  2 Gjjn bj :

þ

Dk uiþk þ gi ðuÞ:

3.2. Boundary conditions of Dirichlet type

 Gjþkjn Dk ðujþk  uj Þ þ ðGjþkjn  Gjjn ÞDk ujþk

k j2B;jþkRB

2un ¼ 

P

k2N

This is a discrete analog of the integration by parts formula. In particular, the last two terms on the right hand side can be viewed as boundary terms. To begin with, we notice that, for each n 2 B

P

fi ¼ 

fjþk g j

j2B;jkRB

un ¼

more precise, we can decompose the force on each atom into a harmonic part and an anharmonic part,

P

Gjjn bj :

ð13Þ

j2B

8 P Gijn Dij ; > < J n;j ¼ i2@Bþ P > Gjjn Dij : : In;i ¼

ð17Þ

j2@B

Here we have used the last property of the force constant (9d), and we only need to keep the atoms near the boundary in the summation. Eq. (13) holds for each n inside B. By restricting n to be at the inner boundary, i.e., n 2 @B , Eq. (13) is closed. As a result, the equation above, to be referred to as the boundary atom equation (BAE), only involves the atoms at the boundary. Therefore, the number of degrees of freedom has been greatly reduced. In order to solve (13), boundary conditions have to be taken into account. Next we discuss how to incorporate boundary conditions of different types. We will assume that no body force is present. For twodimensional square lattices with scalar displacement field and nearest neighbor interactions, equations of this type have been derived by Martinsson [24]. The current formulation is more general. It can be applied to both two and three dimensional models with interactions beyond nearest neighbors. We close this section with some comments on the modeling error. From the linearized molecular statics model (6) to the boundary atom equation (13), the reduction procedure is in principle exact. Hence the modeling error only comes from the linearization of the interatomic forces, i.e., the harmonic approximation. To be

By restricting n to be any atom at the inner boundary, we get a linear system for the displacement of the atoms at the inner boundary @B , given the displacement of the atoms at the outer boundary @Bþ . After Eq. (16) is solved, one can move n to any atom in B to evaluate the displacement of that atom. 3.3. Boundary conditions of Neumann type We now discuss how to impose traction boundary condition. For each atom at the inner boundary, the applied traction is defined as the net force due to the atoms from outside. More precisely, for each j 2 @B ,

tj ¼ 

P

Dij ðui  uj Þ:

ð18Þ

i2@Bþ

We also define the traction for the Green’s functions centered at n 2 @B ,

Sn;j ¼ 

 P Gijn  Gjjn Dij :

i2@Bþ

Now Eq. (13) can be rewritten as: For each n 2 @B ,

ð19Þ

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

un ¼ 

P j2@B

Sn;j uj þ

P j2@B

Gjjn tj :

5

ð20Þ

Assuming that the traction along the boundary is given, this equation may be solved to obtain the displacement at the inner boundary, which can then be used to compute the displacement of any atom inside B. The traction is typically computed from the previous step in a domain decomposition method, which will be discussed in Section 4. It is also possible to impose boundary conditions of mixed type. For this purpose we divide the boundary into two subsets,

@B ¼ @BD [ @BN ; where Dirichlet boundary condition and traction boundary condition are applied to the portion @BD and @BN of the boundary, respectively. In this case, the boundary atom equation is turned into

un ¼ 

P j2@B \@BN

Sn;j uj þ

P



P j2@B \@BN

Gjjn tj þ

P j2@B \@BD

Fig. 3. The traction on atoms at the boundary. The atoms at the inner boundary experience traction force from the atoms outside the boundary.

J n;j uj

In;i ui :

ð21Þ

i2@Bþ \@BD

Combining (24) and (23), the total traction per unit area at the boundary becomes,

3.4. Comparison with the integral equation of the continuum elastostatics model

si ¼

P l

The formulation of the boundary atom equations is quite similar to direct formulations of the boundary integral equations for elastostatics [7]. We will show that in the case when the displacement is smooth, further connections can be made. However, we do not attempt to prove the convergence of the BAE to a boundary integral equation, which requires more assumptions and technicalities. A similar problem can be found in [4] where an analysis on an approximation of a boundary integral using values on regular grid points is given. We recall the boundary integral equation (BIE) for the elastostatics model [7,8],

uðxÞ ¼ 

Z

Tðx; yÞuðyÞdSðyÞ þ

@B

Z

Gðx; yÞsðyÞdSðyÞ:

ð22Þ

@B

Here x 2 B and s is the traction at the boundary,

s ¼ P  n; Pil ¼

P j;k

Eikjl

@uj ; @X k

ð23Þ

where n is the outward normal and Eikjl stands for the elastic constant. G is the elastostatics Green’s function, and T is the traction defined by the Green’s function, and its expression can be found in [7]. The resemblance of the BAE (20) to (22) is apparent. The resemblance of the BAE (13) to the BIE can be noticed when the term 1 ðui þ uj Þ is interpreted as an average approximation of the dis2 placement at the boundary, and the term Dij ðui  uj Þ is viewed as the traction. For the continuum model, if a body force bðxÞ is present, a volume integral will appear in (22) in the form of R Gðx; yÞbðyÞdy. The last term in the BAE (13) is a discrete analog B of this volume integral. To make further connections, we first consider traction along a planar boundary. We recall that the elastic constants in the continuum model can be computed from the force constants of the molecular mechanics model [1],

Eikjl

1 P ¼ Dij ðRÞRk Rl : 2V 0 R

ð24Þ

Here we consider all the neighboring atoms around the origin. DðRÞ is the force constant associated with the point R. V 0 is the volume of the unit cell in the reference coordinate.

Pil nl ¼

P @uj Eikjl nl : @X k

ð25Þ

j;k;l

Therefore, in light of (24) we have,

s¼

1 P DðRÞðR  rÞuðXÞðR  nÞ: 2V 0 R

ð26Þ

Next we turn to the traction (18) for the boundary atom Eq. (13). We will assume that the atoms near the boundary can be arranged into boxes of the same size. This usually can be done when the boundary is a major symmetry plane. In addition, the lattice points can be divided into layers with uniform spacing a. Each box has a face, denoted by A, which overlaps with the boundary. The area is jAj. The volume of the unit cell in this case is given by V 0 ¼ jAja. This is shown in Fig. 3. We now proceed to compute the traction per unit area. Let C be a box at the inner boundary. Using (18) we have,

1 P 1 PP ti ¼  Dij ðuj  ui Þ: jAj i2C jAj i2C j2@Bþ

ð27Þ

We discuss the connection to the BIE when the displacement is smooth. The smoothness of a discrete function can be measured by discrete Sobolev norms, often used in the analysis of difference equations, e.g. [39,6]. An alternative is to assume that there exists a smooth function that agrees with the discrete function at the atom positions, and its gradient is small [11]. More specifically, we assume that there exists a smooth function u such that, uðRj Þ ¼ uj for each j. In addition, jruj  1. A Taylor expansion yields,

  uj  ui ¼ ðRj  Ri Þ  r uðXÞ þ oðjrujÞ; for a point X in A. The traction is approximated by,

  1 P 1 PP ti   Dij ðRj  Ri Þ  r uðXÞ: jAj i2C jAj i2C j We now move the double summation to summations over every pair of the layers. We notice that across the boundary, there is one pair of planes with distance a, two pairs with distance 2a, etc., yielding,

6

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

(

P 1 P 1 P ti   m DðRÞðR  rÞuðXÞ jAj i2C jAj mP0 Rn¼ma ¼

1 P DðRÞðR  rÞuðXÞðR  nÞ: 2jAja R

 @u@ V a ¼ 0; i

ui ¼ ð28Þ

uþi ;

i 2 Xa ;

i 2 @ Xþa ;

ð29Þ

We now return to the full problem on the entire domain X, and discuss how to combine the solution in B with the solutions from the subdomains around local defects. Let Xa be a subdomain in which the nonlinear molecular statics model is used, and inside each Xa , we solve the nonlinear molecular statics equation,

subject to a Dirichlet boundary condition at @ Xþ a . The displacement uþ j is determined from the BAE. The energyPV a is the potential energy restricted to the subdomain Xa ; V a ¼ i2Xa [@ Xþa V i . We will discuss both overlapping Schwarz and non-overlapping domain decomposition methods with Dirichlet–Neumann coupling. In both methods, the displacement of the atoms at the boundary of the atomistic regions are computed from the boundary atom equation, and in the subsequent step, they are held fixed while the rest of the atoms are relaxed. As a result, the net force on each atom in the interior of the atomistic regions is made zero. However, there may be remaining forces on the atoms at the boundary. Therefore, the boundary atom equation needs to be solved again to remove the remaining force, and the procedure will be iterated until some convergence criterion is met. Alternating Schwarz method. A natural approach is to extend Xa to create an overlapping region with B, and we denote the extended domain as Xa . This is illustrated in Fig. 4. The idea is the same as the overlapping Schwarz method in the context of domain decomposition methods. More specifically, within the region B, the linear molecular statics model is solved as described in the previous section. Boundary conditions, either of Dirichlet or Neumann type, are applied at the remote boundary. In addition, a Dirichlet boundary condition is applied to the atoms at the inner boundary of Xa . Solving the boundary atom equations (16) with these given boundary conditions, one obtains the displacement of the atoms at the inner boundary of B. Using the displacement of the atoms at the boundary of B, one can compute the displacement of any atom in B from Eq. (16). In particular, we will compute the displacement of the atoms at the outer boundary of Xa , which in turn provides the boundary condition for the nonlinear molecular statics model in Xa . These steps will be repeated until convergence. In order to quickly examine the convergence, we compare the solutions of the boundary atom equation at successive steps, and stop the iterations if the difference is below certain threshold. Non-overlapping method. For the non-overlapping method, no overlapping region is created to exchange boundary conditions. We first compute the traction at the outer boundaries of Xa . This,

Fig. 4. Overlapping domain decomposition: The domain is partitioned into subdomains that contain local defects, and one surrounding region, the boundary of which is indicated by solid lines. The shaded region is the overlapping region. The position of the atoms that are shown in the figure is fixed to impose Dirichlet boundary conditions.

Fig. 5. Non-overlapping domain decomposition: The domain is partitioned into subdomains that contain local defects, and one surrounding region, the boundary of which is indicated by solid lines. The problem is thus divided into two. In problem I, we solve BAE in the surrounding region. The atoms with arrows indicate the points where Neumann type of boundary condition is applied. For problem II, we solve the molecular statics model with boundary atoms fixed to position that is determined by the BAE in problem I.

In the last step, we have used the symmetry property of the force constant (9a). Finally, we notice that V 0 ¼ jAja. Therefore for a smooth displacement, the traction defined in the molecular statics model (27) agrees with that of the continuum model (26). Meanwhile, for a large distance, i.e. kRj  Rn k  1, the lattice Green’s function Gjjn agrees with the Green’s function GðRj ; Rn Þ, of the continuum equation. Therefore, the first terms of the right hand side of (22) and (13) also agree. At this point, we can observe the multiscale nature at the level of the boundary atom equations. Namely, the coupling between two atoms that are far apart is approximately the same as the coupling in a continuum boundary integral model. But when the distance is comparable to the atomic spacing, the coupling is of atomistic nature. The boundary atom equation will differ from the boundary integral equation in the following cases. 1. Near corners. In this case, the tractions defined in the atomistic and continuum models do not agree in general due to the lack of translational symmetry. 2. When x is close to y in (22), or when n is close to j in (13), especially when the distance is comparable to the atomic spacing. In this case, these two Green’s functions are quite different (see the Appendix and Fig. 12 for the comparison). In fact for the continuum model, the Green’s function becomes singular, and the integral has to be interpreted using Cauchy principal value. But for the atomistic model, the Green’s function remains finite. 4. Coupling method based on domain decomposition

7

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

together with the remote boundary condition, forms a complete set of boundary conditions for the boundary atom equation (20) at the boundary of B. Again the solution of (20) provides the displacement of the atoms at the outer boundary of Xa . The next step is the same as the overlapping method: we use the displacement of the atoms at the outer boundary of Xa , and solve the nonlinear molecular statics model in the subdomains. This coupling method, as demonstrated in Fig. 5, is similar to the Dirichlet–Neumann algorithm in domain decomposition method [41] (Section 1.3). For this coupling, it is often necessary to include a relaxation parameter, c. Let j 2 @ Xþ a be an atom at the outer boundary of a subdomain. Let unj be the displacement at the previous step, and unþ be the displacement computed from the BAE. Then we update j the displacement as follows:

unþ1 j

unþ j

¼c

Þunj :

þ ð1  c

4.1. Further reduction of the boundary atom equations The BAE involves all the atoms at the boundary of the subdomains and the remote boundary. In the case when the remote boundary is large, most boundary atoms are near the remote boundary, and the size of the BAE can be quite large. If the boundary condition at the remote boundary is smooth, further reductions can be achieved. We assume that the remote boundary conditions are of Dirichlet type. We first discuss the reduction in the alternating Schwarz method, which is supplemented with Dirichlet boundary condition from the subdomains. We express the displacement of the atoms at the boundary of B using a smaller set of nodal basis. This is done as follows. Let L ¼ dimð@BÞ. Then the displacement of the atoms is in the space of X ¼ RL . Let Y be a M-dimensional subspace of X, and let um ðRÞ; 1 6 m 6 M, be a set of M basis functions of Y. To make the formulation general, we do not specify the exact form of the basis functions. See the next section for an example. We look for an approximate solution of (13) in Y,

uk ¼

dm um ðRk Þ;

k 2 @B:

ð30Þ

m¼1

M P m¼1

m¼1

Mkm dm ¼

dm um ðRi Þ;

for all i 2 @Bþ :

P

¼

i2@B

ð32Þ

þ

uk ðRi Þum ðRi Þ;

ð33Þ

and

bk ¼

P

Mkm ¼ Kkm

n2@B

P i2@Bþ

uk ðRi Þui :

ð35Þ

uk ðRn Þum ðRn Þ;

P P

¼

n2@B j2@B

P P

¼

n2@B i2@Bþ

J nj um ðRj Þuk ðRn Þ;

ð36Þ

Ini um ðRi Þuk ðRn Þ:

Now let M ¼ M þ Mþ , and K ¼ K  Kþ , which will be referred to as mass and stiffness matrices. Combining Eqs. (35) and (32), we arrive at the linear system, M P

½Mkm  Kkm dm ¼ bk :

ð37Þ

m¼1

Now we turn to the non-overlapping method with Neumann condition applied at the boundaries of the subdomains Xa ; a ¼ 1; 2; . . . ; nd , with nd being the number of atomistic regions. In this case, B ¼ X n [a Xa . As a result, we have

@B ¼ @ X [ @ Xþ1 [    [ @ Xþnd ;

ð38Þ

@Bþ ¼ @ Xþ [ @ X1 [    [ @ Xnd :

The boundary element with mixed type of boundary conditions (20) can now be applied, yielding,

un ¼ 

nd P P

Sn;j uj þ

a¼1j2@ Xþa

nd P P

Gjjn tj þ

a¼1j2@ Xþa

P 

J n;j uj 

j2@ X

ð34Þ

P

In;j ui :

ð39Þ

i2@ Xþ

The displacement of the atoms at the remote boundary is represented by N 0 elements. More specifically, for each k 2 @ Xþ [ @ X , we have

uk ¼

N0 P m¼1

dm um ðRk Þ:

ð40Þ

We first define the load vector,

bk ¼

P

ui uk ðRi Þ;

ð41Þ

which in matrix notation, can be written as,

ð42Þ

where,

Mþkm ¼

P i2@ Xþ

uk ðRi Þum ðRi Þ:

ð43Þ

Next, we multiply Eq. (39) by uk ðRn Þ, sum over @ X , and then add it to (43). This leads to,

where the matrix Mþ is given by,

Mþkm

m¼1

 Kkm  Kþkm dm ;

Mþ d ¼ b; ð31Þ

Multiplying the equation by uk ðRi Þ, and summing over i, we get the following system,

Mþkm dm ¼ bk ;

M  P

in which, the matrices are given by,

i2@ Xþ

Here dm 2 R3 are the nodal values for the displacement. In particular, using the Dirichlet boundary condition, we have,

ui ¼

M P

Kþkm

Remark. For the nonlinear system of equations associated with the atoms in the atomistic region, it is natural to impose Dirichlet boundary conditions. Namely the atoms at the boundary are held at certain given positions. A Neumann boundary condition, however, is not easy to impose, especially for many-body potentials, such as the embedded atom potential [9]. Therefore, for the coupling methods, a Neumann–Neumann coupling idea remains open.

M P

This vector, which is determined by the boundary conditions, will be called a load vector. Next, we consider Eq. (16) for any n 2 @B . Inserting (30), multiplying the equation by uk ðRn Þ, and summing over all n 2 @B , we arrive at,

M0 d ¼

nd P

a¼1

K0;a ua þ K0;0 d þ

nd P

C0;a sa þ b:

ð44Þ

a¼1

Here we have written the system in block forms. ua is the displacement of the atoms at the boundary of the ath atomistic region, expressed in a column vector form. The block matrices are defined as follows:

8

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

P

ðM0 Þkm ¼ K0;0 k;m ¼

P 



n2@ X ;j2@ X

a K0; k;j ¼ 

P

n2@ X

a C0; k;j ¼

uk ðRi Þum ðRi Þ;

i2@ X

P n2@ X

J n;j um ðRj Þuk ðRn Þ 

Sn;j uk ðRn Þ;

Gjjn uk ðRn Þ;

P n2@ X ;i2@ Xþ

In;i um ðRi Þuk ðRn Þ;

j 2 @ Xþa ; j 2 @ Xþa : ð45Þ

Now we move the nth atom to the boundary of Xa . Inserting (40) into (39), we find,

ua ¼

nd P

Ka;b ub þ Ka;0 d þ

b¼1

nd P

Ca;b tb ;

ð46Þ

b¼1

where the block matrices are given by,

8 a;0 P P Kn;j ¼ J n;j um ðRj Þ  In;i um ðRi Þ; > >  > j2@ X i2@ Xþ <

n 2 @ Xþa ; ð47Þ

a;b ¼ Sn;j ; n 2 @ Xþa ; j 2 @ Xþb ; Kn;j > > > : a;b Cn;j ¼ Gjjn ; n 2 @ Xþa ; j 2 @ Xþb :

Here 1 6 a; b 6 nd and the superscripts indicate the atomistic regions or the remote boundary. The subscripts indicate the entries. Finally, we assemble the linear system in a block form,

2 6 6 4

3 2 0;0 d K K0;1 1;0 1 76 u 7 6 K K1;1 76 . 7 ¼ 6 .. 6 5 5 4 .. 4  .  I un d Knd ;0 Knd ;1 2 I C0;1 6 0 C1;1 6 þ6 4  0 Cnd ;1 32

M0 I

  .. .    .. . 

32 3 d K0;nd 1;nd 7 1 u 7 K 76 7 76 4 ... 5 5  un d Knd ;nd 3 2 3 b C0;nd t1 7 C1;nd 7 76 6 74 . 7 5:    5 .. nd nd ;nd t C

Fig. 6. Nodes and elements at the boundaries. At remote boundaries, we use rectangular elements. At the boundaries of the atomistic regions, each atom is a node. The subdomains are defined for a system with two dislocations (see the next section).

5.1. Dislocation dipole

ð48Þ 5. Numerical examples Here we show two numerical examples. In both examples, the molecular statics model is three-dimensional, with periodic boundary conditions in the z direction. In the surrounding region, we reduce the linear molecular statics model to the x  y plane by projecting the force constant matrices. Therefore, the BAE is effectively for a two-dimensional molecular statics model. During the iterations of the domain decomposition methods, we use the average of an atom and its images in the z direction to compute the boundary conditions for the BAE. On the other hand, the boundary condition for the molecular statics model needs to be computed from BAE, and it is computed based on the ðx; yÞ coordinates of an atom. The Green’s functions are computed using Fourier integral approach, and for a large distance, it is replaced by the continuum elasto-statics Green’s function. The calculation of the elasto-statics Green’s function is explained in the Appendix. For the solution of the BAE, we choose rectangular elements around the remote boundary. Each element along the edges is 7 atomic spacings long, and 1 atomic spacing in width. The basis function is bilinear. For the atoms around the atomistic regions, we choose the trivial basis function, i.e., um ðRn Þ ¼ dmn . The nodal points are shown in Fig. 6 for the first test problem. For the molecular statics model, the BFGS optimization method [20] is used.

Our first example is a dislocation dipole in aluminum. The atomic potential used here is the EAM potential [12]. The system studied consists of a 3D rectangular sample, with the three orthog directions respeconal axes being along the h110i; h001i and h110i tively. The system is periodic in z-direction with period equal to 10 times the atomic spacing. The dimension of the entire system is 84:3 nm  84:3 nm, and the full system contains 876,160 atoms. This example has been used as a benchmark problem to compare a number of coupling methods [27]. But in our test, the entire sample is much bigger, and we used two separate atomistic regions for the dislocation cores instead of one. A shear strain is applied at the remote boundary. The solution obtained from the ABEM method is shown in Fig. 7 for e12 ¼ 0:20, where the first component of the displacement is plotted. In Fig. 8, the error is plotted, and the error is within 3% throughout the region. For each loading step, the full atomistic simulation takes about 135 min, and the ABEM with the alternating Schwarz method takes 40 min. For this example, the size of each atomistic region is 11:40 nm  9:68 nm, and each region contains 20,160 atoms. Together, they are about 1/24 the size (in area) of the entire sample. In addition to the atomic degrees of freedom in these two regions, there are 876 elements with 1052 nodes along the interface and the remote boundary (see Fig. 6). We monitored the CPU time within a domain decomposition step. The solution procedure for solve the BAE (13) takes 0.006 (seconds), and the CPU time for the molecular statics model in the atomistic regions is 53.184 (seconds). We see that with the reduction using ABEM, most of the CPU time is spent on computing the solution in the atomistic regions. Therefore the reduction of the computational cost will be much more dramatic if the entire sample is larger. We also compare the number of iterations in the domain decomposition method for different sizes of overlap. The results are shown in Table 1.

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

0.8

250 200

0.6

150

0.4

100

0.2

50 0

0 −50

−0.2

−100

−0.4

−150

−0.6

−200 −0.8

−250 −400

−300

−200

−100

0

100

200

300

400

−1

1 250

0.8

200

0.6

150 0.4

100

0.2

50 0

0

−50

−0.2 −0.4

−150 −0.6

−200

−0.8

−250 −400

−300

−200

−100

0

100

200

300

400

−1

Fig. 7. Dislocation dipole. Top: Solutions obtained from the ABEM method; bottom: solution obtained from the full atomistic simulation.

0.04 250 0.03

200 150

0.02

100 0.01

50 0

distance 8.107 Å, which includes 148 neighbors. Notice that this cut-off distance is larger than that of the electron density, which is 4.095 Å. Next we consider an elliptical crack in the BCC Iron-a system. The system studied is a 3D rectangular sample, with the three  1 0i and h0 0 1i directions orthogonal axes along the h1 1 0i; h1 respectively. The dimension of the entire system is 88:6 nm  88:6 nm  2:13 nm, consisting of 1,520,768 atoms. The crack is chosen to lie on the {1 1 0} plane, and the crack front is parallel to the h0 0 1i direction. The initial position of the atoms are computed from the analytical solution of the anisotropic elasticity [33]. In the analytical solution, the length of the crack is chosen to be 80 times the atomic spacing. The surface energy, which is needed 2 in the anisotropic solution, is taken as 0:089 eV Å [32]. We also choose the stress intensity factor according to the Griffith criterion. Since the crack is of finite length, it tends to heal in the first few steps of the iteration. Therefore, we impose a uniform strain, e22 ¼ 0:02, to the initial displacement field. This is achieved by a continuation method: At the first step, the boundary condition is given by the initial displacement field. Then it is gradually changed to the actual displacement in 20 steps. After the initial preparation, we increase the external strain by 0.1% each step until e22 ¼ 0:03. In the computation, a dumbbell shaped region is selected around the crack. A long thin strip is used in the middle to include a few layers of atoms above and below the crack faces. Then at the two crack tips, we use two square regions. In the third direction, we use 10 layers of atoms with periodic boundary conditions. The total number of atoms in the atomistic region is 278,528, and for the boundary atom equation, we have 3092 elements and 3348 nodes. In Figs. 9 and 10 we show the solutions when e22 ¼ 0:025, including the full atomistic solutions, and the solutions of the ABEM method, and we have found good agreement for the crack tip positions. Notice that the displacement along the crack faces shows discrepancies. This can be attributed to the fact that the dispffiffiffi placement typically behaves like r as the distance r to the crack tip approaches zero, and a small error in the atomistic region can lead to large error in the surrounding region. For the computational cost, we compared the CPU time at loading step 2, and the results are listed in Table 2. To get an idea about the CPU time spent on each Schwarz iteration, we looked at one typical step of the overlapping Schwarz iteration. The solution procedure of the boundary atom equation only takes 0.027 s, and the solution of the problem in the atomistic region takes 42.4 s. For the non-overlapping method, we used c ¼ 0:09. The optimal value of the parameter seems to depend on the size of the surrounding region. For the overlapping method, the overlapping size is 5 times the atomic spacing.

0

−50 −100 −150 −200

−0.01

6. Conclusion and discussion

−0.02

In this paper, we have proposed a new reduction method for large scale molecular statics models of crystalline solids. By removing the atoms in the interior of the surrounding region, we have derived a smaller system of equations that only involve the atoms at the boundaries. We have also discussed several domain decomposition methods to couple the boundary condition with the atomistic regions containing local defects. The method can potentially be used to study mechanical properties of systems with many local defects, e.g. problems in fracture mechanics. Compared to traditional boundary integral approaches (e.g. [35] and many recent works), the current method uses fully atomistic model near crack tips to provide a more accurate description. The present formulation decouples the molecular statics models at the atomistic regions. In fact, after the boundary atom equations

−0.03

−250 −400

9

−300

−200

−100

0

100

200

300

400

−0.04

Fig. 8. Dislocation dipole: the error compared to a full atomistic solution.

5.2. Elliptic crack The following numerical test is conducted on a BCC Iron-a system. For the interatomic potential, we use the embedded atom potential [32]. The lattice parameter is a0 ¼ 2:866 Å with cut-off

10

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

6

Full Atomistic Non−overlapping Overlapping

5

4

3

2

1

0

−1 −500

−400

−300

−200

−100

0

100

200

300

400

500

Fig. 10. Comparison of the solutions along the crack face.

Table 2 Comparison of the computational cost.

CPU time (min)

Full atomistic

Non-overlapping

Overlapping

159

37

82

the applied loading changes. Therefore, an adaptive procedure is needed to identify the center (e.g. crack tip and dislocation cores) and re-define the atomistic regions as the computation proceeds. Here we given an example to illustrate how this can be done. We consider the fracture problem presented in section 5.2. In the top row of Fig. 11, we show the atoms in the atomistic region. We observe that as one increases the strain to the value 0.029, the crack has extended and the crack tip is approaching the boundary of the atomistic region. To allow further extension of the crack, we redefined the atomistic region by moving the square domain further to the right. This is done as follows: If the atom is also inside the previous atomistic region, then its position will be kept. Otherwise, the position can be determined from the evaluation step using the boundary atom equation. The results based on the redefined atomistic region are shown in the bottom row. This allows us to increase the strain to higher values.

Fig. 9. Comparison of solutions u2 on a {0 0 1} plane. From top to bottom: nonoverlapping domain decomposition; overlapping domain decomposition; full atomistic solution. The units are in Å.

are solved, the molecular statics models can be solved separately, which is ideal for parallelizing the numerical procedure. The present paper provides a basic framework to reduce the full atomistic model. There are several remaining issues which are currently being investigated. First the harmonic approximation introduced in the surrounding region can be further improved to include nonlinear terms. This can be done using a modified Newton’s method, and at each step, the linear system in the surrounding region, and the corresponding boundary atom equation, are solved to provide the search direction in the Newton’s method. Secondly, the boundary atom equations have been derived for simple lattices, and it is not yet clear how to treat complex lattices. Finally, in most practical cases, the atomistic regions might move as

Appendix A. Green’s functions In order to solve the boundary atom Eqs. (13) and (19), one needs to compute the lattice Green’s function given the distance between two atoms. A direct approach is to solve (8). Such a procedure has to be carried out within a finite domain with boundary conditions applied at the artificial boundary. Since in practice, the lattice Green’s function has to be evaluated for a wide range of interatomic distances, such a procedure is too expensive to implement. A much more efficient approach is to combine a Fourier integral method for small distances and an asymptotic representation in the far field. By taking Fourier transform, the lattice Green’s function can be expressed as an integral over the first Brillouin zone [38],

Gjjn ¼

1 jBj

Z B

DðnÞ1 cos n  ðRj  Rn Þdn;

ðA:1Þ

11

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

Fig. 11. Projection of atomic positions on the {0 0 1} plane. The figures at the bottom show the position of the atoms after the atomistic region is redefined. The situation for the other crack is similar.

where

DðnÞ ¼

P

Dk e

inRk

ðA:2Þ

k2N

is the dynamic matrix and B is the first Brillouin zone with volume jBj. This integral can be evaluated using the k-point method [28]. Several numerical issues have been addressed in [42]. In particular, for a large distance, i.e. when R ¼ Rj  Rn is large, the integrand is highly oscillatory and the results are not accurate. But in this regime the integral can be expanded, and the leading order approximation coincides with the continuum elastostatics Green’s functions. The error of the approximation is of the order jRj2 in two dimensions, and jRj3 in three dimensions. See [25] for more details on the derivation and high order expansions. The elastostatics Green’s function in 3D can be written as a onedimensional integral [3] and efficiently computed. For a twodimensional system, we will derive an explicit expression of the Green’s function. The following calculations are very similar to that of the Stroh’s formalism [36,40]. But we consider the elastic moduli directly obtained from the force constant of the interatomic potential, which may not have all the symmetry of the standard aniso-

tropic elasticity [1]. To be more specific, we seek the Green’s function of the continuum equation that corresponds to the linearized molecular statics model (10),

P @ 2 uj Eikjl ¼ fi ðXÞ; @X k @X l k;l

in R2 :

ðA:3Þ

Here Eikjl stands for the elastic constant defined in (24). The Green’s function is the fundamental solution of this equation. To quickly see the convergence of the lattice Green’s function to the elastostatics Green’s function, we notice that when R ¼ Rj  Rn becomes large, the integrand in the Fourier integral (A.1) is highly oscillatory. As a result, the main contribution comes from the origin. Near the origin, the dynamic matrix (A.2) can be expanded,

DðnÞ ¼ 

1P Dk ðn  Rk Þ2 þ oðjnj2 Þ ¼ V 0 KðnÞ þ oðjnj2 Þ: 2 k

Here we have used Eq. (24), and defined K to be a matrix with P entries K ij ¼ k;l Eikjl nk nl . We then have,

Gjjn 

1 V 0 jBj

Z

Rd

KðnÞ1 cos n  ðRj  Rn Þdn:

12

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13

0.35

Lattice Green’s function Anisotropic Elasticity Green’s function

0.3

GðX 1 ; X 2 Þ ¼

0

0.2

1

p

Imflnðz1 Þa1 a1 þ lnðz2 Þa2 a2 g:

ðA:8Þ

g/g

For more detailed discussions about derivations of this type, see [36,40].

0.15

References

0.1 0.05 0 0

2

4

6

8

r/a0

10

12

14

16

18

Fig. 12. The elastostatics Green’s function and the lattice Green’s function. Both functions have been normalized by G0j0 for better comparison.

Here V 0 jBj ¼ ð2pÞd . On the other hand, this approximation is exactly the Green’s function of the continuum model (A.3) expressed in a Fourier integral form. As an example, we computed the Green’s function for a body-centered cubic (BCC) system along the h111i direction. The results are plotted in Fig. 12, which shows the lattice Green’s function together with the elastostatics Green’s function. We observe that the lattice Green’s function starts to converge to the elastic Green’s function when the distance is beyond 6–10 atomic spacings. Now we return to (A.3). To express the fundamental solution in an explicit form, we choose a rectangular coordinate ðX 1 ; X 2 Þ, and seek the solution of (A.3) in the form of,

z ¼ X 1 þ lX 2 :

u ¼ lnðzÞa;

ðA:4Þ

For convenience, we define,

T ij ¼ Ei2j2 ; Rij ¼ Ei2j1 ;

ðA:5Þ

Q ij ¼ Ei1j1 : For simple lattices, all the three matrices are symmetric. In order for (A.4) to be a solution of the homogeneous equation, we must have,



l2 T þ 2lR þ Q a ¼ 0: 

ðA:6Þ 

Setting det l2 T þ 2lR þ Q ¼ 0, we find two pairs of roots, which have nonzero imaginary parts to ensure stability. Without loss of generality, we pick the two roots with positive imaginary parts,

l1 ¼ a1 þ ib1 ; l2 ¼ a2 þ ib2 : For the eigenvector a, it is chosen so that the normalization condition holds

aT ðR þ lTÞa ¼

1 2

for each root. This can be achieved if we pick a such that

a a¼

 1Þ d ¼ d0 b2 ðl2  l1 Þðl2  l for the first and second roots, respectively. d0 ¼ detðTÞ. The fundamental solution is then written in a compact form,

0.25



and

"  # T 22 l2 þ 2R22 l þ Q 22  T 12 l2 þ 2R12 l þ Q 12 1   : 2di  T 21 l2 þ 2R21 l þ Q 21 T 11 l2 þ 2R11 l þ Q 11 ðA:7Þ

Here

 2 Þ; d ¼ d0 b1 ðl1  l2 Þðl1  l

[1] N.W. Ashcroft, N.D. Mermin, Solid State Physics, Brooks Cole, 1976. [2] S. Badia, M. Parks, P. Bochev, M. Gunzburger, R. Lehoucq, On atomistic-tocontinuum coupling by blending, SIAM Multiscale Model. Simul. 7 (1) (2008) 381–406. [3] D.M. Barnett, The precise evaluation of derivatives of the anisotropic elastic Green’s functions, Phys. Stat. Sol. 49 (1972) 741. [4] J.T. Beale, M.C. Lai, A method for computing nearly singular integrals, SIAM J. Numer. Anal. 38 (2001) 1902–1925. [5] T. Belytschko, S.P. Xiao, Coupling methods for continuum model with molecular model, Inter. J. Multiscale Comput. Eng. 1 (2003) 1. [6] K.P. Bube, J.C. Strikwerda, Interior regularity estimates for elliptic systems of difference equations, SIAM J. Numer. Anal. 20 (1983) 653–670. [7] T.A. Cruse, Numerical solutions in three dimensional elastostatics, Int. J. Solids Struct. 5 (1969) 1259. [8] T.A. Cruse, An improved boundary-integral equation method for three dimensional elastic stress analysis, Comput. Struct. 4 (1974) 741–754. [9] M.S. Daw, M.I. Baskes, Embedded-atom method: derivation and application to impurities, surfaces, and other defects in metals, Phys. Rev. B 29 (1984) 6443. [10] W. E, J. Lu, J. Yang, Uniform accuracy of the quasicontinuum method, Phys. Rev. B 74 (2006) 214115. [11] W. E, P. Ming, Cauchy-Born rule and the stability of the crystalline solids: static problems, Arch. Rat. Mech. Anal. 183 (2009) 241–297. [12] F. Ercolessi, J.B. Adams, Interatomic potentials from first-principles calculations: the force-matching method, Europhys. Lett. 26 (1994) 583. [13] I. Szlufarska, Atomistic simulations of nanoindentation, Mater. Today 9 (2006) 42. [14] E.G. Karpov, H. Yu, H.S. Park, W.K. Liu, Q.J. Wang, D. Qian, Multiscale boundary conditions in crystalline solids: theory and application to nanoindentation, Int. J. Solids Struct. 43 (2006) 6359–6379. [15] C. Makridakis, C. Ortner, E. Süli, Stress-based atomistic/continuum coupling: a new variant of the quasicontinuum approximation, preprint, 2010. [16] J. Knap, M. Ortiz, An analysis of the quasicontinuum method, J. Mech. Phys. Solids 49 (2001) 1899–1923. [17] J. Knap, M. Ortiz, Effect of indenter-radius size on Au (0 0 1) nanoindentation, Phys. Rev. Lett. 90 (2003) 226102. [18] S. Kohlhoff, P. Gumbsch, H.F. Fischmeister, Crack propagation in b.c.c. crystals studied with a combined finite-element and atomistic model, Philos. Mag. A 64 (1991) 851–878. [19] X. Li, Efficient boundary condition for molecular statics models of solids, Phys. Rev. B 80 (2009) 104112. [20] D.C. Liu, J. Nocedal, On the limited memory BFGS method for large scale optimization, Math. Program. B 45 (1989) 503–528. [21] P. Ming, J.Z. Yang, Analysis of a one-dimensional nonlocal quasicontinuum method, SIAM Multiscale Model. Simul. 7 (2009) 1838–1875. [22] P. Ming, Error estimate of force-based quasicontinuum method, Comm. Math. Sci. 5 (2008) 1089–1095. [23] M. Dobson, M. Luskin, C. Ortner, Stability, instability, and error of the forcebased quasicontinuum approximation, Arch. Rat. Mech. Anal. 197 (2010) 179– 202. [24] P.G. Martinsson, Fast multiscale methods for lattice equations, Ph.D. Thesis, The University of Texas at Austin, 2002. [25] P.G. Martinsson, G.J. Rodin, Asymptotic expansions of lattice Green’s functions, Proc. R. Soc. Lond. A 458 (2002) 2609–2622. [26] S.N. Medyanik, E.G. Karpov, W.K. Liu, Domain reduction method for atomistic simulations, J. Comput. Phys. 218 (2006) 836–859. [27] R.E. Miller, E.B. Tadmor, A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods, Model. Simul. Mater. Sci. Eng. 17 (2009) 053001. [28] H.J. Monkhorst, J.D. Pack, Special points for Brillouin-zone integrations, Phys. Rev. B 13 (1976) 5188. [29] M. Mullins, M.A. Dokainish, Simulation of the (0 0 1) plane crack in a-iron employing a new boundary scheme, Philos. Mag. A 46 (1982) 771–787. [30] M.L. Parks, P.B. Bochev, R.B. Lehoucq, Connecting atomistic-to-continuum coupling and domain decomposition, SIAM Multiscale Model. Simul. 7 (2008) 362–380. [31] S. Rao, C. Hernandez, J.P. Simmons, T.A. Parthasarathy, C. Woodward, Green’s function boundary conditions in two-dimensional and three-dimensional atomistic simulations of dislocations, Philos. Mag. A 77 (1998) 231–256.

X. Li / Comput. Methods Appl. Mech. Engrg. 225–228 (2012) 1–13 [32] V. Shastry, D. Farkas, Molecular statics simulation of crack propagation in alpha-Fe using EAM potentials, in: Materials Research Society Symposium Proceedings, Boston, MA, United States, November 27–December 1, 1995, 1996, pp. 75–80. [33] G.C. Sih, H. Liebowitz, Fracture, An Advanced Treatise, Academic Press, 1968. [34] J.E. Sinclair, P.C. Gehlen, R.G. Hoagland, J.P. Hirth, Flexible boundary conditions and nonlinear geometric effects in atomic dislocation modeling, J. Appl. Phys. 49 (1978) 3890. [35] M.D. Snyder, T.A. Cruse, Boundary-integral equation analysis of cracked anisotropic plate, Int. J. Fract. 11 (1975) 315. [36] A.N. Stroh, Dislocations and cracks in anisotropic elasticity, Philos. Mag. 7 (1958) 625–646. [37] E.B. Tadmor, M. Ortiz, R. Phillips, Quasicontinuum analysis of defects in crystals, Philos. Mag. A 73 (1996) 1529–1563. [38] V.K. Tewary, Green-function method for lattice statics, Adv. Phys. 22 (1973) 757–810.

13

[39] V. Thomée, B. Westergren, Elliptic difference equations and interior regularity, Numer. Math. 11 (1968) 196–210. [40] T.C.T. Ting, Anisotropic elasticity: Theory and applications, Oxford, 1997. [41] A. Toselli, O. Widlund, Domain decomposition method – algorithms and theory, Springer, 2005. [42] D.R. Trinkle, Lattice Green function for extended defect calculations: computation and error estimation with long-range forces, Phys. Rev. B 78 (1) (2008) 014110. [43] C. Woodward, S.I. Rao, Flexible ab initio boundary conditions: simulating isolated dislocations in bcc Mo and Ta, Phys. Rev. Lett. 88 (2002) 216402. [44] A. Yavari, M. Ortiz, K. Bhattacharya, A theory of anharmonic lattice statics for analysis of defective crystals, J. Elasticity 86 (2007) 41–83. [45] S.J. Zhou, P.S. Lomdahl, A.F. Voter, B.L. Holian, Three-dimensional fracture via large-scale molecular dynamics, Engrg. Fract. Mech. 61 (1998) 173–187.