Challenges in Atomistic-to-Continuum Coupling

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.


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][2][3][4][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.

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 Ω  ⊂ Ω 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.

Mathematical Problems in Engineering
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.

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 Ω ⊂ R  is decomposed by 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 −1 (ii) Handshake Region Methods: The atomistic and the continuum part are matched in an overlap region Ξ = Ω MD ∩ Ω CM , Ξ ⊂ R  , and meas  (Ξ) > 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 ).
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.

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.

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, In the engineering literature, these methods are also sometimes termed as domain-decomposition techniques.In general we expect, for an efficient multiscale or atomistic-tocontinuum 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][14][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][19][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][22][23][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][34][35].Recently, coupling methods have been proposed, where the microscale is a canonical ensemble, that is, AtC with finite temperature [24,[36][37][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.

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.

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 point   ∈ R  ,  ∈ 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 where   is the displacement of atom .
The atomic displacements  = (  ) ∈A are assumed to obey Newton's law of motion where  internal and  external are the internal and external forces.
With each atom , we associate the mass   > 0, such that is the mass matrix on the microscale.If the system is conservative a potential  is given; the internal forces acting on a conservative system can be obtained as  internal = −∇ X( X).

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 Ω ⊂ R  .In classical elasticity it is assumed that the motion, due to external and internal forces, is given by the mapping such that the deformed configuration at time  of the body is given by with the displacement () = (, ) − ,  ∈ Ω,  ≥ 0.
In the forthcoming we denote points of the deformed configuration by and write  instead of () and  instead of () whenever possible.Furthermore, it is assumed that  is locally injective in Ω, sufficiently smooth, and orientation preserving, which can be expressed by the pointwise inequality det (     ) (, ) > 0,  ∈ Ω,  ≥ 0, ,  = 1, . . . , 3. ( 9) Let us denote the deformation gradient by ; that is, The mass is given by where (, ) 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 The change of square length of the line segment  =  due to the deformation is given by the Green strain tensor:  = (1/2)(   − ) = (1/2)(∇  + ∇ + ∇  ∇).Moreover, the Cauchy stress is defined by () =  = , where  is the normal and  is the traction on any surface segment .For computational and analytical reasons, we use the Lagrangian description and thus define the first Piola-Kirchhoff tensor P, which is a second-order tensor, given by where  ref 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 where V = u = /.Here,  is the body force density and V ∈  1 ([0, ],  0 (Ω)) and P ∈  0 ([0, ],  1 (Ω)).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  may not be smooth enough to fulfill (13).Thus, we introduce the weak formalism.To do so, let Ω ⊂ R  be measurable.We denote by   (Ω) the space of -Lebesgue-integrable functions on Ω.Moreover, we denote the  2 scalar product by (⋅, ⋅)  2 (Ω) and the norm Let  be a multi-index with ‖‖ 1 := ∑  =1 |  |.Then, for a bounded domain Ω and for  ∈ N, 1 ≤  < ∞, we denote by  , (Ω) the Sobolev space, given by  , (Ω) := { ∈   (Ω) | there exists   ,    ∈   (Ω) ∀ with ‖‖ ≤ } (15) and equipped with the norm We assume that Ω is a polygonal domain and Ω = Γ  ∪ Γ  is the disjoint union of the portions Γ  and Γ  .Moreover, we assume that  ∈  0 ([0, ],  2 (Ω)),  ∈  0 ([0, ],  2 (Γ  )),   ∈  0 ([0, ],  1/2 (Γ  )), and the initial displacements  0 ∈  1 (Ω).Then for the weak formulation of (13) we seek  ∈  2 ([0, ],  1 (Ω)) such that (⋅, 0) =  0 , | Γ  =   , and for all V ∈  0 ([0, ];  1 (Ω)) vanishing on Γ  , where  :  denotes tr(  ).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  2 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  2 (Ξ).

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.
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 the 1 case.

Traveling Waves in Crystals.
Recalling that (, ) denotes the deviation of the atoms from their starting configuration a harmonic wave in the atomistic model is of the form where q is the amplitude,  −  is the phase,  is the wave vector, and  is the angular frequency.We also define the wave number  = || and the wave length  = 1/.In the forthcoming, we consider two different kinds of velocities.
The phase velocity V ph is defined as 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 the group velocity, defined as We can decide between longitudinal and transverse vibrations.An example of a 1 chain in 2 is given by Figure 4.
Remark 1.In 3 the wave vector  has three components and points into the direction of propagation.The atomistic displacements associated with a wave (, ) where  is the equilibrium position of an atom are simultaneously threedimensional vectors.These displacement vectors may be parallel to  (longitudinal), perpendicular to  (transversal), or along a direction, that is, not directly related to the direction of .
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.In 1 each  can be rewritten as  0 , where  ∈ Z and  0 is the atomistic spacing.Then, for each , we have the superposition and we can understand the behavior of the solution by examining the harmonic solutions independently; that is,   = 0 for all but one .
Let us consider a mass spring system with lattice spacing  0 , mass , and spring constant .Assuming that for this harmonic system a Hamiltonian is given, the respective equations of motion of H MD can be stated as which then gives us the dispersion relation for the molecular case (cf.[43,44]) One can see easily that () = (+2/ 0 ) for any  ∈ Z.We thus only consider the case of  ∈ (−/ 0 , / 0 ), which is known as the first Brillouin zone.Note that the open interval is chosen, as for  = / 0 the group velocity V gr 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  =  > 0.

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 sin ( *  0 ) =  *  0 +  ( *  0 ) ; then for all 0 <  <  * we have Moreover, the phase velocity V CM ph is then given by 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 , such that the average distance between the atoms becomes   =  0 (1 − ), where  ≪ 1.
With respect to the harmonic energy, given by  harm,1 = (1/2)  (  −  + ) 2 , the energy of the strained chain becomes Thus, the strain energy, which can be considered as the extra energy per atom, is given by where K :=  2 0 is the elastic constant.From the velocity given by (26) we have Summing up, for the long wavelength limit ( = 1/,  ≪ 1), we obtain The meaning of the different dispersion relations ( 23) and (30) becomes clear when the phase velocity V ph and the group velocity V gr are considered.It can easily be seen that the molecular phase and the group velocity are given by whereas their continuum counterparts are given by Thus, for small  we have In contrast, for large wave numbers, that is,  ≈ /2, we have cos( 0 /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.

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 field , we now employ a finite element discretization of lower order.Let T ℎ denote a mesh with mesh size parameter ℎ > 0, such that the family {T ℎ } ℎ is shape regular.
By using Lagrangian conforming finite elements of first order ( 1 ) for the displacements  and denoting the set of all nodes of T ℎ by N ℎ , the finite element space V ℎ ⊂  1 (Ω) is spanned by the nodal basis The Lagrangian basis functions  ℎ  ∈ V ℎ are uniquely characterized by the Kronecker-delta property where   is the Kronecker-delta.Any function  ℎ ∈ V ℎ (Ω) can uniquely be written as where (  ) ∈N ℎ ∈ R ⋅|N ℎ | ,   ∈ R  , is the coefficient vector.We can identify each element of V ℎ with its coefficient vector (  ) ∈N ℎ .In the forthcoming, we omit the superscript ℎ whenever possible.In space, (17) gives a nonlinear ordinary differential equation in R |N ℎ | 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 in 1.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 Hence, if the mesh size ℎ equals the atomistic spacing  0 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 is 1/ℎ and the maximum wave length is given by  such that Shannon's theorem becomes This shows that the mesh size ℎ 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 where ũ 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 Comparing (23) and ( 40) one can see that both relations only differ by the factor ℎ/ 0 .The phase and the group velocities for the finite element case are given by ) , since we assume that  =  > 0.Moreover, from (41) we can see that for fixed  the speed of a wave decreases as the mesh size ℎ increases.We observe reflections at the interface (cf. ( 41)).
In Figure 5 the dependence of the velocity on the space discretization parameter and the wave number is shown.Here the velocity is the normalized phase velocity One can clearly see that for ℎ ≪  0 the wave becomes slower.
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 û 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 Figure 6, we can see the Fourier transform of the linear finite element basis function (red line).
Obviously, waves with low wave number are better represented (close to  = 0) than high frequencies.
For a given initial amplitude  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].

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.

Figure 2 :
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).

Figure 3 :Figure 4 :
Figure 3: (a) Wave package in 1 in initial position.(b) Propagated initial wave: the waves with lower frequencies move faster.

Figure 5 :
Figure 5: The different dependencies between the space parameters, the normalized velocity, and the wave number.

Figure 6 :
Figure 6: The Fourier transform of an initial wave  and its projection.