Mixed finite element and atomistic formulation for complex crystals

A general formulation for the analysis of complex Bravais crystals using atomic energy functionals embedded within a ﬁnite element framework is presented. The method uses atomistic potentials to determine the constitutive response of the system. Unlike traditional ﬁnite element methods, the nonlinear elastic effects, the symmetries of the underlying crystal, and the possibility of uniform structural phase transformations are naturally included in this formulation. Explicit expressions for empirical energy functionals with separable two-and three-body potentials, and semiempirical tight-binding energy functionals with two-center integrals are presented. A simple application to silicon underscores the importance of including internal relaxation in a ﬁnite element treatment of a complex crystal. In a forthcoming companion paper, the method presented here is applied to the nanoindentation of silicon. (cid:64) S0163-1829 (cid:126) 99 (cid:33) 14301-7 (cid:35)


I. INTRODUCTION
Computer simulations are a valuable tool for understanding the fundamental processes that determine how a material responds to external loads.Broadly speaking, there are two classes of approaches that have proven particularly useful: atomistic and continuum simulations.Atomistic approaches provide detailed information about microscopic processes, and advances in computer design and algorithms have made possible the simulation of ever larger systems.However, even the most ambitious atomic simulations using empirical potentials can only handle systems sizes of order 10 9 atoms, which means that system dimensions generally cannot exceed a few hundred nanometers.Moreover, reliable results for large systems depend on the availability of high-quality, efficient energy functionals.Accurate, transferable empirical potentials are available for many pure metals and a few single-element covalent systems, but very few systems that contain more than one element.An alternative is to use quantum-mechanical energy functionals based on tightbinding or density-functional theory formulations.These are more transferable and can be applied to more complicated crystals, but are several orders of magnitude slower in computational speed than empirical potentials.
Continuum mechanics simulations have a different set of strengths.Since atomic details are coarse grained out of the formulation, much larger systems can be treated.But the constitutive equations that relate the strain and stress in a system are often empirically determined, making their reliability over a wide range of deformations dubious.In addition to this difficulty, the lack of an explicit connection between the continuum fields and atomic degrees of freedom makes it difficult to take into account the important properties of a crystal that depend on its atomic structure, such as crystal symmetries, invariance of a primitive cell with respect to shear by a full Burger's vector, and for complex crystals, the possibility of structural phase transformations.We seek to overcome these deficiencies in traditional continuum approaches by explicitly embedding atomistic totalenergy calculations within a finite element framework.
We focus on a recently introduced approach referred to as the quasicontinuum method. 1 In this approach a continuum finite element formulation is used to characterize the mechanical response of a given system.The difference from standard finite element methodologies lies in that the constitutive response of the system is obtained from an atomistic calculation rather than an empirical phenomenological rule.The basic idea is that every point in a continuum corresponds to a very large region on the atomic scale.Thus, the constitutive response at that point, i.e., the stress-strain relation there, may be obtained by deforming the underlying crystal structure by the local strain to obtain the local state of stress.This approximation will be adequate as long as the variation of the continuum fields is ''slow'' on the atomic scale.When this breaks down, such as near defect cores, the quasicontinuum method includes a nonlocal limit where the inhomogeneous deformation is specifically accounted for.Here, however, we focus on the local limit of the quasicontinuum formulation, leaving the treatment of the nonlocal defect effects to a future publication.
In its original version the quasicontinuum method was formulated with simple Bravais crystals in mind.Many useful materials ͑silicon being one example͒ have more than one atom per primitive unit cell.Furthermore, with a sufficiently large unit cell, a realistic representation of amorphous structures and alloys could be modeled, permitting finite element simulations of these complex materials.Extending the quasicontinuum method to simulate materials with multiple atoms in a unit cell is the focus of this paper.We present a formulation that maintains the advantages of traditional finite element methods, including correct treatment of the far-field boundary conditions.At the same time, the formulation naturally incorporates the anisotropy and symmetries of the material, all the nonlinearities associated with finite deformation, and the possible occurrence of structural phase transformations.Finally, the efficiency of the finite element formulation permits computationally expensive quantummechanical energy functionals to be applied to significantly larger systems than would be possible in atomistic simulations.
This generalized quasicontinuum method is ideal for simulating a variety of complex systems.In a forthcoming companion paper, the method is applied to nanoindentation of silicon. 2 Additional applications include modeling the constitutive response of ferroelectric materials, 3 simulating microelectromechanical systems devices, and modeling highly strained nickel-aluminum alloys.All these systems have multiple atoms in a unit cell, and internal relaxation plays an important role in determining the response of the system to external loads.Furthermore, in some systems large deformation can lead to phase transformations: silicon has been observed to transform into a conducting phase under an indenter, 4 and NiAl can undergo martensitic transformations near a crack tip. 5 The formulation of the quasicontinuum method presented here naturally captures these diverse physical phenomena.
The outline of the paper is as follows: Section II briefly describes the finite element framework.Section III focuses on the atomistic aspects of the formulation.First, the general expressions for the stresses and elastic constants for a general complex crystal are derived.Second, the explicit energy density expressions that link the continuum and atomic descriptions are presented for the case of a separable two-and three-body potential, and for the case of a two-center tightbinding Hamiltonian.Section IV shows how internal relaxation dramatically affects the behavior of system undergoing simple shear.

