{'fraction_non_alphanumeric': 0.03068825106763333, 'fraction_numerical': 0.007175489125037243, 'mean_word_length': 4.555448275862069, 'pattern_counts': {'":': 0, '<': 4, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
An Analysis of the Quasicontinuum Method

The aim of this paper is to present a streamlined and fully three-dimensional version of the quasicontinuum (QC) theory of Tadmor et al. and to analyze its accuracy and convergence characteristics. Specifically, we assess the effect of the summation rules on accuracy; we determine the rate of convergence of the method in the presence of strong singularities, such as point loads; and we assess the effect of the refinement tolerance, which controls the rate at which new nodes are inserted in the model, on the development of dislocation microstructures.


Introduction
The aim of this paper is to present a streamlined and fully three-dimensional version of the quasicontinuum (QC) theory of Tadmor et al. [18,19] and to analyze its accuracy and convergence characteristics. The theory of the quasicontinuum furnishes a computational scheme for seamlessly bridging the atomistic and continuum realms. The chief objective of the theory is to systematically coarsen an atomistic description by-and only by-the judicious introduction of kinematic constraints. These kinematic constraints are selected and designed so as to preserve full atomistic resolution where required, e. g., in the vicinity of lattice defects, and to treat collectively large numbers of atoms in regions where the deformation field varies slowly on the scale of the lattice. Thus, in its purest form all input into the theory concerning the behavior of the material is atomistic, and all approximations are strictly kinematic in nature.
Different variants of the theory have been developed and documented over a series of publications [18,19,12,10,9,11,13,17,14,16], where numerous examples of application have also been presented. Here we endeavor to emphasize the essential building blocks of the static theory, to wit: i) the constrained minimization of the atomistic energy of the solid; ii) the use of summation rules in order to compute the effective equilibrium equations; iii) and the use of adaption criteria in order to tailor the computational mesh to the structure of the deformation field. In particular, we develop a new class of summation rules based on sampling lattice functions over clusters of atoms. In Section 4, we present a detailed numerical analysis that probes various aspects of the accuracy and performance of the method. Specifically, we assess the effect of the cluster size and lumping procedure on accuracy; we determine the rate of convergence of the method in the presence of strong singularities, such as point loads; and we assess the effect of the refinement tolerance, which controls the rate at which new nodes are inserted in the model, on the development of dislocation microstructures.

The quasicontinuum method
We consider a reference configuration of a crystal such that its N atoms occupy a subset of a simple d-dimensional Bravais lattice spanned by basis vectors {a i , i = 1, . . . , d}. The coordinates of the atoms in this reference configuration are: Here, l are the lattice coordinates which designate individual atoms, Z is the set of integer numbers and d is the dimension of space. From the definition, it follows that the set L is the collection of lattice sites occupied by atoms. The coordinates of the atoms in a deformed configuration of the crystal are denoted {q(l), l ∈ L}. For convenience, we shall collect all atomic coordinates in an array q and regard such array as an element of the linear space X ≡ R N d , the 'configuration' space of the crystal. The energy of the crystal is assumed to be expressible as a function E(q), e. g., through the use of empirical interatomic potentials. In addition, the crystals may be subjected to applied loads. We assume these loads to be conservative and to derive from an external potential Φ ext (q). The total potential energy is, therefore: The crystal may also be subject to displacement boundary conditions over part of its boundary, e. g., as a result of the application of a rigid indentor. We wish to determine the stable equilibrium configurations of the crystal. The problem is, therefore, to determine all the local minima of Φ(q) consistent with the displacement, or essential, boundary conditions. We shall, somewhat equivocally, enunciate this problem as: We emphasize, however, that in general the aim is not simply to determine the absolute minimizer of Φ(q), but rather the much richer, and physically more relevant, set of metastable configurations of the crystal. The distinction arises from the fact that the energy function E(q) must be invariant under the class of affine deformations which map the crystal lattice onto itself. This implies, for instance, that E(q) must be periodic under crystallographic slip with period equal to the Burgers vector. Hence, the energy function E(q) is strongly nonconvex and the potential energy Φ(q) in general has vast numbers of local minimizers, or metastable configurations. Techniques for systematically exploring the full set of metastable configurations, and indeed the entire configuration space of the crystal, include statistical mechanics and path-integral methods. An alternative approach is to load the crystal incrementally and proceed by continuation. Thus, at every step in the solution procedure the loads, or the prescribed displacements, are incremented and the system is allowed to relax to a nearby stable configuration. The relaxation process may be governed by inertia, viscosity, kinetics, or some other rate-limiting mechanism; or may be numerical in nature and owe to the use of iterative solvers such as conjugate gradients. This latter approach is adopted in all the numerical tests presented in this paper.

