Eduard Marenić

Jurica Sorić

Adnan Ibrahimbegović

Adaptive modelling in atomistic-to-continuum multiscale methods

Keywords: molecular mechanics, atomistic-to-continuum coupling, quasicontinuum, bridging domain, Arlequin, Cauchy-Born rule, RVE, goal error estimates, goals algorithm

Due to the lack of computational power to perform a fully atomistic simulation of practical, engineering systems, a number of concurrent multiscale methods is developed to limit atomic model to a small cluster of atoms near the hot spot. In this paper the overview of salient features of the main multiscale families is given. The special attention is drawn towards the role of model adaptivity, that is, which part of the problem domain to model by the atomic scale (the hot spot) and which by coarse scale model, as well as where to place the interface of the two models to control the accuracy. Taking Quasicontinuum method as a reference, review of the evolution of the Bridging domain/Arlequin method is given, which parallels the development of a posteriori modeling error estimation.

Introduction

The emphasis of scientific research in material science has shifted from micro-and mesoscale to the study of the behavior of materials at the atomic scale of matter. The first trends of this kind go back to the 1980s. when the scientists and engineers began to include atomistic descriptions into models of materials failure and plasticity [START_REF] Buehler | Atomistic Modeling of Materials Failure[END_REF]. This research is related to the terms nano-technology and nano-mechanics. At nano-scale, the effects of single atoms, individual molecule, or nano-structural features may dominate the material behavior, especially at failure [START_REF] Buehler | Atomistic Modeling of Materials Failure[END_REF]. The classical continuum mechanics, that has been the basis for most theoretical and computational tools in engineering [START_REF] Ibrahimbegovic | Nonlinear Solid Mechanics[END_REF], is not suitable for nano-scale applications. Thus, different kind of computational modeling, in particular atomistic and/or molecular simulation, has become increasingly important in the development of such new technologies [START_REF] Cleland | Foundations of Nanomechanics From Solid-State Theory to Device Applications[END_REF][START_REF] Phillips | Crystals, Defects and Microstructures Modeling Across Scales[END_REF][START_REF] Rapaport | The Art of Molecular Dynmics Simulations[END_REF].

For many engineering application domains, the numerical simulations of this kind are replacing the expensive experimental testing or are being used to complete the experimental observations and to increase the reliability of parameter identification in experiments under heterogeneous stress field. In the case of nano-mechanics it is usually impossible to perform the simple tests (such as the simple tension test), or most of experiments are very expensive and not very reliable. The experimental analysis of nano-mechanical properties at submicrometer scales de facto became possible with the developments of techniques relying upon the atomic force microscope (AFM), nanoindentation, or optical tweezers. These techniques and instrumentation can observe and characterize forces of the order of hundreds of pN, with displacements of the order of nanometers [START_REF] Buehler | Atomistic Modeling of Materials Failure[END_REF]. Atomistic simulation has been used in the research topics like: the atomic-scale effect in fracture and wear, dislocation dynamics in nano-indentation, nano-composites, carbon nano tubes, nano electro mechanical system components, semiconductors and biomechanics [START_REF] Buehler | Atomistic and continuum modeling of mechanical properties of collagen: Elasticity, fracture, and self-assembly[END_REF][START_REF] Buehler | Nanomechanics of collagen fibrils under varying cross-link densities: Atomistic and continuum studies[END_REF].

The main challenge is that atomistic models typically contain extremely large number of particles, even though the actual physical dimension may be quite small. For example, even a crystal with dimensions below a few micrometers sidelength has several tens of billions of atoms. Predicting the behavior of such large particle systems under explicit consideration of the trajectory of each particle is only possible by numerical simulation, and must typically involve the usage of the supercomputers [START_REF] Buehler | Atomistic Modeling of Materials Failure[END_REF]. Even though nanoscale systems and processes are becoming more viable for engineering applications, our ability to model their performance remains limited, since the fully atomistic simulations remain out of reach for engineering systems of practical interest.

Multiscale (MS) modeling methods have recently emerged as the tool of choice to link the mechanical behavior of materials from the smallest scale of atoms to the largest scale of structures [START_REF] Wing | Nano Mechanics and Materials Theory, Multiscale Methods and Applications[END_REF]. MS methods are often classified as either hierarchical or concurrent. Hierarchical methods are the most widely used, for their computationally efficiently. In these methods, the response of a representative volume element (RVE) at the fine scale is first computed, and from this a stress-strain law is extracted. Thus, the computations are performed on each scale separately and the scale coupling is often done by transferring the problem parameters leading to the classical problem of homogenization (e.g. see early work [START_REF] Sanchez-Palencia | Non-homogeneous media and vibration theory[END_REF]). For severely nonlinear problems, hierarchical models become more problematical, particularly if the fine scale response is path dependent. It should be noted that when failure occurs, in many circumstances hierarchical models are invalid and cannot be used [START_REF] Fish | Multiscale Methods Bridging the Scales in Science and Engineering[END_REF].

Concurrent methods, on the other hand, are those in which the fine scale model (e.g. atomistic, treated with molecular mechanics) is embedded in the coarse scale model (usually continuum model treated with FEM) and is directly coupled to it. In the study of fracture, for example, fine scale models can be inserted in hot spots where stresses become large and where there is the biggest risk of failure. These hot spots can be identified on the fly or by a previous run. Molecular mechanics (MM) and/or quantum mechanics (QM) models are required for phenomena such as bond breaking, but the relevant configuration is far too large to permit a completely atomistic description. In order to make such problems computationally tractable, the molecular model must be limited to small clusters of atoms in the vicinity of a domain of interest where such high resolution models are necessary and a continuum method should be used for the rest of the domain [START_REF] Fish | Multiscale Methods Bridging the Scales in Science and Engineering[END_REF]. Here we primarily focus on the concurrent, static (equilibrium), atomistic-to-continuum MS modeling, strongly coupling atomistic and continuum scales.

An overview of current research activities on MS methods can be found in several reviews [START_REF] Broughton | Concurrent coupling of length scales: Methodology and application[END_REF][START_REF] Curtin | Atomistic/continuum coupling in computational materials science[END_REF][START_REF] Kam | An introduction and tutorial on multiple-scale analysis in solids[END_REF][START_REF] Liu | An introduction to computational nanomechanics and materials[END_REF][START_REF] Ronald | A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods[END_REF][START_REF] Rudd | Concurrent coupling of length scales in solid state systems[END_REF][START_REF] Srivastava | Computational nanotechnology: A current perspective[END_REF], each giving a preference to a preferred choice of the method or its particular feature. For that reason, we seek to give a more complete overview covering all the salient features of the main families, each covered with a brief, but pointed discussion. There is also a novel idea to draw attention towards a special role of adaptivity in providing an optimal form of the atomistic-to-continuum coupling based on the overlapping domain decomposition. With such focus on adaptivity, this paper compares quasicontinuum (QC) method and the bridging domain (BD) or Arlequin based coupling. The QC method uses an adaptive coarse graining approach rather than classical coupling, it can be used as a reference for adaptive strategy. More precisely, the question we thus address pertains to which part of the domain should be modeled by the fine (atomic) scale and which by coarse scale model in a particular problem, and where to place the interface of the two models to control the accuracy for any solution stage?

The paper is organized as follows. Following this introduction, the paper starts with the short overview of atomic models in Section 2 and finishes with the motivation for the MS methodology and a short list of leading MS methods. Then the standard QC and BD/Arlequin methods are described in Sections 3 and 4, respectively. Following the general description of the two methods the comparison is summarized in 5. Numerical examples of model adaptivity are presented in Section 6 in order to illustrate these ideas, and provide quantitative comparison between different adaptive modelling strategies. The concluding remarks are stated in Section 7.

2 Atomistic (particle) model