II. FINITE DEFORMATION CONTINUUM FRAMEWORK
The continuum formulation adopted here is that of finite deformation.Thus, no assumptions are made regarding the smallness of the strains.This can be important in microsystems where large mismatch and thermal strains as well as large mechanical operating strains are sometimes encountered.Within a finite deformation framework we differentiate between a solid in its reference or undistorted state ͑also referred to as the material frame͒ and its current or deformed state ͑called the spatial frame͒.Consider a crystal occupying a configuration B 0 in the reference frame which is mapped to its current shape B by the deformation mapping (X) ͑see Fig. 1͒.Thus every point X in the reference configuration is mapped to some point x in the current configuration by

xϭ͑X͒ϭXϩu͑X͒, ͑1͒
where u(X) is the displacement at point X.The deformation of an infinitesimal neighborhood of point X ͑represented as a shaded circle in the figure͒ is completely characterized by the linear part of (X).This defines an affine mapping,

dxϭFdX, ͑2͒
where F is the deformation gradient, where I is the identity tensor.In indicial notation, F iJ ϭ i,J and dx i ϭF iJ dX J where upper case indices refer to the material frame, lower case indices to the spatial frame, (•) ,J indicates differentiation with respect to X J and Einstein's summation convention on repeating indices is observed.
The deformation mapping and its gradient characterize the kinematics of the deformation.The constitutive nature of the material is introduced through the strain energy density function W which relates the energy at a point to the local state of deformation there.From the hypothesis of locality and use of the entropy production inequality it may be shown 6 that W can only be a function of F, i.e., WϭW(F) and not its derivatives or any other variables.Furthermore, by invoking the postulate of material frame indifference it can also be shown 6 that the dependence of W on F can only be through the right Cauchy-Green tensor CϭF T F. However, for the current formulation it is more convenient to work in terms of F.
We may now formulate the general boundary value problem.To this end, we partition the reference boundary ‫ץ‬B 0 into a displacement component ‫ץ‬B 0U and a traction component ‫ץ‬B 0T ͑see Fig. 1͒.The solid is subjected to prescribed displacements u on ‫ץ‬B 0U and prescribed tractions T on ‫ץ‬B 0T .Stable configurations of the crystal are identified with minimizers of the potential energy, where F is a function of u through Eq. ͑3͒, and the trial deformations u satisfy the essential boundary condition u ϭu on ‫ץ‬B 0U .
To solve the problem we discretize our continuum into a set of finite elements bounded by nodes ͑see Fig. 2͒.The continuous displacement field u is now approximated by finite element interpolation from its nodal values u a , where there are N n nodes in the mesh and N a (X) are the standard finite element interpolation, or ''shape'' functions. 7he potential energy in Eq. ͑4͒ may also be rewritten in the discretized form