Interpolation
The essence of the theory of Tadmor et al. [18,19] is to replace (3) by a constrained minimization of Φ(q) over a suitably chosen subspace X h of X. In order to define X h , we begin by selecting a reduced set L h ⊂ L of N h < N 'representative atoms' (Fig. 1). The selection of the representative atoms is based on the local variation of the fields and will be addressed subsequently. In addition, we introduce a triangulation T h of L h . It bears emphasis that the triangulation T h may be unstructured. In particular, L h need not define a Bravais lattice. The positions of the remaining atoms are determined by piecewise linear interpolation of the representative atom coordinates. We shall regard the resulting coordinates {q h (l), l ∈ L} as belonging to a linear space X h of dimension N h d.
Let ϕ h (X|l h ), l h ∈ L h , be a collection of shape functions for T h . Thus, ϕ h (X|l h ) is continuous and piecewise linear, its domain is restricted to the simplices K ∈ T h incident to X(l h ), and it vanishes at all nodes of the triangulation except at X(l h ), where it takes the value 1, i. e., Figure 1: Example of triangulation T h of the crystal.
By construction, where we write Evidently, {ϕ h (l|l h ), l h ∈ L h } constitutes a basis for X h and the fields q h (l) are entirely determined by their values q h (l h ) at the representative atoms. In addition, the basis lattice functions are required to satisfy the identity: i. e., the basis lattice functions must define a partition of unity over L. This requirement ensures that constant fields are interpolated exactly by the basis lattice functions.

Reduced equations
The reduced counterpart of problem (3) is now The minimizers of the reduced problem satisfy the reduced equations of equilibrium Here, are the forces corresponding to q and f (l|q) is the value of f (q) at site l. Thus, the reduced problem entails the solution of the N h d equations (9) in the N h d unknowns q(l h ), l h ∈ L h .

Node-based summation rules
The practicality of the method further hinges on the possibility of avoiding the calculation of the full atomistic force array f and carrying out full lattice sums, as seemingly required in (9). As noted by Tadmor et al. [18,19], this may be accomplished by the introduction of summation rules similar to the conventional quadrature rules of numerical integration. The problem is thus to approximate sums of the general form where g(l) is a lattice function. By analogy to numerical quadrature rules, we begin by exploring summation rules of the form for some suitably chosen collection of summation points S h and weights n h (l), not necessarily integer. Loosely speaking, the weights n h (l) may be regarded as the number of atoms represented by the site l. Proceeding as in the development of numerical quadrature rules, the weights n h (l) may be determined by the requirement that the summation rule (12) be exact for a restricted class of lattice functions [13]. The lowest-order summation rule is obtained by requiring that all lattice functions in X h , i. e., all piecewise linear functions supported by the triangulation T h , be summed exactly. This is tantamount to setting S h = L h and requiring that the summation rule (12) be exact for all shape functions ϕ h (l|l h ), l h ∈ L h . This requirement gives, explicitly, In fine regions of the triangulation, the sums on the right-hand side of (13) may be computed explicitly. By contrast, in coarse regions of the mesh approaching the continuum limit this explicit calculation becomes impractical. However, in such regions the lattice sum (12) ostensibly reduces to an integral and the corresponding weights are those of conventional Lobatto quadrature [2]. In this limit, each simplex K ∈ T h simply contributes (N/V )|K|/(d + 1) atoms to each of its d + 1 nodes, where N/V is the atom density in the undeformed configuration of the crystal. As a simple illustrative example, consider a monoatomic chain discretized into 'elements' each containing L bonds. For an internal node, the summation weight follows as which coincides with the continuum limit. By contrast, for an end node the summation weight is We note that n h tends to the continuum Lobatto limit of L/2 as L → ∞ and reduces to n h = 1 in the full atomistic limit, as required.

