# Challenges in Atomistic-to-Continuum Coupling

Authors: Konstantin Fackeldey
Journal: Mathematical Problems in Engineering
 (2015)

---

## Abstract

This paper is concerned with the design, analysis, and implementation of concurrent coupling approaches where different (atomic and continuous) models are used simultaneously within a single simulation process. Thereby, several problems or pitfalls can happen, for example, the reflection of molecular movements at the “boundary” between the atomic and continuum regions which leads to an unphysical increase in energy in the atomic model. We investigate the problems with the aim of giving an introduction into this field and preventing errors for scientists starting their research towards multiscale methods.


---

## Body

## 1.Introduction

Classical molecular dynamics provides often an adequate description of a physical phenomenon but is computationally expensive. On the other hand, a depiction by a continuum mechanical model combined with the finite elements method is a smaller computational burden but does not go into atomistic details and can be considered less accurate. Multiscale methods which aim at combining fine scale method such as MD with coarse models such as continuum mechanics have been developed. In particular, for elastic systems these methods have reached a considerable interest, since they promise to be a powerful tool in many engineering problems, where only in a highly localized regionΩMD a precise model (fine scale) is needed and on a larger (surrounding) domain ΩCM a less accurate model (coarse scale) suffices.The development of different multiscale methods in different fields started decades ago and has been accelerated in the last few years. Along with this expansion several survey articles have been published in order to classify the multiscale methods by different aspects [1–5].These methods vary not only in scope and the underlying assumptions but also in their approach to broader questions such as a hierarchical and concurrent multiscale approach. In the first class, the computations are performed on each scale separately. Often, the scale coupling is done by transferring problem parameters; that is, the results obtained on one scale determine the parameters for the computational model on another scale [6, 7]. Thus, for instance, a continuum model can be derived from the atomic information [8]. Another approach is pursued in the concurrent coupling techniques, where the behavior at each length scale depends strongly on the others and an appropriate model is solved on each scale simultaneously, while smooth coupling between the scales is introduced. In the following, we confine our considerations to concurrent methods.


## 2.Demands on Multiscale Methods and Domain Decompositions

In the following, we show that Domain Decomposition (DD) methods serve as a good motivation for a classification of multiscale methods. To do so, we briefly explain some basic concepts from DD methods.The term DD [9, 10] is often used to describe a data distribution, in which the local data of each process corresponds topologically to a subdomain ΩS⊂Ω of the whole computational domain Ω.In the following, we use a different approach by defining DD as certain numerical methods that split the computational domain into two or more subdomains. Although DD methods have been developed for the purpose of achieving concurrency, they can be used in sequential as well as concurrent computations.In their origin DD techniques have been developed as a powerful iterative method for solving systems of algebraic equations stemming from the discretization of partial differential equations (PDE), that is, from a continuum description; see [11, 12] for an overview. Therein, DD is considered as a decomposition of the finite element space into a sum of subspaces. Then these subproblems are solved by a direct or iterative method. In a next step, projection operators are developed for the information transfer between the subspaces. As a matter of fact, the quality of the approximation on each subdomain depends on the corresponding properties of the approximation subspace. On the contrary, the DD method allows taking benefit of the presence of the subdomains in order to choose the discretization method, which is best adapted to the local behavior of the solution of the PDE which has to be approximated. Thus, the shape of the subdomains and their magnitude of overlap (interface) can be chosen problem dependent. Summing up, the choice of overlapping or nonoverlapping domains and the choice of the transfer operators deeply influence the performance of the DD method. In particular, the choice of an overlapping or nonoverlapping method directly influences the choice of the transfer operator.Let us now consider the continuum/molecular coupling. Then, the following three points have to be clarified:(i)
Definition of the domain.(ii)
Design of the coupling region.(iii)
Design of suitable transfer operators.Let us elucidate this DD motivated classification.


### 2.1. Definition of the Domain

Depending on the type of problem and on the considered domain, the regions with highly local interest have to be defined. For example, cracks and similar defects can involve a global simulation by finite elements and highly localized regions with strong deformations (e.g., crack tips) which are resolved by a molecular dynamics simulation. More precisely, the domainΩ⊂Rd is decomposed by(1)Ω=ΩMD∪ΩCM,where in ΩMD a fine resolution down to the atomistic scale is employed and in ΩCM a coarser representation by finite elements is applied. As mentioned previously, since the simulation on ΩMD is a higher computational burden, the size of ΩMD plays an important role in the overall performance of any multiscale scheme. As a particular advantage of DD methods the size and shape of the subdomains and thus the interface or handshake region Ξ=ΩMD∩ΩCM can be chosen arbitrarily. For the different types of interfaces we can distinguish between three cases (cf. Figure 1); these are the following:(i)
Interface Methods: The two scales are separated by an interface, which is often a manifold of dimension d-1 (Γ=Ω¯MD∩Ω¯CM).(ii)
Handshake Region Methods: The atomistic and the continuum part are matched in an overlap region Ξ=ΩMD∩ΩCM, Ξ⊂Rd, and measd(Ξ)>0.(iii)
Completely Overlapping Methods: The continuum mechanic description is on the whole domain Ω and the molecular part is a portion of it; that is, the atomistic simulation is everywhere accompanied by the continuum simulation (ΩMD∩ΩCM=ΩMD).Figure 1
Types of DDs: (a) nonoverlapping; (b) overlapping/partial overlapping; (c) overlapping/complete overlapping.
(a)
(b)
(c)The design of the coupling regions has consequences for the coupling between the molecular dynamics and the continuum mechanics. In nonoverlapping methods, the separation of the fine and the coarse scale is defined by the interface. In other words, the interface is the border between coarse and fine scale such that coexistence in some coupling region can be excluded.