͑6͒
where N e is the number of elements and ⍀ e is the domain of element e.All integrals in Eq. ͑6͒ are normally carried out by numerical quadrature, 7 e.g., where Q is the order of the quadrature rule, e q are the quadrature weights for element e and F e q is the deformation gradient evaluated at the quadrature point q of element e.In the limit of linear triangular elements ͑as employed in paper II͒, a single quadrature point is used and the volume integral simply reduces to ⍀ e W ˆ(F e ), where F e is the constant deformation gradient in the element.
To find the equilibrium configuration we need to minimize the potential energy ⌸ ˜(͕u i ͖) in Eq. ͑6͒.We can use either a conjugate gradient algorithm ͑see Ref. 8 for more detail͒ which requires the gradient of ⌸ ˜with respect to u a or quasi-Newton schemes ͑see Ref. 9͒ which have superior convergence properties but also require the second derivative matrix or Hessian of Eq. ͑6͒.Differentiating Eq. ͑6͒ we have where P is the first Piola-Kirchhoff stress tensor defined by Note that the first Piola-Kirchhoff stress is a mixed tensor, i.e., it has one index in the reference configuration and one in the spatial configuration and it is not symmetric.
The Hessian can similarly be obtained where D are the mixed tangential moduli, We now have expressions for the gradient and Hessian of the potential energy in Eqs.͑8͒ and ͑10͒ in terms of the stresses and elastic moduli within the elements.These depend on our choice of a strain energy density function W through an appropriate atomistic model.

III. ATOMISTIC CONSTITUTIVE LAW FOR COMPLEX BRAVAIS LATTICES
In order to formulate a constitutive law for a continuum from the atomic interactions of the discrete crystal that underlies it, we must hypothesize a connection between the continuum displacement field and the motion of atoms.The standard reasoning consistent with the locality approximation of continuum mechanics is that the atomic environment at a continuum point is characterized by the deformation gradient there; this is referred to as the Cauchy-Born hypothesis. 10Thus, each continuum point is taken to represent a large, essentially infinite, region on the atomic scale which is homogeneously distorted according to the deformation gradient at the point.The energy of the distorted crystal and its derivatives can then be obtained from any atomistic model of choice from approximate classical descriptions to rigorous first-principles methods through the use of periodic boundary conditions.
It is worth noting the limitations of the local approximation.First, a local formulation has no internal length scale, thus it is the strain in an element and not its absolute size that enters into the constitutive calculation.This approximation will be valid as long as the continuum displacement field is slowly varying on the atomic scale.If it is not, a nonlocal formulation will be necessary ͑see Ref. 11 for details͒.Second, surfaces and other interfaces will be invisible to a local formulation.It is possible within a local formulation to include regions with different constitutive descriptions, e.g., as in a polycrystal, and of course in a numerical simulation the model is finite and must terminate at surfaces.However, the energetics of the interfaces will not be accounted for since, as stated, the constitutive calculations at each point are performed for an infinite crystal free of boundaries.Third, the details of atomic scale defects such as dislocation cores will not be correctly reproduced.Interestingly, though, due to the nonconvexity of an atomistically calculated strain energy function, a local model can sustain stable dislocations with finite core energies.This is because the energy is computed from the atomic interactions in a crystal and thus all symmetries are automatically reproduced and more importantly lattice invariant shears are recognized.As a result, when a lattice is sheared by a full Burgers vector, as it is in the wake of a dislocation, the energy penalty is nil and the dislocation is stable ͑see Ref. 1 for details͒.The relaxed core structure, however, will only be a crude approximation to the correct structure.Now let us consider the application of the Cauchy-Born rule to complex Bravais lattices.We consider a crystal which may be represented as a Bravais lattice with Bravais vectors A k and N a ϩ1 basis atoms per Bravais site.The coordinates of the atoms making up the lattice in the reference coordinate system will be where N is the Bravais site number running from 0 to N s , if there are N s ϩ1 sites, and 0 corresponds to the origin site ͑i.e., l k 0 ϭ0),n is an index referring to the basis atom number, l N is a triplet of integers locating Bravais site N in space and b n is the position of basis atom n relative to the Bravais site.In order to simplify the notation we bundle the site and basis atom number into a single index denoted by a Greek letter, e.g., ϭ(N,n), which completely identifies an atom.Without loss of generality, we may set b 0 ϭ0 so that one basis atom always lies at the Bravais site.Applying the Cauchy-Born rule, the positions of the atoms after deformation is given by where n are additional inner displacements resulting from energy relaxation with respect to the basis atom positions ͑see Refs.12-14 for details͒.Again, without loss of generality and in order to rule out rigid-body translation we may fix the lattice by setting 0 ϭ0.By defining the inner displacements n in the reference configuration we automatically guarantee the invariance of these measures with respect to rigid-body rotation. 12 general, the energy of a collection of atoms ͕x ͖ will be a function of their coordinates, EϭE͕͑x ͖͒.