Atomistic interaction modeling

The first choice that should be made for any kind of material modeling is the energy function describing the system of interest. Once the energy of the atomic interaction is defined, the essential of material behavior is determined. The main goal of this section is to give a general introduction to atomistic i.e. non-continuum material modeling, to introduce the essential ideas and review the literature. It is worth to mention that the reason for introducing the atomistic modeling is two-folded. First, the atomistic modeling is used as a testing ground for energetics of the system, by using the simplest generic form of the interaction model. Second, the potential function (U ) driving the molecular system can take an extremely complicated form, when the goal is to represent the quantitative predictions for specific material. In this case, an accurate representation of the atomic interactions has to be material specific. The nature of these interactions is due to complicated quantum effects taking place at the subatomic level that are responsible for chemical properties such as valence and bond energy [START_REF] Ercolessi | A molecular dynamics primer[END_REF][START_REF] Griebel | Numerical Simulation in Molecular Dynamics[END_REF][START_REF] Wing | Nano Mechanics and Materials Theory, Multiscale Methods and Applications[END_REF]. However, quantum mechanics based description of atomic interaction is not discussed herein, emphasis is rather on the empirical interaction models that can be derived as the result of such computations. Alternatively, the function U in classical interatomic potential that can be obtained from experimental observations and should accurately account for the quantum effects in the average sense. However, many different expressions can be fit to closely reproduce the energy predicted from quantum mechanics methods (semi-empirical), while retaining computational efficiency [START_REF] Allen | Computer simulation of liquids[END_REF][START_REF] Buehler | Atomistic Modeling of Materials Failure[END_REF]. Needless to say, there is no single approach that is suitable for all materials and for all different phenomena of material behavior that we need to describe. The choice of the interatomic potential depends very strongly on both the particular application and the material.

The general structure of the potential energy function for a system of N atoms is

U (r 1 , r 2 , . . . , r N ) = i V 1 (r i ) + i,j>i V 2 (r i , r j ) + i,j>i,k>i V 3 (r i , r j , r k ) + • • • . ( 1
)
where the function V m , ∀m = 1, 2, . . ., represents the m-body potential and r i the position vector of the atom i. The first term of the equation (1) indicates the effect of an external force field on the system where it is immersed, such as gravitational or electrostatic. This term is usually ignored in practice, [START_REF] Wing | Nano Mechanics and Materials Theory, Multiscale Methods and Applications[END_REF]. The second term V 2 or V ij shows pair-wise interaction depending only on the pair separation r ij = |r i -r j | between atoms i and j. The threebody term involves angle-dependent forces, whereas four-body term includes torsion effects.

In short, m-body potential terms for m > 2 are usually called multi-body potentials. The simplest form used for practical reasons is when the sum in ( 1) is truncated after second term resulting with the pair-wise potential.

Pair-wise potentials

The total energy of the system in pair potentials is given by summing the energy of all atomic bonds 1 (V 2 (r ij )) over all N particles in the system.

U = 1 2 N i =j=1 N j=1 V 2 (r ij ) = N -1 i=1 N j>i V 2 (r ij ). ( 2
)
Note the factor 1/2 accounts for the double counting of atomic bonds. One of the most well known interatomic potentials is the Lennard-Jones (LJ), or yet called 6-12 potential. The potential energy function for the LJ potential is expressed as

V LJ 2 (r ij ) = 4 σ r ij 12 - σ r ij 6 , (3)
where and σ are constants chosen to fit material properties (no relation to continuum stress and strain, see Fig. 1) and r ij is the distance between two atoms i and j. The 1/r 12 ij term is meant to model the repulsion between atoms as they approach each other, and is motivated by the Pauli principle in chemistry. The Pauli principle implies that as the electron clouds of the atoms begin to overlap, the system energy increases dramatically because two interacting electrons cannot occupy the same quantum state. The 1/r 6 ij term adds cohesion to the system, and is meant to mimic van der Waals type forces. The van der Waals interactions are fairly weak in comparison to the repulsion term, hence the lower order exponential is assigned to the term. LJ 6-12 is an example of potential limited to the simulations where a general class of effects is studied, instead of specific physical properties, and a physically reasonable yet simple potential energy function is desired [START_REF] Kam | An introduction and tutorial on multiple-scale analysis in solids[END_REF].

Since the LJ potential is highly nonlinear function of the atom pair distance r ij , it is sometimes useful to use its linearized form in terms of so-called harmonic potential

V H 2 (r ij ) = 1 2 k ij (r ij -r ij,0 ) 2 , (4)
where r ij,0 is the initial (equilibrium) atomic pair distance, and k ij is the bond stiffness. The harmonic potential can describe the atomic system behavior for small atomistic separation (refer Fig. 1). Hence, this potential is usually chosen as the first and simplest description of the atomic interaction, in particular in development of the MS methods where the emphasis is on the coupling and not on the accurate and realistic description of different material mechanisms.

LJ look alike potentials are the Morse and Buckingham potentials. The Morse potential consists of the exponential repulsion and attraction and three adjustable parameters [START_REF] Sunyk | On Aspects of Mixed Continuum-Atomistic Material Modelling[END_REF]. It is originally designed for covalent bond which is strongly space oriented and a description of radial stretching is not sufficient to describe it. The Morse potential is computationally more expensive than the LJ potential due to the exponential term but it models interaction in a more realistic way. The Buckingham potential consists of more physical exponential Born-Meyer repulsion and the van der Waals attraction but at the small inter-atomic separations the potential becomes un-physical (often referred to as Buckingham catastrophe) [START_REF] Sunyk | On Aspects of Mixed Continuum-Atomistic Material Modelling[END_REF]. In the sequel the example of the improved pairwise potential is described. 1 Often the total energy of the system of atoms is represented as the sum of atom energies. That is the energy is given on per atom basis (not per bond) as U = 1 2 N i=1 E i , where E i is the energy of atom i. Lennard-Jones and Harmonic potential (dashed line). Note that the Harmonic potential is a suitable approximation when the particles are around the equilibrium position.

Embedded atom method (EAM) potentials and multi-body potentials

This kind of interaction models are widely used to model metals [START_REF] Daw | Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals[END_REF][START_REF] Doyama | Embedded atom potentials in fcc and bcc metals[END_REF]. A local density dependent contribution is added to the pair potential energy function (2) by embedding energy term

F U = i,j =i V 2 (r ij ) + i F (ρ i ). (5)
The embedding energy F (ρ i ) is related to the environment of the atom i, and ρ i is the local electron density. The main advantage of this kind of potential is the ability to describe surfaces or cracks since they incorporate the variation of bond strength with coordination (density).

A number of different approaches to realistic description of behavior of solids based on the first-principles or quantum mechanical calculations have been developed in recent decades. These models account for the environmental dependence of the bond strength with the force between any two particles depending upon the position of several neighbouring particles. Besides EAM, the other examples of this kind of potentials are Finnis and Sinclair potential describing metallic bonds, and the Brenner potential used for hydrocarbon bonds [START_REF] Griebel | Numerical Simulation in Molecular Dynamics[END_REF].

Solution strategy and motivation for MS methods

We focus in this work upon the mechanics behavior only, in the context of quasi-static loading applications. The equilibrium configuration of solids corresponds to a state of minimum energy. Similar to the FEM, that the positions of all nodes are determined by minimizing the energy in the solid. Thus, for a system of N atoms, the equilibrium configuration is determined by minimizing

Π = U - N i fi • d i , (6)
where U denotes the energy stored in the atomic bonds, and d i and fi denotes displacement and external force on atom i, respectively. Since the continuum FEM and molecular mechanics share a common ground of energy minimization there is a number of contributions regarding so called atomistic FE approach (or AFEM method), inserting molecular mechanics in the context of FEM, see e.g. [39-41, 68, 69].