### 2.2. Design of the Coupling Region

In the coupling region the two descriptions of the matter have to be matched. However, in general, the molecular dynamics is based on a description in some Euclidean space, whereas the continuum mechanics is often given in a weak sense, that is, in a function space. Thus, it is a priori not clear whether the coupling space is a subset of the Euclidean space or of a functions space.


### 2.3. Space and Transfer Operator

Here, two entirely different models (finite elements and molecular dynamics) are coupled. However, the relationship between their parameters is usually not direct and, thus, care must be taken in the construction of suitable transfer operators in order to relate them. Consequently, any operator has to deal with the incapability of the finite element discretization to resolve displacements down to the atomistic scale. Moreover, the choice of the underlying space for seamless coupling has to be made carefully.Summing up, the domainΩ of the problem under consideration is decomposed in a coarse and a fine scale part; that is,(2)Ω¯=ΩMD¯∪ΩCM¯.In the engineering literature, these methods are also sometimes termed as domain-decomposition techniques. In general we expect, for an efficient multiscale or atomistic-to-continuum method (AtC method), that diam⁡(ΩMD)≪diam⁡(ΩCM). One prominent example is the crack propagation, where a zone around the crack tip is simulated by molecular dynamics, whereas in the surrounding domain a finite element simulation is applied.The AtC methods have been applied with success in various engineering applications [13–15]. The different overlapping AtC methods vary in their design of the overlap and the treatment of the interaction between the continuum and the atomistic scale. One of the probably first overlapping AtC methods has been introduced by Mullins and Dokanish, where the stresses in the finite element simulation are computed by an interatomic potential [16]. Kohlhoff et al. [17] developed the FEAt method, where an atomistic model is surrounded by a finite element mesh with a small overlap region enforcing boundary condition on the atomistic and on the continuum domain. In the Quasicontinuum method [18–20], the total domain is considered and the atomistic region is systematically coarsened out by using kinematic constraints for predefined subsets of the domain. In the Bridging Domain method [21–24], the overlapping is defined via a convex combination of the two descriptions, such that in the handshake region the energy terms are not counted twice. Closely related to this approach is the Arlequin method [25], which has been initially introduced in the context of continuum-to-continuum blending (cf. Mortar methods [26]). This approach has been further investigated in the context of atomistic-to-continuum blending [27]. The general blending aspect has been investigated by [28]. In [29] the coupling has been realized in a function space concept allowing for a transfer between the coarse and the fine scale by using less degrees of freedom in the handshake region Ξ. The Bridging Scale method [30, 31], inspired by the variational multiscale method [32], considers the fine scale part as a subset of the coarse domain (ΩMD⊂ΩCM). In order to identify the coarse and the fine scale part in the overlapping region, the total displacement field is decomposed into a fine scale and a coarse scale part. The coarse scale part is the projection of the total displacement onto the finite elements, and the fine scale part is the residue of this projection; that is, the fine scale is that part of the total displacement field which has not been captured by the projection onto the finite element basis. The AtC coupling has also been examined from the mathematical point of view, for example, [6, 33–35]. Recently, coupling methods have been proposed, where the microscale is a canonical ensemble, that is, AtC with finite temperature [24, 36–38].As a matter of fact, coupling of different physical models is not free from artifacts. One well known and well described phenomenon is the spurious reflections on the interface between the pure atomistic and the mixed region [39, 40]. Here, in this paper, we elucidate the problems or challenges, when coupling continuum mechanics with atomistic models. This paper is structured as follows. In Section 3, we introduce the mathematical setting for the fine scale and for the coarse scale; thereby, we work out that the natural space for continuum mechanics is a function space, whereas for the molecular dynamics the underlying space is the Euclidean space. In Section 4, we derive the different dispersion relations for the atomistic and the continuum scale; thereby, we see that they differ for large wave numbers. In Section 5, we show that a larger mesh size of the finite elements, which is not the same as the distance between the particles, also can cause reflections; thereby, we show a relationship with Shannon’s sampling theorem.


## 2.1. Definition of the Domain

Depending on the type of problem and on the considered domain, the regions with highly local interest have to be defined. For example, cracks and similar defects can involve a global simulation by finite elements and highly localized regions with strong deformations (e.g., crack tips) which are resolved by a molecular dynamics simulation. More precisely, the domainΩ⊂Rd is decomposed by(1)Ω=ΩMD∪ΩCM,where in ΩMD a fine resolution down to the atomistic scale is employed and in ΩCM a coarser representation by finite elements is applied. As mentioned previously, since the simulation on ΩMD is a higher computational burden, the size of ΩMD plays an important role in the overall performance of any multiscale scheme. As a particular advantage of DD methods the size and shape of the subdomains and thus the interface or handshake region Ξ=ΩMD∩ΩCM can be chosen arbitrarily. For the different types of interfaces we can distinguish between three cases (cf. Figure 1); these are the following:(i)
Interface Methods: The two scales are separated by an interface, which is often a manifold of dimension d-1 (Γ=Ω¯MD∩Ω¯CM).(ii)
Handshake Region Methods: The atomistic and the continuum part are matched in an overlap region Ξ=ΩMD∩ΩCM, Ξ⊂Rd, and measd(Ξ)>0.(iii)
Completely Overlapping Methods: The continuum mechanic description is on the whole domain Ω and the molecular part is a portion of it; that is, the atomistic simulation is everywhere accompanied by the continuum simulation (ΩMD∩ΩCM=ΩMD).Figure 1
Types of DDs: (a) nonoverlapping; (b) overlapping/partial overlapping; (c) overlapping/complete overlapping.
(a)
(b)
(c)The design of the coupling regions has consequences for the coupling between the molecular dynamics and the continuum mechanics. In nonoverlapping methods, the separation of the fine and the coarse scale is defined by the interface. In other words, the interface is the border between coarse and fine scale such that coexistence in some coupling region can be excluded.