͑14͒
The strain energy density W follows by normalizing this quantity with respect to the atomic volume in the reference configuration ⍀,

͑15͒
We stated earlier that the strain energy function must obey material frame indifference, i.e., invariance with respect to rigid-body rotation.Martin 13 studied the implication of this to atomic energy functions and concluded that the energy can only depend on the atomic coordinates through scalar products of their relative positions.We define R and r as relative position vectors between pairs of atoms in the reference and deformed configurations, respectively, here ϭ(M ,m) and From Eq. ͑17͒ it is clear that scalar dot products of r will depend on F through F T F or C as expected from material frame indifference.For example, the magnitude of r will be given by Other scalar products will have a similar dependence.Thus, in general, we have WϭW(F, 1 ,•••, N a ).For a given macroscopic deformation, the strain energy must be minimized with respect to the inner displacements, 14 thus we require and ‫ץ‬ 2 W/‫ץ‬ m ‫ץ‬ n must be positive definite.The relations in Eq. ͑19͒ form a system of 3N a equations for the 3N a unknowns J n .The dependence of W on n will generally be highly complex and the resulting system of equations coupled and nonlinear.Practically, the solution of Eq. ͑19͒ is best achieved by employing a numerical minimization algorithm, such as the conjugate gradient method, to minimize W with respect to the inner displacements with the deformation gradient F held fixed.This minimization can be strongly history dependent at zero temperature ͑i.e., the local minimum to which it converged will depend on the initial guess͒.Thus, it is necessary to store the inner displacements n at each quadrature point, using the previous solution at that point as the initial guess for the next iteration.
After minimization, we obtain the equilibrium inner displacements ˆnϭ ˆn(F), which implies that the strain energy is only a function of F, WϭW͑F, ˆn͒ϭW"͑F, ˆn͑F͒…ϵW ˆ͑F͒.͑20͒ We may now compute the stress and moduli tensors in the presence of internal relaxation.The first Piola-Kirchhoff stress was defined earlier in Eq. ͑9͒.For a complex Bravais lattice we have

͑21͒
where the subscript ˆindicates the expression is evaluated at the relaxed inner displacements ˆn and repeat indices on vector quantities (K and n on ) implies summation over all components.Under ideal conditions the second term in Eq. ͑21͒ drops out due to Eq. ͑19͒ and we are left with the standard definition for the stress tensor evaluated at the relaxed inner displacements

͑22͒
However, since in practical applications Eq. ͑19͒ is only satisfied approximately to some given precision, it is preferable to work with the more complicated expression in Eq. ͑21͒.If Eq. ͑22͒ is used instead, the energy and its stress derivative would not be compatible and numerical difficulties will ensue at the final stages of convergence where the minimum of the energy surface would not correspond to a stress-free state.To use Eq.͑21͒ we need to compute the nontrivial derivative ‫ץ‬ n /‫ץ‬F.To this end we differentiate Eq. ͑19͒ with respect to F kL and obtain Inverting this relation we find the required derivative To obtain the mixed tangential moduli tensor D in Eq. ͑11͒ we must differentiate Eq. ͑21͒ with respect to F, explicitly accounting for the dependence of ˆn on it.Differentiating, rearranging and simplifying by making use of Eq. ͑19͒, we find where for notational simplicity we have dropped the indication that all terms are evaluated at ˆn, but this is implied.The last two terms in Eq. ͑25͒ can be shown to cancel when Eq. ͑23͒ is substituted into the last term and we are left with the following simpler relation: Finally, by substituting Eq. ͑24͒ into Eq.͑26͒ we can obtain a more symmetric expression for the elastic moduli It is worth noting that the correction due to internal relaxation will be non-zero and may be significant even for an undistorted crystal ͑i.e., FϭI).Thus, internal relaxation may not be negligible even for infinitesimal straining in some cases; this is the case for silicon.As in the case of the stresses, an argument may be made that due to numerical precision considerations the full expression for the elastic moduli which does not assume Eq. ͑19͒ ͑not given here͒ should be used in lieu of the simplified expression in Eq. ͑27͒.However, since the elastic moduli are only used to obtain search directions and are not used as a convergence criterion the discrepancy is less significant and Eq.͑27͒ may be used.