The task of minimizing [START_REF] Paul | On the application of the arlequin method to the coupling of particle and continuum models[END_REF] quickly becomes intractable for large number of particles (atoms). This, together with the assumption that the calculation of specific quantities of the solution can be accurately approximated by replacing the particle model by a coarser model (i.e. continuum model), is the basis for multiscale (MS) modeling. Extensive work has been done in the development of atomistic-to-continuum MS modelling approaches, starting with early works by [START_REF] Mullins | Simulation of the (001) plane crack in alpha-iron employing a new boundary scheme[END_REF] [START_REF] Mullins | Simulation of the (001) plane crack in alpha-iron employing a new boundary scheme[END_REF] and Kohlhoff et al. (1990.) [START_REF] Kohlhoff | Crack propagation in b.c.c. crystals studied with a combined finite-element and atomistic model[END_REF]. Mullins simulated 2D cracks in α-iron with the atomic scale models, and due to the restrictions of the computational power the question was how to connect the atomic model and surrounding continuum. Kohlhoff et al. proposed a new method for combined FE and atomistic analysis of crystal defects, called FEAt. Both papers dealt with the problem of proper treatment of the transition between the lattice and continuum. A number of MS methods developed recently from theoretical standpoint of view appear very different. However, as shown in [START_REF] Ronald | A unified framework and performance benchmark of fourteen multiscale atomistic/continuum coupling methods[END_REF], at the implementation level all these methods are in fact very similar. The performance of a number of most frequently used methods is compared in a linear framework on a common benchmark test. The ways in which various multiscale methods differ are the formulation (energy or force based model), the coupling boundary conditions, the existence of the handshake region (overlapped or non-overlapped domain decomposition), and the choice of the continuum model. The list containing: quasicontinuum (QC) method (in Section 3), bridging domain (BD) method (in Section 4), coupling of length scales (CLS), bridging scale (BS) method [START_REF] Kam | An introduction and tutorial on multiple-scale analysis in solids[END_REF][START_REF] Karpov | Multiscale boundary conditions in crystalline solids: Theory and application to nanoindentation[END_REF][START_REF] Qian | A multiscale projection method for the analysis of carbon nanotubes[END_REF], coupled atomistics and discrete dislocations (CADD), Atomistic-to-continuum coupling (AtC) [START_REF] Badia | On atomistic-to-continuum coupling by blending[END_REF][START_REF] Badia | Nuggehally, and Max Gunzburger. A force-based blending model foratomistic-to-continuum coupling[END_REF][START_REF] Fish | Concurrent AtC coupling based on a blend of the continuum stress and the atomistic force[END_REF], etc. is also not exhaustive but the unified framework, available computer code, and a quantitative comparison between the methods offer a good overview. Note that there is also a very recent effort of coupling non-local to local continuum [START_REF] Han | Coupling of nonlocal and local continuum models by the arlequin approach[END_REF] in the Arlequin framework (see Section 4). An alternative to discrete modeling of atomic/particle systems is the use of non-local continuum mechanics models (NLCM) [START_REF] Lubineau | A morphing strategy to couple non-local to local continuum mechanics[END_REF]. NLCM reduces the computational costs but has the ability to capture non-local interactions. However since the simulation using NLCM is also costly due to assembly operation of the discretized model where each interaction point interacts with multiple neighbours, and the fact that this reduces the sparsity of the matrices, similar principle as in BD and QC method of coupling discrete, non-local particle model with local continuum is used. The key challenge is then again the gluing of non-local continuum model with the local one.

The standard approach in these models is to a priori identify the atomistic and continuum regions and tie them together with some appropriate boundary conditions. In addition to the disadvantage of introducing artificial numerical interfaces into the problem a further drawback of many of these models is their inability to adapt to changes in loading and an evolving state of deformation. Take for example the problem of nanoindentation. As the loading progresses and dislocations are emitted under the indenter the computational model must be able to adapt and change in accordance with these new circumstances.

In the sequel the QC and the BD/Arlequin methods are described in more detail. The goal is, however, to show the evolution of the BD/Arlequin coupling approach and to compare the features regarding ability to adapt.

Quasicontinuum method

The Quasicontinuum (QC) method is originally proposed in late 90's by Tadmor, Ortiz and Phillips [START_REF] Tadmor | Quasicontinuum analysis of defects in solids[END_REF]. Since then it has seen a great deal of development and application by a number of researchers. The QC method has been used to study a variety of fundamental aspects of deformation in crystalline solids, including fracture [START_REF] Hai | Deformation twinning at aluminum crack tips[END_REF][START_REF] Miller | Quasicontinuum models of fracture and plasticity[END_REF][START_REF] Miller | Quasicontinuum simulation of fracture at the atomic scale[END_REF] 2 , grain boundary slip and deformation [START_REF] Shenoy | An adaptive finite element approach to atomic-scale mechanics-the quasicontinuum method[END_REF]. The nano-indentation [START_REF] Vivek | Nucleation of dislocations beneath a plane strain indenter[END_REF] and similar applications are examples where neither atomistic simulation nor continuum mechanics alone were appropriate, whereas the QC was able to effectively combine the advantages of both models. The main goal of the QC method is the provide a seamless link of the atomistic and continuum scales. This goal is achieved by the three main building blocks [START_REF] Eidel | Atomistic simulation methods and their application on fracture[END_REF][START_REF] Miller | The quasicontinuum method: Overview, applications and current directions[END_REF]:

1. Reduction of degrees of freedom (DOF) by coarse-graining of fully atomistic resolution via kinematic constraints. The fully atomistic description is retained only in the regions of interest.

2. An approximation of the energy in the coarse grained region via numerical quadrature.

The main idea is to avoid the need to calculate the energy of all the atoms, but retain only a few so-called rep-atoms.

3. Ability of the fully refined, atomistic region to evolve with deformation, where adaptivity is directed by suitable refinement indicator.

DOF reduction or coarse graining

If the deformation changes gradually on the atomistic scale, it is not necessary to explicitly track the displacement of every atom in the region. Instead it is sufficient to consider some selected atoms, often called representative atoms or rep-atoms. This process is in essence the upscaling via coarse graining. Only rep-atoms have independent DOF while all other atoms are forced to follow the interpolated motion of the rep-atoms. The QC incorporates such a scheme by means of the interpolation functions of the FE method, and thus the FE triangulation has to be performed with rep-atoms as FE mesh nodes. This way continuum assumption is implicitly introduced in QC method. Thus, if the potential Π is given as a function of displacement u (similarly as in ( 6))

Π(u) = E tot (u) - N i=1 fi u i , ( 7
)
where fi is the external force on the atom i and E tot is an atomistic internal energy

E tot = N i=1 E i (u), (8)
the kinematic constraint described above is performed by replacing E tot with E tot,h

E tot,h = N i=1 E i (u h ). ( 9
)
In the above equation the displacement approximation is given via standard FE interpolation

u h = Nrep i=1 N i (u i ), (10)
where N i is a shape function for the node/rep-atom i. The density of rep-atoms vary in space according to the considered problem. In the vicinity of region of interest every atom is considered as rep-atom and in region of more slowly varying deformation gradient, only a few atoms are chosen.

Efficient energy calculation via Cauchy-Born rule