## 2.2. Design of the Coupling Region

In the coupling region the two descriptions of the matter have to be matched. However, in general, the molecular dynamics is based on a description in some Euclidean space, whereas the continuum mechanics is often given in a weak sense, that is, in a function space. Thus, it is a priori not clear whether the coupling space is a subset of the Euclidean space or of a functions space.


## 2.3. Space and Transfer Operator

Here, two entirely different models (finite elements and molecular dynamics) are coupled. However, the relationship between their parameters is usually not direct and, thus, care must be taken in the construction of suitable transfer operators in order to relate them. Consequently, any operator has to deal with the incapability of the finite element discretization to resolve displacements down to the atomistic scale. Moreover, the choice of the underlying space for seamless coupling has to be made carefully.Summing up, the domainΩ of the problem under consideration is decomposed in a coarse and a fine scale part; that is,(2)Ω¯=ΩMD¯∪ΩCM¯.In the engineering literature, these methods are also sometimes termed as domain-decomposition techniques. In general we expect, for an efficient multiscale or atomistic-to-continuum method (AtC method), that diam⁡(ΩMD)≪diam⁡(ΩCM). One prominent example is the crack propagation, where a zone around the crack tip is simulated by molecular dynamics, whereas in the surrounding domain a finite element simulation is applied.The AtC methods have been applied with success in various engineering applications [13–15]. The different overlapping AtC methods vary in their design of the overlap and the treatment of the interaction between the continuum and the atomistic scale. One of the probably first overlapping AtC methods has been introduced by Mullins and Dokanish, where the stresses in the finite element simulation are computed by an interatomic potential [16]. Kohlhoff et al. [17] developed the FEAt method, where an atomistic model is surrounded by a finite element mesh with a small overlap region enforcing boundary condition on the atomistic and on the continuum domain. In the Quasicontinuum method [18–20], the total domain is considered and the atomistic region is systematically coarsened out by using kinematic constraints for predefined subsets of the domain. In the Bridging Domain method [21–24], the overlapping is defined via a convex combination of the two descriptions, such that in the handshake region the energy terms are not counted twice. Closely related to this approach is the Arlequin method [25], which has been initially introduced in the context of continuum-to-continuum blending (cf. Mortar methods [26]). This approach has been further investigated in the context of atomistic-to-continuum blending [27]. The general blending aspect has been investigated by [28]. In [29] the coupling has been realized in a function space concept allowing for a transfer between the coarse and the fine scale by using less degrees of freedom in the handshake region Ξ. The Bridging Scale method [30, 31], inspired by the variational multiscale method [32], considers the fine scale part as a subset of the coarse domain (ΩMD⊂ΩCM). In order to identify the coarse and the fine scale part in the overlapping region, the total displacement field is decomposed into a fine scale and a coarse scale part. The coarse scale part is the projection of the total displacement onto the finite elements, and the fine scale part is the residue of this projection; that is, the fine scale is that part of the total displacement field which has not been captured by the projection onto the finite element basis. The AtC coupling has also been examined from the mathematical point of view, for example, [6, 33–35]. Recently, coupling methods have been proposed, where the microscale is a canonical ensemble, that is, AtC with finite temperature [24, 36–38].As a matter of fact, coupling of different physical models is not free from artifacts. One well known and well described phenomenon is the spurious reflections on the interface between the pure atomistic and the mixed region [39, 40]. Here, in this paper, we elucidate the problems or challenges, when coupling continuum mechanics with atomistic models. This paper is structured as follows. In Section 3, we introduce the mathematical setting for the fine scale and for the coarse scale; thereby, we work out that the natural space for continuum mechanics is a function space, whereas for the molecular dynamics the underlying space is the Euclidean space. In Section 4, we derive the different dispersion relations for the atomistic and the continuum scale; thereby, we see that they differ for large wave numbers. In Section 5, we show that a larger mesh size of the finite elements, which is not the same as the distance between the particles, also can cause reflections; thereby, we show a relationship with Shannon’s sampling theorem.


## 3.Challenge 1: Atomistic versus Continuum Description

When coupling continuum mechanics and molecular dynamics one has to account for the different fundamentals of each theory. Whereas in the atomistic theory every individual particle is considered, in the continuum mechanics an average is considered. In this section we briefly introduce each theory.


### 3.1. Molecular Dynamics

The material behavior on the microscale is modeled by means of an isolated system of atoms or molecules of a crystalline solid. We identify each of the atoms in their reference position with a pointXα∈Rd, α∈A, where A is an index set. Under the influence of external and internal forces, the atoms are displaced in space. The position X^α of the αth atom in a deformed configuration is then given as(3)X^α=Xα+qα,where qα is the displacement of atom α.The atomic displacementsq=(qα)α∈A are assumed to obey Newton’s law of motion(4)MAq¨=Finternal+Fexternal,where Finternal and Fexternal are the internal and external forces. With each atom α, we associate the mass mα>0, such that(5)MA=diag⁡mαIdRd×dα∈Ais the mass matrix on the microscale. If the system is conservative a potential V is given; the internal forces acting on a conservative system can be obtained as Finternal=-∇X^V(X^).


### 3.2. The Continuum Mechanics Setting