Two-body and three-body separable potentials
To this point the derivation was general without reference to a particular energy function.We now assume that the energy of a collection of identical atoms may be expressed in the form of a two-body and three-body separable potential,

͑28͒
where r is the distance between atoms and and is the angle defined by the triplet of atoms , , and .In Eq. ͑28͒ a sum over a Greek index implies a double summation over the Bravais sites and basis atoms associated with the index.In this case the indices , , and run over all Bravais sites and associated basis atoms in the solid.Note that in Eq. ͑28͒ we have assumed a trigonometric dependence of the energy on the angle .This is normally the case and greatly simplifies the expressions derived later in Appendix A.
For a homogeneously distorted crystal, the same basis atoms at different lattice sites will experience identical environments and thus the calculation of Eq. ͑28͒ may be limited to a single unit cell in order to obtain the correct energy density.For convenience we choose to compute the energy with respect to the origin Bravais site oϭ(0,a), Wϭ  1

͑29͒
where ⍀ is the primitive Bravais cell volume, r o is the distance of atom from basis atom a of the origin Bravais site, and o is the angle measured relative to that atom.The ͚ ˆsymbol indicates summation only over the second component of a Greek index ͑i.e., over a in this case͒.
Given the specific potential description of Eq. ͑29͒ we may explicitly obtain the derivatives of the strain energy density appearing in the stress ͑21͒ and moduli ͑27͒ expressions.The first derivatives are given

͑31͒
In the above expressions ( ) , indicates differentiation with respect to r o and ( ) , indicates differentiation with respect to cos o .The second derivatives can easily be obtained in similar fashion and are not given here for reasons of brevity.The derivatives in Eqs.͑30͒ and ͑31͒ and the second derivatives contain partial derivatives of the form ‫,ץ/‪r‬ץ,‪F‬ץ/‪r‬ץ‬ When the energy is computed classically as in Eq. ͑29͒, it is necessary to locate all neighboring atoms displaced to within the computed atom's cutoff radius upon deformation.However, obtaining an atom's neighbors in a deformed crystal is computationally intensive.A convenient approach, is to construct a crystal template.This is a list of the atomic coordinates of all atoms within some large radius of the unit cell sorted by distance from the basis atoms.It is then possible to compute a conservative estimate for the outer radius from which atoms could move inside the cutoff sphere.This is referred to as an influence radius.A derivation of this measure for complex Bravais lattices is given in Appendix C.
In simulations, this influence radius serves two purposes.First, it is used to ensure that the crystal template stored for the energy calculations is sufficiently large.Second, the influence radius is used to identify the smallest subset of atoms which have a possibility of entering into the cutoff radius, thereby making the energy calculations as efficient as possible.

Tight-binding models
A completely different method for computing the energy of a collection of atoms uses the tight-binding method, a simple quantum-mechanical model. 15In this method the atoms are treated as classical particles that interact in part through an effective potential exerted by the electrons that are treated quantum mechanically.Hypothetical basis orbitals with the angular symmetries of single atom eigenstates are centered around each atom.The details of the basis functions do not enter into the energy calculation, but only the interactions between basis elements ͉ i ͘ that form the over- lap matrix S i j ϭ͗ i ͉ j ͘, ͑32͒ and the Hamiltonian matrix

͑33͒
The indices i and j run over the all the atoms and all the basis elements centered on each atom.The energy eigenvalues are found by solving the generalized eigenvalue equation where ⑀ n is the nth eigenvalue and n j is the jth component of the nth eigenvector.
In the particular tight-binding model used here, the energy density, W, is made up of a classical two-body pair repulsion term and a band-structure term 16 Wϭ 1 where N k is the number of k points, and the sum ͚ n is over both bands and k points.The distribution f n is taken to be the Fermi function, where ⑀ f is the Fermi energy, and is a smoothing parameter.The smoothing parameter is used to maintain numerical stability when configurations with no band gap are treated with a finite k point mesh in the course of a simulation.It is not related to the physical temperature of the system.Ideally simulation results should not depend on the value of .Da Vita 17 has shown that while the energy and freeenergy density both vary as 2 , the average of the two varies as 4 ͑see Appendix B͒.Thus for the tight-binding model we use an average of the energy and free-energy density, denoted W ˜, which has an entropic term added to the expression in Eq. ͑35͒, For the derivatives,