Described kinematic constraint on most of the atoms in the body will achieve the goal of reducing the number of degrees of freedom in the problem. However, for the purpose of energy minimization the energy of all the atoms (not just rep-atoms) has to be computed. The way to avoid visiting every atom is the Cauchy-Born (CB) rule [START_REF] Ericksen | The cauchy and born hypotheses for crystals. Phase transformation and material instabilities in solids -from book[END_REF][START_REF] Ericksen | On the cauchyborn rule[END_REF][START_REF] Zanzotto | The cauchy-born hypothesis, nonlinear elasticity and mechanical twining in crystals[END_REF]. The CB rule postulates that when a simple, mono-atomic crystal is subjected to small displacement on its boundary then all the atoms will follow this displacement. In QC this rule is implemented in that every atom in a region subject to a uniform deformation gradient will be energetically equivalent. Thus, energy within an element can be estimated by computing the energy of one, single atom in the deformed state. The estimation is performed simply by multiplying the single atom energy by the number of atoms in the specific element. The strain energy density (SED) of the element can be expressed as:

W (F) = E 0 (F) Ω 0 , (11)
where E 0 is s the energy of the unit cell when its lattice vectors are distorted according to deformation gradient F, and Ω 0 is the volume of unit cell. The sum in eq. ( 9) is reduced to number of FEs (N elem )

E tot,h ≈ E tot,h = N elem e=1 Ω e W (F e ), (12)
where the element volume and unit cell volume are related as n e Ω 0 = Ω e , and n e is the number of atoms contained in element e. Using the CB rule, the QC can be thought of as a purely continuum formulation (local QC), but with a constitutive law that is based on atomistic model rather than on an assumed phenomenological form [START_REF] Miller | The quasicontinuum method: Overview, applications and current directions[END_REF]. Within QC framework, the calculation of CB energy is done separately in a subroutine. For a given deformation gradient F the lattice vectors in a unit cell are deformed according to given F and the SED is obtained according to eq. ( 11).

Non-local QC and local/non-local coupling

In settings where the deformation is varying slowly and the FE size is adequate with respect to the variations of the deformation, the local QC is sufficiently accurate and very effective.

In the non-local regions, which can be eventually refined to fully atomistic resolution, the energy in ( 9) can be calculated by explicitly computing only the energy of the rep-atoms by numerical quadrature

E tot,h ≈ E tot,h = Nrep i=1 n i E i (u h ) ( 13
)
where n i is the weight function for rep-atom i and is high for rep-atoms in regions of low rep-atom density and low for high density. Thus, n i is the number of the atoms represented by the i-th rep-atom with the limiting case of n i = 1 for fully atomistic case and consistency requirement

Nrep i=1 n i = N. (14)
The main advantage of the non-local QC is that when it is refined down to the atomic scale, it reduces exactly to lattice statics. High accuracy of non-local formulation can be combined with the high efficiency of the local formulation. In order to do that non-local formulation is employed in the region where atomic scale accuracy is needed, and local where the deformation is changing relatively slow. In this coupling approach the rep-atom can be chosen as local or non-local depending on its deformation environment giving N rep = N loc + N nonloc . Total energy is approximated as

E tot,h = N nonloc i=1 n i E i (u h ) + N loc i=1 n i E i (u h ), (15)
where the weights n i are determined from the Voronoi tessellation i.e. by means of the cells around each rep-atom. The cell of atom i contains n i atoms, and of these atoms n e i reside in FE e adjacent to rep-atom i. The weighted energy contribution of rep-atom i is then found by applying the CB rule within each element adjacent to i such that

n i E i = N i el e n i Ω 0 W (F e ), n i = N i el e n e i , (16)
where Ω 0 is the cell volume for single atom, and N i el is the number of FE adjacent to atom i.

Local/non-local criterion

The criterion to trigger the non-local treatment is based on the significant variation of deformation gradient 3 . Precisely, we say that the state of deformation near a representative atom is nearly homogeneous if the deformation gradients that it senses from the different surrounding elements are nearly equal. The non-locality criterion is then:

max a,b,k |λ a k -λ b k | < , (17)
where λ a k is the k-th eigenvalue of the right stretch tensor for element a, k = 1 . . .

Adaptivity

Without a priori knowledge of where the deformation field will require fine-scale resolution, it is necessary that the method should have a built-in, automatic way to adapt the finite element mesh through the addition or removal of rep-atoms. This is a feature that is in QC inherent from the FE literature, where considerable attention has been given to adaptive meshing techniques for many years. Typically in FE techniques, a scalar measure is defined to quantify the error introduced into the solution by the current density of nodes (or repatoms in the QC). Elements in which this error estimator is higher than some prescribed tolerance are targeted for adaptation, while at the same time the error estimator can be used to remove unnecessary nodes from the model.

The error estimator in terms of deformation gradient is defined as the difference between the actual solution and the estimate of the higher order solution (see [START_REF] Miller | The quasicontinuum method: Overview, applications and current directions[END_REF]). If this error is small, it implies that the higher order solution is well represented by the lower order elements in the region, and thus no refinement is required. Needles to say, elements for which the error is greater than some prescribed error tolerance are targeted for refinement. Refinement then proceeds by adding three new rep-atoms at the atomic sites closest to the mid-sides of the targeted elements 4 . If the nearest atomic sites to the mid-sides of the elements are the atoms at the element corners, the region is fully refined and no new rep-atoms can be added. The same error estimator is used in the QC to remove unnecessary rep-atoms from the mesh. In this process, a rep-atom is temporarily removed from the mesh and the surrounding region is locally re-meshed. If all of the elements produced by this re-meshing process have a value of the error estimator below the threshold, the rep-atom can be eliminated. Essentially, the idea is to examine the necessity of each node. To prevent excessive coarsening of the mesh far from defects the nodes corresponding to the initial mesh are usually protected from deletion [START_REF] Shenoy | An adaptive finite element approach to atomic-scale mechanics-the quasicontinuum method[END_REF].

Bridging domain and Arlequin-based coupling