The most common way to classical elasticity is the axiomatic introduction of the relevant quantities and the equations of motion (see, e.g., [41]). Therein, the atomistic structure of the solid is neglected and the motion of the body is described by functions. Let us consider Ω⊂Rd. In classical elasticity it is assumed that the motion, due to external and internal forces, is given by the mapping(6)ϕ:Ω¯×0,T⟶Rd,such that the deformed configuration at time t of the body is given by(7)Ωt=ϕΩ,tt≥0,with the displacement u(t)=ϕ(X,t)-X, X∈Ω¯, t≥0.In the forthcoming we denote points of the deformed configuration by(8)xt≔ϕX,tand write x instead of x(t) and u instead of u(t) whenever possible. Furthermore, it is assumed that ϕ is locally injective in Ω, sufficiently smooth, and orientation preserving, which can be expressed by the pointwise inequality(9)det⁡∂ϕi∂XjX,t>0,X∈Ω¯,t≥0,i,j=1,…,3.Let us denote the deformation gradient by F; that is,(10)F≔∂ϕi∂Xji,j=1,…,d.The mass is given by(11)mΩt=∫Ωtϱx,tdx,where ϱ(x,t) is the mass density. Here, we assume that the mass is constant in time, which means that there are no material flows through the boundary of a material subdomain and we do not consider mass to energy conversions. For the sake of simplicity the density is given as a function on Rd×[0,T] by extending ϱ(·,t)∈L∞(Ω,[0,∞]) by zero outside of Ωt.The change of square length of the line segmentdx=FdX due to the deformation is given by the Green strain tensor: E=1/2(FTF-I)=1/2(∇uT+∇u+∇uT∇u). Moreover, the Cauchy stress is defined by σn(x)dS=dF=fdS, where n is the normal and f is the traction on any surface segment dS. For computational and analytical reasons, we use the Lagrangian description and thus define the firstPiola-Kirchhoff tensorP, which is a second-order tensor, given by(12)P:0,T×Ω⟶R3×3,PXndSref=dF=fdS,where dSref is the surface segment in the undeformed (reference) configuration. The governing equations are then derived by conservation of mass, momentum, and angular momentum. They are given by(13)ϱ∂v∂t=div⁡P+ϱb,Pun=f,onΓN,where v=u˙=∂u/∂t. Here, ϱb is the body force density and v∈C1([0,T],C0(Ω)) and P∈C0([0,T],C1(Ω)). Let us remark that, for the sake of simplicity, we neglected the boundary conditions in our introduction. However, depending on the material model and the pointwise smoothness of the boundary ∂Ω of Ω, the solutions for u may not be smooth enough to fulfill (13).Thus, we introduce the weak formalism. To do so, letΩ⊂Rd be measurable. We denote by Lp(Ω) the space of p-Lebesgue-integrable functions on Ω. Moreover, we denote the L2 scalar product by (·,·)L2(Ω) and the norm(14)uL2Ω≔∫Ωu2dΩ1/2.Let α be a multi-index with α1≔∑i=1d‍|αi|. Then, for a bounded domain Ω and for m∈N, 1≤p<∞, we denote by Hm,p(Ω) the Sobolev space, given by(15)
H
m
,
p
Ω
≔
u
∈
L
p
Ω
∣
there  exists
∂
α
u
,
∂
α
u
∈
L
p
Ω
∀
α
with
α
≤
mand equipped with the norm(16)uHm,pΩ≔∑α1≤m∫Ω∂αupdΩ1/p.We assume that Ω is a polygonal domain and ∂Ω=ΓN∪˙ΓD is the disjoint union of the portions ΓD and ΓN. Moreover, we assume that b∈C0([0,T],L2(Ω)), f∈C0([0,T],L2(ΓN)), uD∈C0([0,T],H1/2(ΓD)), and the initial displacements u0∈H1(Ω). Then for the weak formulation of (13) we seek u∈C2([0,T],H1(Ω)) such that u(·,0)=u0, u|ΓD=uD, and(17)∫Ωϱu¨·dΩ+∫ΩPu:∇vdΩ=∫Ωϱb·vdΩ+∫ΓNf·vdSfor all v∈C0([0,T];H1(Ω)) vanishing on ΓD, where A:B denotes tr(ATB). Comparing (4) and (17) it becomes clear that there are different methods to combine these two expressions. The interpretation of the molecular description as points in the function space L2 can be realized in different manners. Many of the existing AtC methods constrain the atomistic displacement in the handshake region Ξ by an interpolation of the continuum displacement (e.g., [21, 22, 30]). More precisely, the continuum displacement is approximated by finite elements and the atomistic displacement is projected onto the finite element basis functions. Thus, a matching between the two scales is realized in a pointwise sense. In [29, 42], this matching in the handshake region is realized in a function space oriented sense, by embedding the molecular displacement into the L2(Ξ).


## 3.1. Molecular Dynamics

The material behavior on the microscale is modeled by means of an isolated system of atoms or molecules of a crystalline solid. We identify each of the atoms in their reference position with a pointXα∈Rd, α∈A, where A is an index set. Under the influence of external and internal forces, the atoms are displaced in space. The position X^α of the αth atom in a deformed configuration is then given as(3)X^α=Xα+qα,where qα is the displacement of atom α.The atomic displacementsq=(qα)α∈A are assumed to obey Newton’s law of motion(4)MAq¨=Finternal+Fexternal,where Finternal and Fexternal are the internal and external forces. With each atom α, we associate the mass mα>0, such that(5)MA=diag⁡mαIdRd×dα∈Ais the mass matrix on the microscale. If the system is conservative a potential V is given; the internal forces acting on a conservative system can be obtained as Finternal=-∇X^V(X^).


## 3.2. The Continuum Mechanics Setting