Reduced equations based on summation rules
The application of a summation rule to the evaluation of the equilibrium equations (9) leads to the system of equations As desired, the calculation of the effective forces in (16) is of complexity O(N h ) provided that the number of sampling sites in S h is of order N h and the atomic interactions are short-ranged. The formulation of the method takes a somewhat unexpected turn at this point in that the node-based summation rule discussed in the foregoing suffers from a rank-deficiency problem. Thus, in the linear case the node-based summation rule leads to rank-deficient systems of equations. This is in analogy to the rank-deficiency of finite-element stiffness matrices computed by Lobatto quadrature, and manifests itself as the existence of socalled zero-energy modes which pollute the displacement field or render the system of equations singular altogether. By way of illustration, we consider a Lennard-Jones crystal in the form of a cube containing 32 × 32 × 32 fcc unit cells. In this example, eight representative atoms are introduced at the corners of the cube, and the cube is triangulated as shown in Fig. 2a. The stability of the crystal in its reference configuration may be assessed by examining the spectral properties of the reduced stiffness matrix. A direct calculation of its eigenvalues reveals that seven of them are identically equal to zero, which in turn denotes the presence of one zero-energy deformation mode, Fig. 2b. However, the zero-energy deformation mode is absent when the lattice sums are carried out in full, i. e., when the test is repeated with the equilibrium equations (16) replaced by (9).
The numerical test just described shows that the node-based summation rules are indeed rank-deficient, but it also suggests that the reduced system is stable provided that a sufficient number of sampling points is included in the summation rules. A new set of summation rules which offers this added flexibility is presented next.

Cluster summation rules
A class of summation rules which generalizes the node-based rule may be obtained by sampling the lattice function over neighborhoods of the representative atoms. We shall refer to these neighborhoods as clusters, and the resulting summation rules as cluster summation rules. Each cluster may be regarded as a representative crystallite where the state and the behavior of the crystal are sampled.
More specifically, let C(l h ) = {l : |X(l) − X(l h )| ≤ r(l h )} be the cluster of lattice sites located within a sphere of radius r(l h ) centered on the representative atom l h (Fig. 3). Assume for now that the clusters thus defined do not overlap. The cluster summation rule induced by the set of clusters where S(l h ) denotes the sum over all atoms in the cluster C(l h ), i. e., As before, the cluster weights n h (l h ), l h ∈ L h are computed by requiring that the cluster summation rule (17) be exact for all basis functions. It should be carefully noted that, as the triangulation size approaches the atomic length scale, the spherical clusters defined previously in general overlap, Fig. 4. In this case, we truncate the clusters as follows. A site l belongs to cluster C(l h ) if its distance to X(l h ) is less than a prescribed value r(l h ) and less than its distance to any other representative atom. Thus, the cluster C(l h ) is the intersection of the Voronoi cell containing X(l h ) and sphere of radius r(l h ) centered at X(l h ). Ambiguous cases corresponding to sites which are equidistant from two representative atoms are resolved randomly. In the limit of full atomistic resolution, each cluster contains exactly one atom and the corresponding weight is 1. Thus, the cluster summation rules reduce seamlessly to the exact lattice sum when the fully atomistic limit is attained.