The Bridging domain (BD) method is in essence a partially overlapping domain decomposition scheme used for atomistic-to-continuum coupling developed by Belytschko and Xiao in 2003 [START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF] for the static, and [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF] for dynamical problems. The compatibility in the overlapping domain is enforced by Lagrange multipliers. The evolution of the aforementioned method (see also [START_REF] Anciaux | Ghost force reduction and spectral analysis of the 1d bridging method[END_REF][START_REF] Belytschko | A continuum-to-atomistic bridging domain method for composite lattices[END_REF][START_REF] Zhang | A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids[END_REF]) has much in common with recent works in the finite element (FE) community on the coupling of nonconforming meshes in the overlapping subdomain. This approach is known as Arlequin method developed by Ben Dhia [START_REF] Hachmi Ben Dhia | Multimodeling of multi-alterated structures in the Arlequin framework. solution with a domain-decomposition solver[END_REF][START_REF] Hashmi | The Arlequin method as a flexible engineering design tool[END_REF][START_REF] Guidault | On the l2 and the h1 couplings for an overlapping domain decomposition method using lagrange multipliers[END_REF]. The same Arlequin approach is lately also aplied for atomistic-to-continuum coupling like in [START_REF] Paul | On the application of the arlequin method to the coupling of particle and continuum models[END_REF][START_REF] Paul | Adaptive multiscale modeling of polymeric materials with arlequin coupling and goals algorithms[END_REF][START_REF] Chamoin | Ghost forces and spurious effects in atomic-to-continuum coupling methods by the arlequin approach[END_REF][START_REF] Dhia | A new adaptive modeling strategy based on optimal control for atomic-to-continuum coupling simulations[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with l 2 or h 1 couplings[END_REF][START_REF] Prudhomme | Computational analysis of modeling error for the coupling of particle and continuum models by the Arlequin method[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF]. More precisely, the domain with models coupling Ω is divided in three subdomains as shown in Fig. 2. The atomistic domain Ω a is discretized with molecular dynamics or rather molecular mechanics (MM), whereas the continuum mechanics domain Ω c discretization is carried out by FEs. The atomistic and continuum domains overlap in the domain Ω b = Ω a ∩ Ω c . This overlapping domain is also called bridging, handshake or coupling domain. The role of the continuum model is to replace the molecular model with a coarser, and thus computationally cheaper, model in Ω c ⊂ Ω away from the region of interest (e.g. lattice defect). Initially emphasis of the research was to make the coupling of the two different models as seamless as possible. No special attention was devoted to the question how to adaptively refine the model around the region of interest and where to position the handshake zone i.e. how far from the region of interest.

Construction of surrogate model

The role of the surrogate model is to propagate only the large-scale information. The choice of this model depends on the nature of the material but it should be selected as the most 'compatible' model with the atomistic or particle model in the sense of homogenization [START_REF] Chamoin | Ghost forces and spurious effects in atomic-to-continuum coupling methods by the arlequin approach[END_REF]. Thus, the material parameters of the surrogate continuum constitutive model should be calibrated accordingly. To that end, there are two approaches that appear in the BD/Arlequin literature.

The first one is related to the construction of constitutive equation via the Cauchy-Born rule introduced by the quasi-continuum (QC) approach as proposed e.g. in [START_REF] Xiao | A bridging domain method for coupling continua with molecular dynamics[END_REF][START_REF] Zhang | A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids[END_REF]. The Cauchy-Born rule is described in detail in section 3.2. The second approach pertains to computing the equivalent continuum model parameters through homogenization. Simple illustration for 1D case is given in [START_REF] Paul | On the application of the arlequin method to the coupling of particle and continuum models[END_REF] for the case of linear elastic continuum. More systematic approach to calibrating the continuum model parameters exploits the virtual experiments on the representative volume element (RVE) as suggested in [START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF]. The continuum model is based on plane stress linear elasticity and the constitutive relation is defined by Hooke's law. A RVE is considered to be a piece of material (atomistic lattice in this case) whose dimension is increased in iteratively until the consistent homogenized medium is obtained where the material parameters do not vary with further size increase. The lattice samples are larger than the effective RVEs in order to avoid boundary effects. The choice of the continuum model is, naturally, problem dependent (see [START_REF] Paul | Adaptive multiscale modeling of polymeric materials with arlequin coupling and goals algorithms[END_REF] for nonlinear hyperelastic material model suitable for polymeric materials).

An example of the (virtual) experiment for the uniaxial tension case is described in subsequently. Let the lattice sample size be l 1 × l 2 with constraints on the left and bottom atom layers, and imposed force (to obtain traction t) or displacement ū on every right-most atom Fig. 3. Without loss of generality we take the pair-wise interaction, where the internal energy . . . of the RVE is calculated as in (2) with ∀i, j ∈ RV E. The strain energy density is calculated by dividing the atomic energy U RV E by the initial volume of the RVE (V i )

W a = U RVE V i . ( 18
)
The variation of the atomistic SED with the size of the RVE is shown in Fig. 4. It is assumed that the energy densities obtained from the continuum and atomistic models are identical, so that W = W a 5 , where W is continuum strain energy density. Such value of W is further used to calculate the continuum part of energy (as indicated in equations ( 19) and ( 20)). This kind of approach obviates the need for the CB hypothesis to be fulfilled, which is used for composite lattices e.g. in 1D setting [START_REF] Chamoin | Ghost forces and spurious effects in atomic-to-continuum coupling methods by the arlequin approach[END_REF].

Governing equations and coupling

Total potential energy of the system may be written as

Π = E tot,w = E c w (u) + E a w (d) -E ext (u, d), (19)
where u i d are displacement vectors in the continuum and atomistic domains, respectively. Furthermore, E ext is the work of external forces, while E c w and E a w are weighted continuum and atomistic energies, defined as

E c w (u) = Ω c w c (X)W (F) dΩ c , E a w (d) = 1 2 i,j∈Ω a w a V ij . (20)
In the bridging domain the two models overlap, and the weighting functions w c and w a in (20) partition the energy. The weighting function serves to blend the behavior from the continuum model (w c ) and the atomistic model (w a ) and to avoid the double counting of the energy in the bridging domain. Furthermore, the use of an overlapping subdomain obviates the need for the FE nodes of the continuum model to correspond to the atomic positions. The weighting functions w c and w a define a partition of unity of the energy in the bridging domain as follows:

w c (X) = 1 for X ∈ Ω c \ Ω b , w a (X) = 1 for X ∈ Ω a \ Ω b , (21)
w c (X) + w a (X) = 1 for X ∈ Ω b .
The energy weighting functions can be taken with constant value (e.g. 0.5), linear (ramp) and cubic functions of X in Ω b for 1D case see Fig. 5). We note that all three possible distributions for the atomistic weighting function w a are depicted in Fig. 5, but only the linear continuum weighting function for continuum energy is indicated in Fig. 5) by dashed line in order not to overburden that illustration. Already in early contributions by [START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF][START_REF] Zhang | A bridging domain and strain computation method for coupled atomistic-continuum modelling of solids[END_REF] the discrete coupling of the atomistic and continuum models was achieved by forcing displacement compatibility in the bridging domain as u(X = X i ) = d i , ∀i ∈ Ω b . There, the Lagrange multiplier (LM) method was used to convert the problem of constrained minimization into finding the minimum of the larger, unconstrained problem. Hence, LMs, denoted with λ, are used to enforce the compatibility between the discrete atomistic displacement (discretely for each atom within the coupling zone) and the continuum displacement field. This results with the following Lagrangian

W L = E tot,w + C = E tot,w + i∈Ω b Ω b λ(X) • [u(X) -d i ] δ(X -X i )dΩ, ( 22
)
where C is the constraint in terms of the energy. In numerical implementation, displacement field in Ω c and LM fields in Ω b are approximated by using, respectively, the shape functions N i (X) and N λ k (X) as