The most common way to classical elasticity is the axiomatic introduction of the relevant quantities and the equations of motion (see, e.g., [41]). Therein, the atomistic structure of the solid is neglected and the motion of the body is described by functions. Let us consider Ω⊂Rd. In classical elasticity it is assumed that the motion, due to external and internal forces, is given by the mapping(6)ϕ:Ω¯×0,T⟶Rd,such that the deformed configuration at time t of the body is given by(7)Ωt=ϕΩ,tt≥0,with the displacement u(t)=ϕ(X,t)-X, X∈Ω¯, t≥0.In the forthcoming we denote points of the deformed configuration by(8)xt≔ϕX,tand write x instead of x(t) and u instead of u(t) whenever possible. Furthermore, it is assumed that ϕ is locally injective in Ω, sufficiently smooth, and orientation preserving, which can be expressed by the pointwise inequality(9)det⁡∂ϕi∂XjX,t>0,X∈Ω¯,t≥0,i,j=1,…,3.Let us denote the deformation gradient by F; that is,(10)F≔∂ϕi∂Xji,j=1,…,d.The mass is given by(11)mΩt=∫Ωtϱx,tdx,where ϱ(x,t) is the mass density. Here, we assume that the mass is constant in time, which means that there are no material flows through the boundary of a material subdomain and we do not consider mass to energy conversions. For the sake of simplicity the density is given as a function on Rd×[0,T] by extending ϱ(·,t)∈L∞(Ω,[0,∞]) by zero outside of Ωt.The change of square length of the line segmentdx=FdX due to the deformation is given by the Green strain tensor: E=1/2(FTF-I)=1/2(∇uT+∇u+∇uT∇u). Moreover, the Cauchy stress is defined by σn(x)dS=dF=fdS, where n is the normal and f is the traction on any surface segment dS. For computational and analytical reasons, we use the Lagrangian description and thus define the firstPiola-Kirchhoff tensorP, which is a second-order tensor, given by(12)P:0,T×Ω⟶R3×3,PXndSref=dF=fdS,where dSref is the surface segment in the undeformed (reference) configuration. The governing equations are then derived by conservation of mass, momentum, and angular momentum. They are given by(13)ϱ∂v∂t=div⁡P+ϱb,Pun=f,onΓN,where v=u˙=∂u/∂t. Here, ϱb is the body force density and v∈C1([0,T],C0(Ω)) and P∈C0([0,T],C1(Ω)). Let us remark that, for the sake of simplicity, we neglected the boundary conditions in our introduction. However, depending on the material model and the pointwise smoothness of the boundary ∂Ω of Ω, the solutions for u may not be smooth enough to fulfill (13).Thus, we introduce the weak formalism. To do so, letΩ⊂Rd be measurable. We denote by Lp(Ω) the space of p-Lebesgue-integrable functions on Ω. Moreover, we denote the L2 scalar product by (·,·)L2(Ω) and the norm(14)uL2Ω≔∫Ωu2dΩ1/2.Let α be a multi-index with α1≔∑i=1d‍|αi|. Then, for a bounded domain Ω and for m∈N, 1≤p<∞, we denote by Hm,p(Ω) the Sobolev space, given by(15)
H
m
,
p
Ω
≔
u
∈
L
p
Ω
∣
there  exists
∂
α
u
,
∂
α
u
∈
L
p
Ω
∀
α
with
α
≤
mand equipped with the norm(16)uHm,pΩ≔∑α1≤m∫Ω∂αupdΩ1/p.We assume that Ω is a polygonal domain and ∂Ω=ΓN∪˙ΓD is the disjoint union of the portions ΓD and ΓN. Moreover, we assume that b∈C0([0,T],L2(Ω)), f∈C0([0,T],L2(ΓN)), uD∈C0([0,T],H1/2(ΓD)), and the initial displacements u0∈H1(Ω). Then for the weak formulation of (13) we seek u∈C2([0,T],H1(Ω)) such that u(·,0)=u0, u|ΓD=uD, and(17)∫Ωϱu¨·dΩ+∫ΩPu:∇vdΩ=∫Ωϱb·vdΩ+∫ΓNf·vdSfor all v∈C0([0,T];H1(Ω)) vanishing on ΓD, where A:B denotes tr(ATB). Comparing (4) and (17) it becomes clear that there are different methods to combine these two expressions. The interpretation of the molecular description as points in the function space L2 can be realized in different manners. Many of the existing AtC methods constrain the atomistic displacement in the handshake region Ξ by an interpolation of the continuum displacement (e.g., [21, 22, 30]). More precisely, the continuum displacement is approximated by finite elements and the atomistic displacement is projected onto the finite element basis functions. Thus, a matching between the two scales is realized in a pointwise sense. In [29, 42], this matching in the handshake region is realized in a function space oriented sense, by embedding the molecular displacement into the L2(Ξ).


## 4.Challenge II: The Dispersion Relation

In the foregoing section, we introduced the molecular description (fine scale) and the continuum description (coarse scale). Here, in this section we consider the velocity at which waves propagate fromΩMD to ΩCM or from ΩCM to ΩMD. In detail, we want to guarantee that this velocity is conserved when the model and consequently the discretization change (i.e., cross ΩMD∩ΩCM), since otherwise artefacts at the common boundary can spoil the simulation. In Figure 2, an example of a wave entering a region of slower (faster) wave propagation speed due to a change of the spatial discretization size is shown: The sketch reveals that the difference in the dispersion relation leads to a partial reflection of the wave. To do so, we derive the dispersion relation, giving the dependence of the frequency on the wave number. In this context, we then see that continuum mechanics, as introduced here, with a specific constitutive relation, only represents the linear part of this relation. This difference has fundamental impact for atomistic-to-continuum coupling schemes. More precisely, in an atomistic setting, the dispersion relation, giving the dependence of the frequency on the wave number, is nonlinear whereas in the continuum setting this relation is linear. As a consequence, in the atomistic region the waves with higher frequency move slower than waves with a lower frequency. A typical behavior of such a nonlinear dispersion relation is shown in Figure 3. On the contrary, in the continuum setting, both types of waves have the same velocity. When coupling continuum mechanics and molecular dynamics, this fact has to be taken into account.Figure 2
(a) Smooth wave entering from a faster wave speed region (black line) into a slower wave speed region (blue line). (b) Smooth wave entering from a slower wave speed region (black line) into a faster wave speed region (blue line).
(a)
(b)Figure 3
(a) Wave package in1d in initial position. (b) Propagated initial wave: the waves with lower frequencies move faster.
(a)
(b)In order to derive the dispersion relation, we first briefly introduce the traveling waves in crystals and show then the dispersion relation for the molecular dynamics setting as well as for the continuum setting. For the sake of simplicity, we confine our discussion to the1d case.