Summation-rule error analysis for monatomic chains
The accuracy of cluster summation rules can be analyzed simply in the case of a onedimensional monatomic chain. In this case, L ≡ Z is the set of all sites in the chain. Without loss of generality, we take the coordinate of a site l ∈ L in the undeformed configuration of the chain to be X(l) = l. For simplicity, we consider a uniform triangulation of L of size h. Thus, the coordinates of the representative atoms are where h is the simplex size. We consider clusters of uniform radius r < h/2, Fig. 5. The summation weights follow from the requirement that the basis functions ϕ h (l|l h ) be summed exactly, with the result, Thus, for the infinite monatomic chain, n h (l h ) is simply the ratio of the number of atoms in a simplex to the number of atoms in a cluster.
In order to obtain the summation error to leading order, we consider the sum of a quadratic polynomial supported on two adjacent simplices in the reduced chain, namely, The exact sum of g(l) over all sites of the chain is The calculation of the approximate sum S h requires three local cluster sums The use of the cluster summation rules leads to the approximate lattice sum The summation error is |S h −S|. An upper bound on the error is obtained by setting r = 0. Then, for h ≫ 1 we have |S h − S| r=0 ∼ (4/3)h 3 , which shows that the summation rule is of third order. The effect of the cluster size on the accuracy of the summation rule may be understood as follows. For simplicity, we consider the case of h ≫ 1, whereupon (21) simplifies to and (25) to where we write ξ = r/h. In the regime of ξ ≪ 1, (27) simplifies further and the corresponding relative error behaves as which is a decreasing function of r. Thus, as expected, increasing r from 0 increases the accuracy of the summation rule. In the regime of r ≫ 1 and r ∼ h, the relative error behaves as which vanishes for ξ = 1/2, or r = h/2. Thus, the exact sum is obtained when the clusters encompass the entire lattice, as required.

Reduced equations based on cluster summation rules
The application of the cluster based summation rule to the reduced equilibrium equations (9) yields The computational complexity of (30) is O(N h N r ), where N r is the number of lattice sites in a cluster of radius r. For N r > 1, the computational effort is in excess of that corresponding to the node-based summation rules. Thus, the optimal value of the cluster size r(l h ) is subject two opposing requirements. On one hand, the cluster size has to contain a sufficient number of atoms to ensure the appropriate accuracy and stability of the summation rule. On the other, however, such number should remain small for the method to retain its computational efficiency. The effect of the cluster size on the accuracy is investigated subsequently by way of numerical testing.

Adaptive selection of representative atoms
The third key component of the quasicontinuum method is the use of mesh adaption in order to tailor the computational mesh to the structure of the deformation field. Ideally, the adaptivity itself should be driven by the energetics of the system, i. e., the mesh should return the least possible potential energy for a fixed number of representative atoms. However, in order to optimize the mesh in this manner it is necessary to know the relation between the energy, or suitable energy bounds thereof, and the mesh size. Unfortunately, unlike the continuum case, an approximation theory for discrete systems such as considered here appears to be entirely lacking at present, and such energy bounds are presently unavailable.
Under these conditions, it becomes necessary to resort to empirical adaption indicators. In calculations we adopt as adaption indicator ǫ(K) for simplex K the variation of the displacement field over K modulo rotations. This variation should be well-defined since the displacement field of the crystal may reasonably be expected to be of bounded variation. Specifically, we measure the variation of the displacement field as where II E (K) denotes the second invariant of the Lagrangian strain tensor in simplex K, and h(K) is the size of K (cf [18,19]). It follows from its definition that ǫ(K) is invariant under rotations. The element K is deemed acceptable if for some prescribed tolerance T OL < 1, and is targeted for refinement otherwise. In (32) b denotes the magnitude of the smallest Burgers vector of the crystal. In our implementation of the method, the refinement of the simplices which violate (32) is accomplished by the application of subdivision rules. In particular, we choose to bisect the longest edge of the simplices which are targeted for refinement. The new representative atoms are inserted in the lattice site nearest to the geometrical midpoint of the longest edge. In order to preserve the quality of the computational mesh, we apply standard mesh-improvement operations such as edge-face or octahedral swapping [3,4,1] and smoothing [15]. Each new mesh is equilibrated and the remeshing criterion (32) is reevaluated for the new solution. The process is repeated until the mesh remains unchanged.
The significance of (32) becomes apparent by considering a processes of crystallographic slip across K. To this end, imagine that a pair of representative atoms in K undergo a relative sliding displacement of magnitude b across a slip plane of the crystal. For this deformation one has It follows, therefore, that the adaption criterion (32) ensures that the mesh size h(K) < b under the conditions just described. Thus, the adaption criterion is designed so that full atomistic resolution is attained when a simplex slips by a full Burgers vector. Interestingly, we have found that the dislocation patterns predicted in calculations are sensitive to the choice of tolerance, and generally T OL must be chosen much smaller than 1 in order not to inhibit dislocation nucleation.