͑38͒
where With the aid of the Hellman-Feynman theorem 18,19 one can show, ͑40͒ where r i o is the ith component of the the position vector between the atoms o and .

IV. SIGNIFICANCE OF INTERNAL RELAXATION
So far the derivation has been given in general terms for arbitrary crystal structures.To demonstrate the benefits and implications of implementing an atomistically based constitutive law we turn to the specific example of silicon.In particular, we will investigate the effect of internal atomic relaxation on observable macroscopic response.We compare two different atomistic descriptions, the classical Stillinger-Weber ͑SW͒ model 20 and a nonorthogonal tight-binding Hamiltonian ͑TB͒ developed by Bernstein and Kaxiras. 16e focus on a simple shear deformation.We choose a ͓1 ¯10͔ shear direction with a ͓111͔ slip plane normal.The corresponding deformation gradient has the form

FϭIϩ␥s n, ͑41͒
where ␥ is the shear parameter, s is the slip direction, and n is the slip plane normal.Figure 3 presents the strain energy density of the crystal as it is sheared.In frame ͑a͒ no internal relaxation is allowed, while in frame ͑b͒ relaxation is taken into account.The difference in the energies are startling.Allowing the second basis atom to relax reduces the maximum energy by a factor of nearly 30.This is emphasized in the figure by the dashed inset in frame ͑a͒ which corresponds to the domain of frame ͑b͒.
In addition to the dramatic reduction in energy, and perhaps more importantly, the periodic nature of the response is also affected by relaxation.Whereas in the unrelaxed case the perfect crystal structure is restored after shearing by ␥ ϭ2ͱ6, when the second basis atom is allowed to relax, the period is reduced by a factor of four, i.e., ␥ϭͱ6/2.The reason for this becomes apparent by studying a single tetrahedron as it is sheared with and without internal relaxation as seen in Fig. 4. The lighter atoms are the Bravais site atoms that are constrained to move according to the applied deformation field.The dark atom is the second basis atom which in frame ͑a͒ is constrained along with the Bravais site atoms and in ͑b͒ is free to move relative to them.Each frame contains a series of three snapshots labeled by letters as indicated in Fig. 3.
In the unrelaxed case, we see that the second basis atom lags further and further behind the triplet of atoms above it.The tetrahedron gets more and more distorted as the crystal is sheared.At point C the crystal structure is far from the ideal diamond structure with a correspondingly high energy.The perfect structure will be restored only by shearing four times further.When relaxation is taken into account the second basis atom displaces to minimize the distortion.The initial structure in D is, of course, the same as A in the unrelaxed case, but then as the upper triangle moves to the right, the second basis atom moves with it, remaining underneath.In this manner, the bond angles between the second basis atom and the upper triplet of atoms, remain closer to the ideal 109.47°.As the shear progresses, the lower left atom falls further behind causing the second basis atom to move down in the ͓1 ¯1 ¯1 ¯͔ direction and slightly back in the ͓11 ¯0͔ direction to reduce the strain on that bond ͑see snapshot E͒.As the lower right atom draws closer, the second basis atom is pulled towards it, now moving a little ahead of the upper triangle, and eventually the perfect crystal structure is restored in F.
Finally, it is of interest to compare the predictions of the tight-binding model with the Stillinger-Weber potential.The tight-binding approach is attractive because its quantummechanics basis promises for a more transferable model, however it is computationally far more intensive than classical approaches ͓e.g., computing the relaxed curves in Fig. 3͑b͒ takes 300 times longer using the TB method than with the SW potential͔.For the simple shear deformation studied here both approaches yield rather similar results.The overall energetics are close and the structures observed along the shear trajectory were nearly identical.The main concern associated with the SW potential is the appearance of a nonphysical local minimum corresponding to the midshear position ͓see Fig. 3͑b͔͒.This appears to be a general feature of classical potentials.In fact, it appears that the more elaborate the potential and the more it is fitted to given structures, the more undesirable local minima appear between these desired states.This may not be of concern in finite temperature simulations where the system has sufficient thermal energy to escape the local wells, however at zero temperature these minima can be problematic.Thus, in addition to transferability, the smoothness of the energy surface of more rigorous quantum-mechanical approaches is of significant benefit.These issues will be pursued in more detail in a companion paper where more complex deformation pathways will be studied. 2