u(X) = i∈S N i (X)u i , λ(X) = k∈S λ N λ k (X) λk , (23)
with u i and λk as the corresponding nodal values. Two limit cases regarding the LM field approximation are: i) the strict (or so called non-interpolated or atomic/particle) coupling with the LMs defined with atoms in the bridging zone and where i.e. N λ k (X i ) = δ ki ; ii) the interpolated (or continuum) coupling where the λ-nodes coincide with FE nodes and the LM shape functions N λ k correspond to the FE shape functions N k . The distribution of the λ-nodes for the two cases is illustrated in Fig. 6 for 1D case. In recent works on BD/Arlequin method e.g. [START_REF] Paul | On the application of the arlequin method to the coupling of particle and continuum models[END_REF][START_REF] Chamoin | Ghost forces and spurious effects in atomic-to-continuum coupling methods by the arlequin approach[END_REF][START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with l 2 or h 1 couplings[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF][START_REF] Qiao | Implementation of the Arlequin method into ABAQUS: Basic formulations and applications[END_REF] the displacement compatibility is given as u(X) = d b (X), where d b (X) is the regularized atomistic displacement field in Ω b that can be interpolated. For example, the atomistic displacement field is based on an MLS approximation in [START_REF] Guidault | Bridging domain methods for coupled atomistic-continuum models with l 2 or h 1 couplings[END_REF], or on a linear polynomial interpolations in [START_REF] Paul | On the application of the arlequin method to the coupling of particle and continuum models[END_REF]. As in the discrete case, the Lagrange multiplier field λ is used to enforce the displacement continuity in a weak sense defined through the scalar product (λ, u cd b ). Two kind of scalar products are considered (λ, u cd b ) L 2 and (λ, u cd b ) H 1 together with the performance and applicability in atomistic-to-continuum coupling. The coupling term C (as shown in [START_REF] Eidel | Atomistic simulation methods and their application on fracture[END_REF]) is now defined as follows

C L 2 = Ω b λ • (u c -d b )dΩ, ( 24
)
C H 1 = Ω b λ • (u c -d b ) + l 2 ε(λ) : ε(u c -d b )dΩ, ( 25
)
where l has the unit of a length, and corresponds to a characteristic dimension of the coupling zone. Apart from the advances in the coupling itself which is mostly related to the development of the Arlequin method advocated in initial work by Ben Dhia [START_REF] Hashmi | The Arlequin method as a flexible engineering design tool[END_REF] and its further application to the atomistic-to-continuum coupling, this method is acquiring the ability to accommodate the model and decrease the error in chosen quantity of interest. That is, the adaptivity described above for the QC method was included in the BD/Arlequin. This evolution parallels recent development in goal oriented error estimate theory as discussed in forthcoming section.

Adaptivity and error estimate

In computer simulations of physical models there are two major sources of error. Approximation error due to the the discretization of mathematical models, and modeling error related to the simplification or in general to the natural imperfections in abstract models of actual physical phenomena. The focus is here on the estimation and control of modeling error. This subject has been introduced in recent years and was initially devoted to estimating global modeling error e.g. [START_REF] Ainsworth | A posteriori error estimation in finite element analysis[END_REF]. Since then, extensions to error estimates in specific quantities of interest (estimate upper and lower bounds of error in linear functionals) have been proposed [START_REF] Oden | Estimation of modeling error in computational mechanics[END_REF][START_REF] Oden | Estimation of local modeling error and goal-oriented adaptive modeling of heterogeneous materials: I. error estimates and adaptive algorithms[END_REF][START_REF] Prudhomme | Practical methods for a posteriori error estimation in engineering applications[END_REF]. As an example Oden and Vemaganti [START_REF] Oden | Estimation of local modeling error and goal-oriented adaptive modeling of heterogeneous materials: I. error estimates and adaptive algorithms[END_REF] proposed an extension of a posteriori modeling error estimation for heterogeneous materials to the quantities of interest so-called goal-oriented error estimates. Many candidates for local quantities of interest are de facto quantities that one actually measures when assessing mechanical response e.g. average stresses on material interfaces, displacement, etc. Mathematically, a quantity of interest is any feature of the fine-scale solution that can be characterized as a continuous linear functional on the space of functions to which the fine-scale solution belongs. Analogously as in [START_REF] Oden | Estimation of local modeling error and goal-oriented adaptive modeling of heterogeneous materials: I. error estimates and adaptive algorithms[END_REF], where the error estimates are related to the error between fine-scale and regularized (homogenized) model, goal-oriented error estimation is extended to the case of discrete models (lattice) in [START_REF] Oden | On the extension of goal-oriented error estimation and hierarchical modeling to discrete lattice models[END_REF]. That is, this approach is used to estimate the modeling error between the atomic structure (lattice) and the surrogate, continuum model (i.e. FE discretization of the continuum model).

The idea behind the goal oriented adaptive modeling (as shown in the mentioned references) is to start from coarse, regularized model and to adaptively proceed towards fine model. Hence, the model is adopted to deliver local quantity of interest to within preset accuracy. The general process of adapting the surrogate model in order to decrease the modeling error in specific quantities of interest is referred to as the Goals Algorithms.

The basis of goal oriented error estimates is furthermore employed in the coupling of atomic and continuum models. The difficulty in the use of such coupling methods is to decide where to locate the overlap region between the two models so as to control the accuracy of their solutions with respect to the fully atomistic model. The convergence study of the modeling error in the context of atomistic-to-continuum (Arlequin type) coupling approaches is performed in [START_REF] Dhia | A new adaptive modeling strategy based on optimal control for atomic-to-continuum coupling simulations[END_REF][START_REF] Prudhomme | Computational analysis of modeling error for the coupling of particle and continuum models by the Arlequin method[END_REF]. The study is performed on a simple 1D problem that consists of chain of springs (as a fully atomic model), with a local defect modeled by a sudden change in the spring stiffness, and the coupled atomic-continuum model. The errors that are quantified between the two models are defined in terms of atomic displacement in Ω a . The exact displacement is the one obtained by the fully atomic model (d) and the approximation (d h ) obtained by the coupled model. Two quantities of interest are studied: firstly (Q 1 ) defined as atomic displacement, and the second one (Q 2 ) as average force on atom. The associated modeling errors are

ε i = Q i (d) -Q i (d h )
, where i = 1, 2. This convergence study shows decrease of the mentioned modeling errors with the increase of the Ω a . This analysis was a first step 6 , and a basis for the development of the adaptive strategy.

Finally, the Goals algorithm is extended to the Arlequin based coupled atomistic-tocontinuum modeling [START_REF] Dhia | A new adaptive modeling strategy based on optimal control for atomic-to-continuum coupling simulations[END_REF][START_REF] Prudhomme | An adaptive strategy for the control of modeling error in two-dimensional atomic-to-continuum coupling simulations[END_REF]. The adaptive procedure that controls the error is obtained by generating a sequence of surrogate problems so that the modeling error satisfies:

e = Q(d) -Q(d h ) ≤ γ tol , (26)
where γ tol is predefined tolerance. Reduction of the modeling error at each iteration is done by locally enriching the surrogate model, i.e. by locally switching on the atomic model in the subregions where the continuum model is not accurate enough 7 . Modeling error is defined globally over the whole domain but can be decomposed into local contributions (subdomains). Naturally, the elements of the finite element mesh used to discretize the continuum model are chosen as subdomains or cells. Now, some user-defined parameter (subdomain tolerance similar as γ tol ) is chosen to decide when can the subdomain be switched from the continuum model to the particle model.

MS methods comparison

In the foregoing, we have given an overview of the mainstream MS methods in terms of QC method and BD/Arlequin based coupling. The latter currently attracts the greatest attention with many BD/Arlequin recent developments, but it still has not been fully completed as the simpler, but well known QC method. Therefore, we try to present in this section the comparison of these two methods, hoping to be able to draw lessons on further improvements to the present practise. This comparison is carried out regarding: 1. coupling algorithm, 2. continuum modeling, 3. applicability, and finally 4. adaptivity strategy.

Coupling algorithm

The coupling algorithm of these two methods are drastically different. QC method seeks to provide a gradual transition, where the mesh composed of repatoms as nodes is gradually refined starting from the local towards the non-local description. This gradual transition is numerically more convenient regarding its capability to reduce the ill-conditioning. However, it also has a few drawbacks. First of all, an enormous refinement has to be performed in going from the FE continuum representation to the atomistic lattice size. Furthermore, the FE nodes and the atoms have to coincide. Contrary to that, BD/A method couples the two models only in the zone of partial overlap. Neither gradual transition nor coincidence between the nodes and elements are needed. However, atomistic and continuum DOFs are completely separated and additional unknowns in terms of Lagrange multipliers that enforce the coupling need to be accounted for. In addition, in order to avoid double counting the blending of the energy in overlapping domain is done by weighting functions, which also have to be chosen appropriately.

Continuum modeling

QC method works with large deformation intrinsically. That is, CB rule is used for continuum constitutive relation thus constitutive law is based on atomistics rather than on an assumed phenomenological form. On the other hand CB hypothesis is satisfied only for simple lattice structures8 . On the other hand, due to use of classical coupling of atomic and continuum domains, The BD/A method offers another approach to defining the surrogate (homogenized) continuum model. Namely, fitting the material parameters by virtual experiments on RVE, the atomistic part is reduced to the corresponding continuum. Thus, there is no need for CB hypothesis to be satisfied.

Applicability

During development period of the QC method, it served both as a key vehicle for understanding of the nature of atomistic-continuum coupling, and as a practical tool for investigating problems requiring coupled atomistic-continuum solution procedure. Nowadays, there is a unified web site qcmethod.org as the original source of information, with publications and the most important download section. Under the download section the QC code is available written in Fortran90 by the Tadmor and Miller. On the other hand, BD/A method was not used that much as a practical tool (apart some application to carbon nano tubes by Xiao and Belytschko), it was more used for the theoretical testing of different aspects of coupling and MS modeling in general. However the method is from the very beginning extended to dynamics, dealing with spurious wave reflections in the transition from the atomic to continuum domain. There is no unified web site as for QC method, but there are examples like libmultiscale.gforge.inria.fr.

Adaptivity strategy

Original QC method is in essence an adaptive FE approach, and adaptivity is intrinsically in the formulation in QC method. BD/A method was initially assumed as approach to couple two different models. However, the described evolution associated with the goal oriented error estimate theory, with the strong mathematical foundations, improved the method so that it is equal if not better compared to the QC method in the sense of model adaptivity. In particular, the choice where to place the fine and where to remain with coarse scale model, and how to provide the appropriate evolution of that region is still the most important question. More precisely, the BD/A method adaptivity is driven by the goals algorithm, controlling the model refinement with respect to the any chosen quantity of interest. In QC method non-locality criterion is based on a significant variation in the deformation gradient (no other criteria). In very recent contributions the model adaptivity is being combined with optimal control theory and shape optimization allowing size and shape of the zone of interest to be automatically determined (or controlled in the sense of the error in the quantity of interest).

Numerical examples with model adaptivity

The idea of model adaptivity described in the previous chapters is shown schematically in the following table for the simplest 1D case. Even though this procedure is similar to a mesh refinement (especially for the QC approach), the main goal is to address the model comparison that allows us to substitute the continuum model with the atomic one.

QC BD/A

As described in QC section adapting process advances by selecting new atoms as repatoms/nodes in the area where deformation gradient changes severely. This procedure is illustrated for the 1D scheme below.

The idea in the BD/A based methods is to adapt i.e. switch from continuum to atomistic model cell by cell in order to deliver accurate results regarding the selected QOI.

The overlap region has to be reconfigured. In the examples the follow, the quantities are selected for the sake of overview: Q1displacement of the rightmost node, Q2 -L2 norm of displacement error in overlapping zone, Q3 -mean strain in the overlapping zone, Q4 -L2 norm of strain error in overlapping zone, Q5 -stress difference between neighbouring bonds. Next to QOIs, some other parameters of the model ought to be selected selected, and properly adapted. These parameters are divided in two groups. The first one pertains to the shape i.e. size of the overlap along with the size of the FEs for continuum domain. The second group concerns the size of the fine scale (particle) domain i.e. to issue where to place the overlap.

Adapting model topology -FE and overlap size

In this example, we demonstrate the influence of the model topology on the accuracy of QOIs Q1 . . . Q4. The accuracy improvement is performed in the simple, chainlike (1D) problem as in Fig. 7 (see Section 4 and [START_REF] Marenic | Nano-submodelling technique based on overlapping domain decomposition method[END_REF]). Two parameters are taken into consideration for the ... topology adaptation: the size of the FE (l e ), and the size of the overlapping zone (l b ). In extension, local (only k 1 ) and nonlocal (k 1 and k 2 , see Fig. 7) types of interaction are selected in the atomistic domain and taken as the third parameter.

Parameter 1: the size of the FE l e

For the local, non-interpolated case of the model problem the solution of the full molecular and coupled model show no error in the mentioned QOIs with the variation of the parameter l e . On the other hand, for the case of interpolated coupling within the local interaction in Ω a , the error in the QOIs exists and it is shown in the Fig. 8 with respect the variation of the parameter l e . In the analysis of the influence of parameter l e upon the QOIs, the other parameter (the size of the overlapping zone) is kept constant (l b = 16l 0 = cst.). On the

e r,d = (u n -u ex n )/u ex n , (27)
where u n and u ex n is the displacement of the rightmost node and the exact value for the displacement, respectively. Note that in the examples presented herein, the exact values refer to the fully molecular (or particle) solution. On the Fig. 8 b) relative L 2 norm of the displacement error in the overlapping zone is given as