Details of the computer implementation
In order for the quasicontinuum method to perform optimally, some implementational issues require careful attention. Several potential performance bottlenecks are discussed in this section.

Site-element mapping
The numerical solution of the quasicontinuum equilibrium equations (16), requires repeated computation of atomic coordinates q h (l), l ∈ L which depend, through the interpolation constraints (5), on the coordinates of representative atoms q h (l h ), l h ∈ L h . It is therefore important to devise an efficient algorithm for locating the simplex K ∈ T h which contains a given site l of the crystal. The simplex K ∈ T h containing site l ∈ L may be found from elementary geometry and an exhaustive search over T h . Such method carries, however, substantial computational cost, due to a large number of required floating-point operations. Our approach is based on the fact that once the pair {l, K} has been established, it remains valid for the lifespan of T h . Thus, the following two stage strategy may be adopted: 1. A general search routine, which for a given lattice site l ∈ L and triangulation T h returns K ∈ T h containing l. 2. One must establish a strategy for dealing with cases of two distinct keys k 1 and k 2 for which the resulting values of the hash function are the same, i. e. h(k 1 ) = h(k 2 ) (collision resolution).

A look-up
In his monograph, Knuth [6] provides an excellent description of hashing, including some possible choices of hash functions. For a simple three-dimensional Bravais lattice (1), the lattice coordinates of a site l = {l 1 , l 2 , l 3 } ∈ Z 3 may be considered as a key. The hash function acting on this key produces an index into a table of simplices where the K, containing l, can be found. In cases in which a key contains more than one word (each of l i can be considered as a separate word) Knuth [6] suggests using the following hash function where h i (k) denotes a hash function for l i and M is an integer parameter. The function a mod p computes the remainder of a/p. In general, the optimal choice of h i depends on T h and cannot easily be established. We have found, however, that the particular hash function works well for a wide range of triangulations. The meaning of lshift(l, m) is simply "shift all bits in l left m positions". The operations involved in the calculation of h i are elementary, and the computational cost of each such calculation is small. The integer parameter M in (34) may be chosen arbitrarily, but Knuth [6] points out that, when it is taken to be some power of 2, i. e., M = 2 p , a mod M is equivalent to masking the low p-bits from a (on most of currently available computer hardware). The application of (34) to a lattice site l results in an index into the hash table from the interval [0, 2 p ). Accordingly, the size of the table is limited to 2 p , which renders collisions likely. The collisions are handled with the help of additional short tables which store {l, K} pairs for all l such that the h(l) produces the same value. Therefore, the search through the cache for a simplex K containing site l becomes essentially a two stage process, in which the application of the hash function is followed by a linear search. Since, the number of lattice sites having the same index is small, and the additional cost due to the linear search phase is not significant.

Computation of cluster weights
As explained earlier, the cluster weights n h (l h ), l ∈ L h are obtained requiring that the shape functions are summed exactly by the cluster summation rule (17). The application of the summation rule to shape functions ϕ h (l|l h ) leads to the system of N h linear algebraic equations to be solved for n h (l h ), l ∈ L h . Here the matrix and C(l ′ h ) denotes the cluster of lattice sites centered at l ′ h . The array b(l h ), l h is simply The calculation of the array b requires the evaluation of the shape functions at all lattice sites within a crystal. While this may be regarded as acceptable for small crystals, it becomes prohibitively expensive when applied to large samples. However, as we have already shown, in the continuum limit summation may be replaced by integration, in which case simplex K contributes N (K)/(d + 1) to each of its vertices. N (K) = (N/V )|K| is the approximate number of sites within K and N/V is the atom density in the undeformed configuration of the crystal. Specifically, one may introduce a cutoff N c and restrict the direct calculation of the sum in (38) to simplices K for which N (K) ≤ N c . It is apparent from (37) that A has a sparse structure, which suggests the use of specialized solvers for sparse linear systems. However, a simpler alternative route is to resort to lumping in order to replace A by a diagonal matrix (see, e. g., [2]). A lumping technique which is widely used in finite elements is the row-sum technique, which gives the diagonal entries of the lumped matrix as with all other entries in A set to zero. Once the matrix A is lumped, the solution of (37) is trivial and gives In the atomistic limit, the values of cluster weights computed from (40) become identical with obtained exactly.