### 4.1. Traveling Waves in Crystals

Recalling thatq(X,t) denotes the deviation of the atoms from their starting configuration a harmonic wave in the atomistic model is of the form(18)qX,t=q~eikX-ωt,where q~ is the amplitude, kX-ωt is the phase, k is the wave vector, and ω is the angular frequency. We also define the wave number κ=k and the wave length λ=1/κ. In the forthcoming, we consider two different kinds of velocities. Thephase velocity  vph is defined as(19)vph=ωκ.Moreover, for the other kind of velocity, we define the wave package as a wave whose amplitude is only in a bounded domain nonzero. The velocity of a wave package is thegroup velocity, defined as(20)vgr=∂ω∂κ.We can decide between longitudinal and transverse vibrations. An example of a1d chain in 2d is given by Figure 4.Figure 4
An example of (a) longitudinal and (b) transversal waves of a1d chain in 2d. The arrows show the displacement direction of the particles.
(a)
(b)Remark 1.
In3d the wave vector k has three components and points into the direction of propagation. The atomistic displacements associated with a wave q(X,t) where X is the equilibrium position of an atom are simultaneously three-dimensional vectors. These displacement vectors may be parallel to k (longitudinal), perpendicular to k (transversal), or along a direction, that is, not directly related to the direction of k.
For the sake of simplicity, we assume that the time dependent motion of the atoms is a linear superposition of harmonic waves. From linearized equation of motion, necessary conditions for nonlinear systems such as typical MD simulations can then be given. In1d each X can be rewritten as αr0, where α∈Z and r0 is the atomistic spacing. Then, for each α, we have the superposition(21)qαt=∑kq~keikαr0-ωtand we can understand the behavior of the solution by examining the harmonic solutions independently; that is, qk=0 for all but one k.Let us consider a mass spring system with lattice spacingr0, mass m, and spring constant K. Assuming that for this harmonic system a Hamiltonian is given, the respective equations of motion of HMD can be stated as(22)md2qαdt2=Kqα+1+qα-1-2qαwhich then gives us the dispersion relation for the molecular case (cf. [43, 44])(23)ω2=2Kmsin2⁡kr02.One can see easily thatω(k)=ωk+2πn/r0 for any n∈Z. We thus only consider the case of k∈(-π/r0,π/r0), which is known as the first Brillouin zone. Note that the open interval is chosen, as for k=π/r0 the group velocity vgr equals zero, meaning that the solution is a standing wave. This physical phenomenon is also known as Bragg-Reflection [43, 44]. For symmetry reasons it is sufficient to be restricted to k=κ>0.


### 4.2. The Dispersion Relation in the Continuum Setting

The continuum mechanics and the molecular dynamics stem from the same basic physical laws, namely, the Hamiltonian or Lagrangian principles. In this section, we consider the dispersion relation for the continuum mechanics. Let us therefore reconsider the molecular dispersion relation given by (23) for small wave numbers κ. More precisely let κ∗ be such that(24)sin⁡κ∗r0=κ∗r0+oκ∗r0;then for all 0<κ<κ∗ we have(25)ωκ→0=12r0Km1/2κ.Moreover, the phase velocity vphCM is then given by(26)vphCM=ωκ=r0Km1/2.A connection to the macroscopic elastic properties can be given by compressing the one-dimensional chain from Section 4.1 giving rise to the strain e, such that the average distance between the atoms becomes rc=r0(1-e), where e≪1.With respect to the harmonic energy, given byEharm,1d=1/2Csqα-qα+s2, the energy of the strained chain becomes(27)Eharm,1d=Vr+12NKr0-rc2.Thus, the strain energy, which can be considered as the extra energy per atom, is given by(28)Estrain,1d=12Kr0-rc2=12Kr02e2=12K^e2,where K^≔Kr02 is the elastic constant. From the velocity given by (26) we have(29)vphCM2=ω2κ2=Kr02m=K^m.Summing up, for the long wavelength limit (λ=1/κ, κ≪1), we obtain(30)ω2=K^κ22m=Kκ2r022m.The meaning of the different dispersion relations (23) and (30) becomes clear when the phase velocity vph and the group velocity vgr are considered. It can easily be seen that the molecular phase and the group velocity are given by(31)vphMD=2K/msin⁡κr0/2κ,vgrMD=K2mr0cos⁡κr02,whereas their continuum counterparts are given by(32)vphCM=vgrCM=K2mr0.Thus, for small κ we have(33)vphMD≈vphCM,vgrMD≈vgrCM.In contrast, for large wave numbers, that is, κ≈π/2, we have cos⁡(κr0/2)≈0, which implies that in the atomistic setting waves with high frequencies propagate slower than waves with low frequencies, which is not the case in the continuum since its dispersion relation is linear (cf. (30)). This difference in the velocities of the respective models leads to spurious reflections on the interface between the atomistic and continuum description. We have thus motivated that the spurious reflection are not only an artefact of the discretization but a consequence of the difference in the models.Therefore, each simulation or discretization has to account for this fact. To overcome this difference, absorbing boundary conditions or nonreflecting boundary conditions, like the time history kernel method combined with a scale decomposition [34], the variational boundary condition [45], blending methods [46], and other methods, for example, [24, 29], have been introduced.