e L 2 ,d = i∈Ω b (d i -d ex i ) 2 i∈Ω b (d ex i ) 2 , (28)
where d i and d ex i are atom displacement solution (∀i ∈ Ω b ) for the coupled model and the exact solution, respectively. The relative error in Q3, the mean strain in the overlapping zone, is given on Fig. 8 c) as

e r,ε = (ε -εex ) εex , (29)
where ε and εex are the mean strain in overlapping zone and exact mean strain, respectively. Likewise, on the Fig. 8 d) relative L 2 norm of the strain error in the overlapping zone is given as

e L 2 ,ε = i∈Ω b (ε i -ε ex i ) 2 i∈Ω b (ε ex i ) 2 , ( 30
)
where ε i and ε ex i are strain solution (∀i ∈ Ω b ) for the coupled model and the exact solution, respectively.

Similarly, the same analysis is performed with the non-local interaction in atomistic domain. The results are shown in the Fig. 9 indicating the same tendency as for local interaction with the bigger error. Note that for all the plots in Figs. [START_REF] Belytschko | Coupling methods for continuum model with molecular model[END_REF] and Fig. 9 the errors in QOIs drops down to zero as the size of the FE decreases to lattice constant being l e = l 0 . This result is rather logical, since decreasing the FE size for the interpolated coupling case we approach the non-interpolated case (see Fig. 6) where no error occurs, as already mentioned above.

Parameter 2: the size of the bridging zone l b

The FE size is varied here together with the size of the bridging zone, as in Fig. 10 a), keeping FE size equal to overlap size (l e = l b as in Fig. 10 a)). This is because if the size of the FE is kept constant with the variation of the l b (see Fig. 10 b)) then the influence of the number of the FE in the bridging zone is notable (as studied in the section above). For the local, non-interpolated case of the model problem, the solution of the full molecular and coupled model show no error in the mentioned QOIs with the variation of the parameter l b . The same goes for FE size as parameter. For the case of interpolated coupling, the error in the QOIs is shown in the Fig. 11 with respect to the variation of the parameter l b , only for the simpler case of local interaction. Note that the diagrams a) and c) on Fig. 11 show that the error in displacement of the end node and the error in mean strain drops with the increase of the l b .

Adapting the position of the overlap

In this section the parameter for adapting the coupled model is the position of the overlap zone with respect to the strain gradient caused by distributed load or by hypothetical defect. Actually, the position of the overlap is directly related to the fine scale model size, that is to the question pertaining to the size of the fine-scale model.

Model with distributed load

A model with the distributed load spreading in the particle domain is chosen to analyse the influence of the overlap position on the accuracy of the certain QOIs. Three different configurations are considered as shown in Fig. 12. The two limiting cases, one where the distributed load is completely in the atomistic domain (spreading also in the overlap called case 3)) and the other with distributed load only in atomistic domain but not in overlap (denoted as case 1)). Case 2) concerns configuration(s) between. The error in QOIs Q1 and Q2 versus mentioned three cases is plotted on Fig. 13. Not quite surprising, the presented results show better accuracy in terms of selected QOIs as the particle size is increased (i.e. as the distributed load is further from the overlap). Furthermore, QOI Q5 representing stress difference between neighbouring bonds is taken as the control variable to adapt the fine scale size. The relative error in Q5 is defined as where ∆σ ex and ∆σ ex is the exact stress difference and the one obtained from coupled model, respectively (see Fig. 14 a)). Stress difference is obtained as ∆σ = |σ i -σ i+1 |, i.e. the difference of stress (piece-wise constant) in the neighbouring bonds. Results of the relative error in stress difference of the leftmost atom in the overlap versus the position of overlap is presented in Fig. 14 b). The results show that the error in terms of stress QOI decreases with the increase of the size of fine scale model. Clearly, when the strain gradient, caused by the distributed load, is in fine scale model completely the error in stress QOI does not exist. This QOI provides a very good local refinement criterion. We note in passing that such a QOI, apart from being a good refinement criteria, can be related to the often mentioned ghost forces problem. Thus, choosing this QOI Goals algorithm can be used to iteratively adapt coupled model to increase the coupling quality i.e. ghost forces. This is more important for complex interatomic potential with non-local interaction, which is the subject of forthcoming research. rather, to illustrate several useful choices for the parameters to adapt. Likewise, the choice of functional defining quantities of interest is not fully exhaustive. Different choices of the parameters and the quantities of interest made herein are used in order to illustrate that they remain problem dependent. Further perspectives of the BD/A method development is related to the implementation of the complex multi-body potential which enables more realistic description of the discrete model. In many such case, the complex potential is even more computationally demanding as elaborated in the Section 2, which additionally justifies using the MS strategy. The use of this kind of potentials enables modeling of inelastic behavior and localized failure at the nano-scale. Equally, those potentials are able to describe the behavior of 'living' material in life science. Both of this issues are presently very important in the technical and biomechanics research domain (e.g. see Kojic et al. [START_REF] Kojic | On the application of discrete particle methods and their coupling to the continuum based methods within a multiscale scheme[END_REF][START_REF] Kojic | Comuter Modelling in Bioengineering Theory, Examples and Software[END_REF]). 3)