Numerical tests
The accuracy of the quasicontinuum method is largely determined by three factors: • The value of the cluster-size parameter r, which controls the accuracy of the cluster summation rules.
• The approximations introduced in the computation of the cluster weights, namely: the approximate calculation of vector b, controlled by the cuttoff parameter N c ; and the lumping procedure for constructing the diagonal matrix A. • The value of parameter T OL in (32), which controls the process of representative atom insertion.
In this section we present an analysis of the influence of each of these factors on the accuracy of the quasicontinuum method. We measure the quality of the solution in energy terms. Specifically, we identify the error in the approximate solution q h with where q is the solution of the full system obtained, e. g., by relaxing q h . Since the potential energy decreases during this relaxation, the energy error is greater or equal to zero. For a well-supported harmonic crystal the energy error defined in (41) defines a proper norm of the error lattice function q h − q.

Test problem definition
We take nanoindentation as a convenient test problem for assessing the performance of the method. A salient feature of nanoindentation is the presence of a highly nonuniform state of deformation resulting in a sharp mesh gradation away from the indentor. This feature effectively tests the adaptivity of the method. The test sample used in simulations is a cube of an fcc nearest-neighbor Lennard-Jones crystal containing 64 × 64 × 64 fcc unit cells, or a total of 1 073 345 atoms. The size of the cube is limited by the need to compute the solution of the full problem by direct atomistic simulation in order to evaluate the error. The Lennard-Jones potential is [7,8] where σ and ε are parameters which set the length and energy scales, respectively. The choice of the Lennard-Jones potential is motivated by the desire to eliminate surface effects which would otherwise compound the interpretation of results.
The surfaces of the sample are aligned with the cube directions, Fig. 6. Atoms in the bottom surface are constrained to remain at their initial positions throughout the test, whereas in all side surfaces the atoms are allowed to move in the z-direction only. No displacement constraints are introduced on the top surface of the cube.
In calculations we use a model of a spherical indentor proposed by Kelchner et al. [5]. In this model, the indentor is regarded as an additional external potential Φ ext interacting with atoms in the substrate. The potential has the particular functional form where R is the radius of the indentor, r denotes distance between a site and the center of the indentor, A is a force constant and H(r) is the step function. In calculations we adopt the following values of the parameters: The initial triangulation of the cube is specifically tailored to the nanoindentation geometry, Fig. 7. Thus, in the region of the crystal located directly underneath the indentor, full atomistic resolution is introduced from the outset. Away from this region, the triangulation becomes gradually coarser. The resulting number of representative atoms in the initial mesh is 888, a significant reduction from the total number of atoms (1 073 345) in the sample. A cross section of the initial mesh through the center of the cube with a plane y = const is show in Fig. 7. All solutions are computed using a nonlinear version of the Conjugate Gradient method.

Effect of the cluster size
The cluster radius r controls the accuracy of the cluster summation rules. In order to isolate the effect of the cluster size from other factors, in the present tests the summation weights are computed exactly and the mesh is kept fixed. The indentor is pushed through a distance 0.1 [σ] into the crystal along the z-axis. This value ensures that the material immediately under the indentor becomes highly deformed and is pushed well into the nonlinear regime. However, the induced deformation is not sufficient for any defects to arise. This greatly simplifies the interpretation of results and the error analysis, since under the stated conditions the full atomistic solution q to compare with is well-defined.
The computed dependence of the error e on the cluster size r is shown in Fig. 8. It is evident from this plot that the inclusion of a single shell of neighbors in the clusters suffices to eliminate the rank-deficiency of the node-based summation rule. As suggested by our previous analysis, the error is greatest for the smallest cluster size and decreases rapidly as the cluster size increases. For r above √ 8 [σ], corresponding to the radius of the eighth shell of neighbors, the solution becomes insensitive to the cluster size. Thus, the use of relative small clusters results in high accuracy while preserving the computational efficiency of the method.