## 4.1. Traveling Waves in Crystals

Recalling thatq(X,t) denotes the deviation of the atoms from their starting configuration a harmonic wave in the atomistic model is of the form(18)qX,t=q~eikX-ωt,where q~ is the amplitude, kX-ωt is the phase, k is the wave vector, and ω is the angular frequency. We also define the wave number κ=k and the wave length λ=1/κ. In the forthcoming, we consider two different kinds of velocities. Thephase velocity  vph is defined as(19)vph=ωκ.Moreover, for the other kind of velocity, we define the wave package as a wave whose amplitude is only in a bounded domain nonzero. The velocity of a wave package is thegroup velocity, defined as(20)vgr=∂ω∂κ.We can decide between longitudinal and transverse vibrations. An example of a1d chain in 2d is given by Figure 4.Figure 4
An example of (a) longitudinal and (b) transversal waves of a1d chain in 2d. The arrows show the displacement direction of the particles.
(a)
(b)Remark 1.
In3d the wave vector k has three components and points into the direction of propagation. The atomistic displacements associated with a wave q(X,t) where X is the equilibrium position of an atom are simultaneously three-dimensional vectors. These displacement vectors may be parallel to k (longitudinal), perpendicular to k (transversal), or along a direction, that is, not directly related to the direction of k.
For the sake of simplicity, we assume that the time dependent motion of the atoms is a linear superposition of harmonic waves. From linearized equation of motion, necessary conditions for nonlinear systems such as typical MD simulations can then be given. In1d each X can be rewritten as αr0, where α∈Z and r0 is the atomistic spacing. Then, for each α, we have the superposition(21)qαt=∑kq~keikαr0-ωtand we can understand the behavior of the solution by examining the harmonic solutions independently; that is, qk=0 for all but one k.Let us consider a mass spring system with lattice spacingr0, mass m, and spring constant K. Assuming that for this harmonic system a Hamiltonian is given, the respective equations of motion of HMD can be stated as(22)md2qαdt2=Kqα+1+qα-1-2qαwhich then gives us the dispersion relation for the molecular case (cf. [43, 44])(23)ω2=2Kmsin2⁡kr02.One can see easily thatω(k)=ωk+2πn/r0 for any n∈Z. We thus only consider the case of k∈(-π/r0,π/r0), which is known as the first Brillouin zone. Note that the open interval is chosen, as for k=π/r0 the group velocity vgr equals zero, meaning that the solution is a standing wave. This physical phenomenon is also known as Bragg-Reflection [43, 44]. For symmetry reasons it is sufficient to be restricted to k=κ>0.


## 4.2. The Dispersion Relation in the Continuum Setting

The continuum mechanics and the molecular dynamics stem from the same basic physical laws, namely, the Hamiltonian or Lagrangian principles. In this section, we consider the dispersion relation for the continuum mechanics. Let us therefore reconsider the molecular dispersion relation given by (23) for small wave numbers κ. More precisely let κ∗ be such that(24)sin⁡κ∗r0=κ∗r0+oκ∗r0;then for all 0<κ<κ∗ we have(25)ωκ→0=12r0Km1/2κ.Moreover, the phase velocity vphCM is then given by(26)vphCM=ωκ=r0Km1/2.A connection to the macroscopic elastic properties can be given by compressing the one-dimensional chain from Section 4.1 giving rise to the strain e, such that the average distance between the atoms becomes rc=r0(1-e), where e≪1.With respect to the harmonic energy, given byEharm,1d=1/2Csqα-qα+s2, the energy of the strained chain becomes(27)Eharm,1d=Vr+12NKr0-rc2.Thus, the strain energy, which can be considered as the extra energy per atom, is given by(28)Estrain,1d=12Kr0-rc2=12Kr02e2=12K^e2,where K^≔Kr02 is the elastic constant. From the velocity given by (26) we have(29)vphCM2=ω2κ2=Kr02m=K^m.Summing up, for the long wavelength limit (λ=1/κ, κ≪1), we obtain(30)ω2=K^κ22m=Kκ2r022m.The meaning of the different dispersion relations (23) and (30) becomes clear when the phase velocity vph and the group velocity vgr are considered. It can easily be seen that the molecular phase and the group velocity are given by(31)vphMD=2K/msin⁡κr0/2κ,vgrMD=K2mr0cos⁡κr02,whereas their continuum counterparts are given by(32)vphCM=vgrCM=K2mr0.Thus, for small κ we have(33)vphMD≈vphCM,vgrMD≈vgrCM.In contrast, for large wave numbers, that is, κ≈π/2, we have cos⁡(κr0/2)≈0, which implies that in the atomistic setting waves with high frequencies propagate slower than waves with low frequencies, which is not the case in the continuum since its dispersion relation is linear (cf. (30)). This difference in the velocities of the respective models leads to spurious reflections on the interface between the atomistic and continuum description. We have thus motivated that the spurious reflection are not only an artefact of the discretization but a consequence of the difference in the models.Therefore, each simulation or discretization has to account for this fact. To overcome this difference, absorbing boundary conditions or nonreflecting boundary conditions, like the time history kernel method combined with a scale decomposition [34], the variational boundary condition [45], blending methods [46], and other methods, for example, [24, 29], have been introduced.