e r,σ = (∆σ ex -∆σ)/∆σ ex , (31)
Figure 12: Three cases of the position of the bridging zone with respect to the distributed load 1) distributed load (q) not in overlap, 2) q partially in overlap and 3) q on all atoms, completely covering the overlap.

Figure 1 :

 1 Figure1: Lennard-Jones and Harmonic potential (dashed line). Note that the Harmonic potential is a suitable approximation when the particles are around the equilibrium position.

  3 and indices a and b (a = b) refers to the neighboring elements of rep-atom. The rep-atom will be made local if this inequality is satisfied, and non-local otherwise, depending on the empirical constant .

Figure 2 :

 2 Figure 2: Scheme of the coupled model and the nSM technique.

Figure 3 :

 3 Figure 3: Scheme of the lattice sample and increasing square shaped RVEs (1, 2, . . . , i, . . . , n) increasing until material parameters convergence. Due to simplicity, only RVE 1 is shown as a lattice structure.

Figure 4 :

 4 Figure 4: Atomistic strain energy density (W a ) convergence with increasing size of RVE.

Figure 5 :

 5 Figure 5: Energy weighting function distribution in the bridging zone.

Figure 6 :

 6 Figure 6: Scheme of the distribution of the LM nodes for strict and interpolated coupling.

Figure 7 :

 7 Figure 7: 1D coupling BD/A based model scheme with the symmetry BC on the left end of the atomistic domain.

4 FE 3 FEFE 25 FE 2 ,εFigure 8 :

 432528 Figure 8: Local interaction in Ω a with FE size l e as a parameter. Quantities of interest Q1, Q2, Q3 and Q4 are shown on plots denoted as a), b), c) and d), respectively.

4 FE 3 FEFE 35 FEFigure 9 :

 43359 Figure 9: Nonlocal interaction in Ω a with FE size l e as a parameter. Quantities of interest Q1, Q2, Q3 and Q4 are shown on plots a), b), c) and d), respectively.

Figure 10 :

 10 Figure 10: Options for the study of the influence of the bridging zone size parameter. a) l e = l b and b) l e = cst.

2 ,εFigure 11 :

 211 Figure 11: Local interaction in Ω a with size of Ω b (l b ) as a parameter. Quantities of interest Q1, Q2, Q3 and Q4 are shown on subplots a), b), c) and d), respectively.

Figure 13 :

 13 Figure 13: Local interaction in Ω a with position of distributed load as a parameter (for L2 and H1 coupling, see eq. 24). Quantities of interest Q1 and Q2 are shown on plots a) and b), respectively.

Figure 15 :

 15 Figure 15: Modeling of defect by the sudden spring stiffness drop located on the left end a), and characteristic cases regarding the overlap position (d 0 ) with respect to the defect radius (R def ) used to illustrate adaptive process b).

Figure 16 :

 16 Figure 16: Local interaction in Ω a with position of the defect (d 0 ) as parameter. QOI Q2 is shown for the four variants of coupling (strict, interpolated, L2 and H1)

  We present further some numerical examples that can clearly demonstrate the model adaptivity for the BD/A based coupled model. The accuracy of chosen quantities of interest (QOI) is used as the measure of the model adaptivity performance. Contrary to QC method, there are many candidates for local QOIs, and the best choice certainly depends on the problem on hands.

	-rep-atom	-atom driven	-selected to	-fully atomistic	-atom in	-FE node
		by continuum	be rep-atom	model	overlap
		. . .			. . .

In[START_REF] Miller | Quasicontinuum models of fracture and plasticity[END_REF][START_REF] Miller | Quasicontinuum simulation of fracture at the atomic scale[END_REF] the QC method has been applied to crack tip deformation and was shown to account for both the brittle fracture and ductile crack tip dislocation emission.

Note that simply having a large deformation in a region does not in itself require a non-local rep-atom, as the CB rule of the local formulation will exactly describe the energy of any uniform deformation, regardless of its size.

In the QC formulation the constant strain triangle (CST) elements are used with the linear shape functions to interpolate the displacement field within each element. In this case the deformation gradient will be uniform.

Taking continuum model behavior to be governed by Hooke's law, it is possible the calibrate Young's modulus (only the converged value is taken into further consideration).

This fact is rather obvious as increasing the atomic domain corresponds to setting the overlap region away from the defect, that is, to couple the two models in a region in which they are supposed to be compatible.

No switching back to the coarse model is described in the mentioned references, whereas in QC method the procedure of removing unnecessary repatoms from the mesh is studied.

CB hypothesis is satisfied only for simple lattice structures. However there is an extension of CB rule[START_REF] Sunyk | On higher gradients in continuum-atomistic modelling[END_REF][START_REF] Sunyk | On Aspects of Mixed Continuum-Atomistic Material Modelling[END_REF].

Model with defect

Next, a model with the defect is analyzed. This defect is modeled as the sudden stiffness change (see Fig. 15 a)) which occurs inside the particle domain. Problem is similar as the distributed load but with a more severe strain gradient. According to the adaptive scheme in Fig. 15 b) the fine-scale model size is increased. Not surprisingly, adapting the model in a way that the defect (severe strain gradient) is included in fine scale model, reduces the error in QOIs as can be seen in Fig. 16.

Conclusions and perspectives

In the consideration where nano-scale effects are important, the reference solution can be obtained by using the full atomic model relying upon interatomic potentials to provide the results of interest. However, due to the complexity of engineering problems and the corresponding scales that appear in realistic problems, explicit modeling of all with only the atomic degrees of freedom will very likely never be feasible. Thus, one must reduce the size of the problem by multiscale strategies (MS) that selectively removes most of the degrees of freedom by homogenized continuum formulation in order to make the problem solution tractable.

In this survey, we have discussed the salient features, similarities and some recent developments in the most frequently used strategies that selective reduce the number of degrees of freedom, BD/A and QC methods. In the references related to the BD/A MS methods authors usually concentrate on the atomistic-to-continuum coupling performance. Thus, selective removing of degrees of freedom i.e. model adaptivity is neglected. The QC method is based on adaptive approach, thus being an exception and the reference for comparison in this review survey. The evolution of the BD/A method from atomistic-to-continum coupling to adaptive MS method is presented, as well as the literature regarding error estimation theory which parallels the development.

The general idea of model adaptivity is demonstrated on few numerical examples. The presented examples deal with the simplest 1D case, and they should not be used to quantify computational efficiency or the limits of adaptive criteria (tolerances). The idea was,