Effect of the lumping procedure
The error in the calculation of the cluster weights originates from two sources: the approximate computation of the array b in (38), controlled by the value of parameter N c ; and the lumping of matrix A in (37). In order to appraise this error, we repeat the calculations described in the preceding section with the summation weights computed approximately. The cutoff N c is set to 2000, i. e. the contributions to the array b from simplices which contain fewer than 2000 atoms are computed explicitly. The dependence of the error e on the cluster size parameter r for both the approximate (squares) and exact (circles) weights is shown in Fig. 9. Small differences do arise for very small clusters, and virtu- ally disappear for large clusters. In view of these results we may conclude that the use of approximate summation weights does not result in a significant loss of accuracy.

Convergence properties of the quasicontinuum method
We proceed to assess numerically the rate at which the quasicontinuum solution converges towards the fully atomistic solution. The test case is as in the preceding sections. In all cases, the cluster size r is set to 2 √ 2 [σ] and the cluster summation weights are computed approximately with a cutoff parameter N c = 2000. The coarsest mesh is shown in Fig. 10a. Subsequent finer meshes are constructed by regular refinement, Fig. 10b, resulting in increasing numbers of representative atoms N h , with the exception of the fully resolved case, h = √ 2 [σ], in which the representative atoms are placed on all the crystal lattice sites. In these calculations, the crystal is loaded simply by imparting a downward displacement −0.1 [σ] to the central atom on the surface labeled A in Fig. 6.  These results demonstrate the convergence characteristics of the quasicontinuum method in the extreme case of an applied point load. When full lattice sums are performed the calculated convergence rate of 0.55 approaches the optimal rate of convergence of linear finite elements applied to linear elasticity, which is 2/3. The introduction of cluster summation rules results in a degradation of the convergence rate. It should be carefully noted, however, that some of these conclusions may well depend on the loading geometry. As is well known, the linear elasticity solution corresponding to a point load diverges under the point of application of the load and does not possess finite energy. The finiteness of the energy of atomistic solution owes entirely to the discreteness of the atomic lattice. In this sense, the test case considered here is, therefore, a worse case, and it is possible that the converge of the quasicontinuum method is more robust for smooth loading.

Effect of the remeshing-indicator tolerance
The lack of convexity of the energy functional for crystalline materials allows for lattice defects to develop. Among these defects, dislocations play a crucial role, as the carriers of the plastic deformation. One of the appealing features of the quasicontinuum method is its ability to follow the nucleation and motion of dislocations inside the crystal. This is achieved by adaptively providing full atomistic resolution in the highly deformed regions of the crystal. In the nanoindentation test, dislocations nucleate under the indentor upon the attainment of a critical load and subsequently propagate into the crystal.
The insertion of new representative atoms, leading to mesh refinement, is controlled by criterion (32) and, specifically, by the tolerance T OL. A small value of T OL promotes refinement, whereas a large value of T OL inhibits refinement. Most importantly, an excessively large value of T OL may have the undesirable effect of inhibiting the nucleation of dislocations under the indentor altogether. On the other hand, an unduly small value of T OL results in excessive refinement and a prohibitive computational burden.
In this section we explore this trade-off by way of numerical testing, with a view to bracketing the optimal value of the tolerance parameter T OL. To this end, we fix the cluster size to r = √ 2 [σ], and we use the approximate cluster summation weights with a cutoff N c = 2000. The indentor is driven into the crystal at increments of 0.1 [σ] up to a maximum indentation depth of 2.0 [σ]. The first loading step does not result in the nucleation of dislocations, even in the fully atomistic model, and is identical to the calculations described in the preceding sections. By contrast, at the indentation depth of 2.0 [σ] the presence of dislocation structures under the indentor is clearly evident in the fully atomistic model, Fig. 12.
The energy error e is plotted in Fig. 13 as a function of indentation depth δ for different values of the tolerance parameter T OL. Initially, the error is insensitive to T OL, and all curves ostensibly follow the same linear growth pattern. However, the error curves diverge after a penetration depth δ = 0.5 [σ]. For large values of T OL, the error continues to  The final mesh at δ = 2.0 [σ] corresponding to a tolerance of T OL = 0.004 is shown in Fig. 14. The total number of representative atoms in this configuration is N h = 141, 469, which corresponds approximately to 1/8 of the total number N of atoms in the sample. It should be stressed, however, that since all of new representative atoms are inserted in the vicinity of the indentor, the ratio N h /N can be made arbitrarily small by considering increasingly large samples.
The dislocation structures predicted by the quasicontinuum method for a tolerance T OL = 0.004, Fig. 15, should be contrasted with those predicted by the fully atomistic model, Fig. 12. The figures display the energetic atoms (red) underneath the surface of the crystal (blue). Evidently, these structures differ in detail. This is expected, owing to the massive lack of uniqueness which characterizes the problem. Thus, as dislocation are nucleated and the lack of convexity of the energy comes into play, numerous deformation paths become available to the crystal. Many of these deformation paths differ only slightly in their energy content and are, therefore, indistinguishable at the macroscale,   e. g., in terms of the corresponding force-depth of penetration curves. This is specially so when the crystal is loaded along a direction of high symmetry, which results in the possible activation of a large number of competing slip systems. Under these conditions, a pointwise comparison of solutions becomes mathematically meaningless, and the sole meaningful criterion to measure the quality of a solution is its energy contents. By this criterion, the quasicontinuum and fully atomistic solutions are of comparable quality, as their energies differ by less than 1%. This contention is confirmed by Fig. 16, which compares the force on the indentor f as a function of penetration depth δ predicted by the fully atomistic and quasicontinuum (T OL = 0.004) models. Indeed, the two curves are ostensibly indistinguishable even beyond the onset of dislocation nucleation. On this basis, in practice choices of the tolerance T OL in the range of 0.001-0.004 would appear to strike an adequate balance between accuracy and performance demands.