## 5.Challenge III: Discretization of the Coarse Scale

So far, we have shown that the continuum and the molecular scale have different dispersion relations, which carries over to a difference in the velocities of waves. However, for long wavelengths these differences are insignificant. For the numerical dispersion, that is, the case when the continuum is discretized, the mesh size and the shape of the basis functions have to be accounted for.In order to approximate the continuous displacement fieldu, we now employ a finite element discretization of lower order. Let Th denote a mesh with mesh size parameter h>0, such that the family Thh is shape regular.By using Lagrangian conforming finite elements of first order (P1) for the displacements u and denoting the set of all nodes of Th by Nh, the finite element space Vh⊂H1(Ω) is spanned by the nodal basis(34)Vh=spanp∈Nhψph.The Lagrangian basis functions ψph∈Vh are uniquely characterized by the Kronecker-delta property(35)ψphq=δpq,p,q∈Nh,where δpq is the Kronecker-delta. Any function uh∈Vh(Ω) can uniquely be written as(36)uh=∑p∈Nhu¯pψph,where (up)p∈Nh∈Rd·Nh, u¯p∈Rd, is the coefficient vector. We can identify each element of Vh with its coefficient vector (up)p∈Nh. In the forthcoming, we omit the superscript h whenever possible. In space, (17) gives a nonlinear ordinary differential equation in RdNh which can be solved by, for example, the Verlet algorithm (e.g., [47]). In other words, we do a spatial discretization followed by a time discretization (method of lines).Analogously to the atomistic case we consider a finite element approximation for the continuum Hamiltonian of a harmonic system in1d. By assuming that the material of the atomistic as well as of the continuum model has the same parameters, the respective equations of motion with a spatial discretization can be stated as (37)md2updt2=Kr02h2up+1+up-1-2up.Hence, if the mesh size h equals the atomistic spacing r0 we have that (22) and (37) coincide.The harmonic wave given by (18) can be interpreted as a signal and the discretization of it as a sampling, which leads us directly to Shannon’s sampling theorem. It states that the sampling frequency of a signal should be at least twice the largest frequency contained in the signal.In our context, the sampling rate is1/h and the maximum wave length is given by λ such that Shannon’s theorem becomes(38)1h>21λ⟺h<12λ.This shows that the mesh size h has to be chosen small enough to capture the total atomistic displacement. However, doing so would dramatically increase the computational cost.Thus, in most AtC methods, smallλ waves on the continuum scale are not represented by finite elements but by other techniques (e.g., [30, 45]).Let us assume that the solution of (37) is given componentwise by(39)upt=u~eikxp-ωt,where u~ is the amplitude and ω is the frequency.Closely related to the Shannon sampling condition, the mesh size also influences the velocities of the wave. To see this let us consider the numerical dispersion relation for the finite element discretization(40)ω2κ=2lmsin2⁡κh2l≔Kr02h2.Comparing (23) and (40) one can see that both relations only differ by the factor h/r0. The phase and the group velocities for the finite element case are given by(41)vphFE=2Kmr0hκsin⁡κh2,vgrFE=K2mr0cos⁡κh2,since we assume that k=κ>0.Moreover, from (41) we can see that for fixed k the speed of a wave decreases as the mesh size h increases. We observe reflections at the interface (cf. (41)).In Figure5 the dependence of the velocity on the space discretization parameter and the wave number is shown. Here the velocity is the normalized phase velocity(42)v¯≔vphr02Km-1.One can clearly see that for h≪r0 the wave becomes slower.Figure 5
The different dependencies between the space parameters, the normalized velocity, and the wave number.Besides the size of the elements, the shape of its basis functions also plays a crucial rule. To see this, let us consider the piecewise linear basis functions.In molecular dynamics the spectrum of the displacement is large in the sense that there are low and high frequent parts. This fact can cause trouble when applying finite element basis functions based on polynomials. To see this, we observe that the Fourier transformation of our basis function isɛ-band limited; that is, its Fourier transform u^ is less than ɛ outside some compact interval Σ. As a consequence high frequent parts of a molecular spectrum being outside of the interval Σ are only represented poorly. In the literature, methods can be found, which attempt to capture these high frequent parts by enrichment basis functions.In Figure6, we can see the Fourier transform of the linear finite element basis function (red line).Figure 6
The Fourier transform of an initial waveq and its projection.Obviously, waves with low wave number are better represented (close toκ=0) than high frequencies.For a given initial amplitudeq with high and low frequencies the Fourier spectrum is given in Figure 6 (blue, dashed). Moreover, the Fourier spectrum of the corresponding coarse displacement (finite element) is given (red). This problem has been treated in the existing literature differently. Either damping techniques have been applied such that no or damped high frequent waves enter the coarse scale [29, 30, 45] or constraints (Lagrange multiplier) in the handshake region hinder the high frequent waves to enter the coarse scale [21, 29] or the discretization on the coarse scale is enriched [48].


## 6.Conclusion

In this paper, we show that coupling molecular dynamics simulations with continuum mechanics is a challenging task. We describe the respective equations for the two scales and showed that it is not clear at all how to combine them in a sound mathematical way. The spurious reflections which have been observed by several AtC articles are described and their origin is explained: Firstly, the reflections occur due to the different physical models; that is, the dispersion relation, being linear in continuum mechancis and not linear in the atomistic model. Secondly, the different length scale of the two discretizations, finite element mesh size on the one hand and the distance of the atoms on the other hand, causes a reflection which can be explained by Shannon’s sampling theorem. This paper is directed to scientists working for the first time on atomistic-to-continuum methods and it intends to explain the main challenges.



---

*Source: 834517-2015-06-28.xml*