V. SUMMARY AND CONCLUSIONS
We have formulated a finite element method for materials with multiple atoms in a unit cell that uses atomistic energy functionals to determine the constitutive relations.The method presented here naturally incorporates crystal anisot-FIG.3. Strain energy density of a silicon crystal as it is sheared along the ͓1 ¯10͔/(111) slip system.Frame ͑a͒ presents the case where no internal relaxation is allowed.In frame ͑b͒, which occupies the dashed inset in ͑a͒, relaxation is taken into account.Solid lines were computed using the TB formalism, dashed lines with the SW potential.FIG. 4. Shearing of a single silicon tetrahedron ͑a͒ without internal relaxation and ͑b͒ with internal relaxation.The dark atom is the second basis atom which in frame ͑b͒ is allowed to relax at each step.The letters correspond to points on the graphs in Fig. 3.All structures were computed using the TB formalism.
ropy, nonlinear contributions to the energy density, and the possibility of structural phase transitions within the finite element framework.We have derived explicit expressions for an empirical potential with two-and three-body terms, and for a two-center tight-binding Hamiltonian.Throughout, special care has be taken to correctly treat the internal degrees of freedom.We have shown that correct treatment of these degrees of freedom dramatically affects the behavior of even a simple system such as sheared bulk silicon.

͑A4͒
The second derivatives are given by

͑A10͒ APPENDIX B: AVERAGED ENERGY DEPENDENCE ON SMOOTHING PARAMETER
In this section we show that in a tight-binding-model the average of the energy and free-energy density W ˜and the corresponding forces vary as 4 .We assume that we have a smooth, continuous density of states, g(⑀).Thus our starting expression for W ˜is the same as the expression in Eq. ͑37͒, except that (2/⍀N k ) ͚ n is replaced by ͐g(⑀)d⑀.We consider W ˜to be a function of and ⑀ f , and require that the number density of electrons N be independent of the smoothing parameter.Then, Determining the zero-order contributions to the derivatives with respect to the Fermi energy is immediate, and additional terms in the expansion are of order 2 , The derivatives with respect to can be evaluated using Sommerfield's lemma 21 Combining Eqs.͑B1͒-͑B5͒, we obtain Taking the derivative of the expansion of W ˜with respect to atom positions confirms that the forces have the same temperature dependence as W ˜.
In numerical tests of typical highly compressed silicon unit cells, we find that W ˜is significantly less sensitive to the smoothing parameter value than the energy or free energy alone.Up to values of a few tenths of an eV for ,W ˜deviates from the zero-temperature energy by about 0.001%, while the energy and free energy deviate by about 0.01%.However, the behavior of the forces is worse by at least an order of magnitude.We keep the temperature as low as possible during simulations, between 10 Ϫ2 and 10 Ϫ3 eV.

APPENDIX C: INFLUENCE RADIUS FOR COMPLEX LATTICES
The energy and derivative computations in the main body of the paper are computationally intensive.It is of great practical importance to minimize the number of atoms that are involved in the computation.To this end we compute an influence radius. 1 This is the distance of the furthest atom from the origin that will come within the cutoff radii of the computed Bravais site atoms.We obtain a continuum estimate for this measure.
Following the notation of Sec.III, a generic point X in the reference configuration is mapped to the deformed configuration by