Summary and Discussion
We have developed a streamlined and fully three-dimensional version of the quasicontinuum method of Tadmor et al. [18,19] and we have presented a numerical analysis of its accuracy and convergence characteristics. As a new addition to the theory, we have formulated a new class of summation rules in which the lattice function being summed is sampled over clusters of atoms. The size of these clusters is an adjustable parameter and controls the accuracy of the summation rule. We have also presented an efficient method for computing the requisite summation weights in linear time. Beyond its usefulness as a numerical scheme for approximating lattice sums, the cluster approach provides a systematic means of sampling the behavior of small representative crystallites, and thus opens a possible avenue for incorporating additional physics to quasicontinuum models such as diffusion of solute atoms and thermal lattice vibrations.
We have presented a suite of numerical tests which demonstrate the accuracy and performance of the method. As expected, the accuracy of cluster summation rules increases with cluster size. Furthermore, our numerical tests suggest that the addition of a single shell of neighbors suffices to stabilize the rank-deficiency which afflicts node-based summation rules. The computed convergence rate of the method in problems involving the application of point loads to crystals is close to linear when the lattice sums are performed exactly, and decreases somewhat when the sums are approximated using a cluster summation rule. It is worth noting that, contrary to continuum problems for which a well-developed approximation theory exists, a similar approximation theory for lattice problems appears to be missing at present, even for linear problems. The development of rigorous error bounds for finite-element approximations to lattice problems is a clear worthwhile direction of future research.
Perhaps the most interesting issue among those addressed in this paper concerns the ability of the quasicontinuum method to simulate microstructural evolution, which owes largely to mesh adaptivity. The same lack of convexity which allows for defects and microstructures to arise in the first place renders solutions massively nonunique. For a given loading or prescribed deformation, it is in general possible to find multiple equilibrium states of a crystal with defects possessing equal or nearly equal energies. These states are indistinguishable from a macroscopic point of view, e. g., in the sense of yielding identical force-displacement curves. Whether the solution follows one deformation path or another depends sensitively on small perturbations of the system, including details of the mesh design. Under these conditions, the pointwise comparison of solutions is not particularly meaningful. Our numerical tests suggest that, with sufficient mesh adaptivity, the quasicontinuum method is capable of simulating evolving microstructures comparable, in energetic terms, to those obtained from a full atomistic calculation.