xϭFXϩz, ͑C1͒
where zϭF is the inner relaxation in the deformed configuration.We wish to maximize X T X subject to the constraint x T xϭr c 2 , where r c is the cutoff radius.We define the functional where is a Lagrange multiplier.Inverting Eq. ͑C1͒ and substituting into Eq.͑C2͒, we have where B Ϫ1 ϭF ϪT F Ϫ1 is the inverse of the left Cauchy-Green deformation tensor.Taking a variation of Eq. ͑C3͒ with respect to x, ␦⌳ϭ0, leads to the system of equations ͑ B Ϫ1 ϪI͒xϭB Ϫ1 z. ͑C4͒ In the absence of internal atoms, i.e., zϭ0, this is a simple eigenvalue problem with the solution R inf ϭr c ͱ max , where PRB 59 243 MIXED FINITE ELEMENT AND ATOMISTIC . . .max is the maximum eigenvalue of B Ϫ1 .When z is not zero, Eq. ͑C4͒ is no longer an eigenvalue problem.Instead we need to solve for the value of by solving Eq. ͑C4͒, subject to the constraint x T xϭr c 2 .To facilitate this, we introduce the eigenvalues and eigenvectors of matrix B Ϫ1 , i and i , respectively, and express x and B Ϫ1 z in this basis Substituting Eqs.͑C5͒ and ͑C6͒ into the system ͑C4͒ and making use of the eigenvalue identity, we obtain We can then solve for ␣ i and thus obtain x from Eq. ͑C5͒, The coefficients ␤ i in Eq. ͑C10͒ can readily be obtained from Eq. ͑C6͒.Applying j to both sides of Eq. ͑C6͒, making use of the symmetry of B Ϫ1 and the orthogonality of its eigenvectors and recalling relation ͑C7͒, we find ␤ i ϭ i ͑z• i ͒͑no sum͒.͑C11͒ The remaining unknowns in Eq. ͑C10͒ are the values of the Lagrange multiplier .We obtain an implicit equation for the values of by enforcing the constraint x T xϭr c 2 , where we have again made use of the orthogonality of the eigenvectors i .Equation ͑C12͒ may have as many as six distinct roots.These will appear in pairs around the eigenvalues of B Ϫ1 , i , where the left-hand side of Eq. ͑C12͒ is infinite.This can be seen in Fig. 5 where a typical curve for silicon is drawn.
Once the roots k are found, we may compute the corresponding influence radii where x k is the solution ͑C10͒ with ϭ k .After substituting in the appropriate expressions and following some algebra we obtain The influence radius R inf follows as

͑C15͒
Equations ͑C12͒ and ͑C14͒ above represent a complete solution to the problem.However, solving for the roots k and obtaining the eigenvalues and eigenvectors of B Ϫ1 may prove to be too time consuming for this to be worthwhile.It is thus of interest to obtain an approximate solution to the problem that will not underestimate the exact influence radius, and will be significantly faster.
Empirically, we observe that the maximum influence radius is always associated with the maximum root, i.e., the one to the right of the maximum eigenvalue 1 ͑see Fig. 5͒.We can approximately compute max by assuming that the term associated with 1 in Eq. ͑C12͒ dominates ͓since the denominator ( 1 Ϫ max ) will be very small͔, thus Substituting Eq. ͑C16͒ into Eq.͑C14͒, we have

͑C17͒
The term z• 1 is bounded by zϭ͉z͉, since 1 is a unit vector.
In addition we may assume the jϭ1 term in Eq. ͑C17͒ dominates since the denominator is smallest and numerator greater there.Applying these assumptions and rearranging we find the final simplified estimate for the influence radius R inf ϭͱ 1 r c ͩ 1ϩ z r c ͪ .

͑C18͒
Interestingly, this expression is simply (1ϩz/r c )R inf 0 , where R inf 0 is the influence radius for the shuffle-free case.Although, we have not established it rigorously, in practice this estimate has always behaved as a conservative bound on the actual influence radius.*Present address: Faculty of Mechanical Engineering, Technion- Israel Institute of Technology, 32000 Haifa, Israel.

FIG. 1 .
FIG. 1.The deformation mapping relating the material (B 0 ) and spatial ͑B͒ configurations.Displacements are specified on the dark dashed part of the boundary (‫ץ‬B 0U ) and tractions on the remaining boundary (‫ץ‬B 0T ).

FIG. 2 .
FIG. 2. Finite element mesh of reference solid B 0 showing boundary conditions and nodal degrees of freedom.

FIG. 5 .
FIG. 5. A typical solution for the values of the Lagrange multiplier .The curve is drawn for Stillinger-Weber silicon (r c ϭ3.77 Å͒ for a randomly selected deformation gradient and inner displacement.The six roots k are indicated by black dots.