# A deformation gradient tensor and strain tensors for atomistic simulations

## Abstract

A kinematical algorithm is proposed for the construction of strain tensors from atomistic simulation data. Local strain tensors such as the Almansi and Green strain tensors suitable for use in large deformation molecular dynamics/statics simulations are computed directly from a discrete form of the deformation gradient. The discrete, incremental form of the deformation gradient emerges from a weighted least squares minimization that includes a length scale relating the distance from the atom in question with a particular radius. This region defines the nonlocal domain of the strain at that atom. The local atomic strain tensors are then computed using continuum definitions in terms of the deformation gradient. The results of molecular dynamics simulations are presented that compare the Almansi and Green strain tensors under inhomogeneous deformation and indicate that the small-strain approximation should not be used to determine large atomic strains.

Authors: Philip M. Gullett, Mark F. Horstemeyer, M I Baskes, H Fang

## Strain formulation

In this section we briefly consider the continuum deformation gradient, introduce a discrete version of the deformation gradient, and finally present the strain tensors.

### Weight function

In this paper, weighting factors are computed as follows: neighbors are grouped according to their proximity to the current atom, i.e. first nearest neighbors, second nearest and so forth. Next, weights are assigned to the groups using the cubic spline

$r=(r_{gi}-r_{g1})/r_{cut}$

Here, $r_{g1}$ and $r_{gi}$ are the distances between the current atom and the closest neighbors in the first and ith neighbor groups; $r_{cut}$ is the cutoff radius that specifies the domain over which the weight function is nonzero. The effect of changing the size of the cutoff radius is examined in section 2.

## Example strain calculation

This simulation provides an example of the computed strain for a block of nickel atoms subjected to fixed-end, simple shear-like boundary conditions. The embedded atom method (EAM) potential described by [13, 14] was used to compute atomic motions. The EAM utilizes a pair potential augmented by an effective embedding energy term to account for many-body interactions.

The geometry of the nickel atom block, shown in figure 3(a), is 35Å high, 140Å long and 14Å thick. The lattice crystal directions [1 0 0], [0 1 1] and [0 1 -1] correspond to the x, y, and z axes, respectively. The atomic velocities were initialized using a Boltzmann distribution at 300 K, and the system was equilibrated before the shear velocity was applied.

The simple shear-like deformation was achieved by imposing a constant x-direction velocity of vx = 0.035Å/ps on the +y surface atoms, while holding the −y surface atoms fixed as indicated in figure 3(b). The x and z surfaces were modeled with periodic boundary conditions. The imposed boundary conditions differ from a rigorous continuum description of simple shear in that the x faces are not constrained to remain parallel during the simulation, and the prescribed shear velocity is only applied to the atoms on the +y surface. As will be seen, the equivalence of the actual boundary conditions with those of true simple shear is approximately satisfied until the onset of plastic deformation. A block subjected to a true simple shear motion is described by a spatially constant deformation gradient

### Bulk strain response

The average computed shear strain results are summarized in figure 4. Qualitatively, figure 4(a) shows that the average computed Green and Almansi shear strains initially increase at the same rate as the assumed simple shear response until the maximum stress is reached at approximately 6.4% applied strain. After 6.4% applied strain, plasticity initiates via dislocation nucleation and motion, and the average computed Green and Almansi strains begin to differ. The computed Almansi strain increases at a rate that is slightly higher than the applied strain, and the Green strain increases a slightly lower rate.

Figure 4(b) shows that during the initial, elastic deformation, the difference between the computed Almansi and Green strains is negligible, averaging 0.23%. When the system reaches its maximum stress value at 6.4% applied strain, the difference in the Almansi and Green strains jumps to 5.0% and then increases steadily. After 10% applied strain, the difference between the computed Green and Almansi strains continues to increase, but at a lower rate.

The initial deformation of the system is consistent with the Cauchy–Born hypothesis with the motion of the boundary atoms being accommodated through stretching of atomic bonds and the uniform tilting of the initially vertical atomic planes as shown in figures 5(a) and (b). As expected, before the onset of plastic deformation, the computed Green and Almansi shear strains are consistent with a homogeneous deformation. At 3.0% applied strain, figure 5, the average local computed shears are 3.0% with a standard deviation of 0.6%. As the system approaches its peak stress, the variation in local computed strain increases. At 6.0% applied strain, figure 6, the bulk average computed shear strain is 6.0% with a standard deviation of 0.8%.

The system reaches the peak stress at 6.4% applied strain, and although bulk average computed shear strain is equal to the applied strain, the system undergoes a significant redistribution of deformation. Figure 7 shows that the shear strain is highly localized in the (1 0 0) slip planes. Near the dislocation nucleation sites, the strain levels vary from 0% to 20%. Away from the highly strained region, the lattice has uniformly unloaded to approximately 3% strain. This deformation pattern represents a transition from the homogeneous simple shear motion to an inhomogeneous one.

At 7.0% applied strain (figure 8) the deformation is characterized by organization and motion of dislocations concentrated along the (1 0 0), (1 1 1), and (-1 1 1) crystallographic slip planes. The computed Almansi and Green shear strains provide equivalent results for (1 0 0) slip but differ for (1 1 1) and (-1 1 1). This difference is consistent with the expected behavior in the post-plastic flow regime. In particular, the Almansi strain shows most shear activity along the (-1 1 1) slip plane, while the Green strain shows most activity along the (1 1 1) planes. This difference is accounted for by considering the orientation of the active slip systems. For example, the relative motion along adjacent (1 1 1) slip planes can be approximated locally by simple shear.

### Local atomic strains

To examine the local strain computations at individual atoms, we compare the computed results for atoms 764 and 2129 with positions indicated in figure 9. Atom 764 is located in a region below the active slip plains, while atom 2129 is subject to significant slip activity. The computed shear strain time histories for these atoms are shown in figures 9(b) and 9(c). Until the localization of strain occurs at 6.4% applied strain, the magnitudes of the computed strains match those of the the applied strain. After localization the computed strain rates correspond with the applied strain rate, but the magnitudes do not.

Atom 764 lies in a region of homogeneous and elastic deformation. The computed strain history follows a cyclic loading and unloading. During the loading cycles, the strain increases at a rate consistent with the applied strain rate. The local strain is then partially or completely relieved at a much higher rate due to localization of the strain along the slip planes. Atom 2129 is located in a region of significant inhomogeneous deformation, and the strain history follows a somewhat different pattern from that of atom 764. Instead of relaxing at the maximum stress, the strain jumps from 7.6% to 17.6%, followed by a linear increase in strain at a rate equal to the applied strain rate. At 8% applied strain, the strain relaxes from 18% to 14%. This decrease is then followed by a steady increase in computed strain and a jump at 14% applied strain from 21% to 25%. One observes that the jumps in strain at both atoms occur at the same applied strains, and that the strain at the atom in the elastic region decreased while the strain at the atom in the slip plane increased.

### Weight function effects on computed strain

The relevant length scale for the atomic strain examined by changing both the weight function form and cutoff radius. The bulk average computed shear strains using three different weight functions are plotted in figure 10. The weight functions used in computing the atomic Green strains were: a Heaviside step function in which all atoms within a cutoff radius were given a weight of 1.0, a cubic spline as shown in figure 2(a), and the multi-step function shown in 2(b). When considering the weight function form, the cutoff radius was fixed at 8Å resulting in a total of approximately 200 atoms contributing to the deformation gradient optimization. When examining the cutoff radius the size varied from 4 to 10Å.

During homogeneous deformation, the weight function has little effect on the computed strains with specific values varying from the applied strain by less than 1.0%. Once plastic flow begins, the magnitude of the computed strains exhibits some sensitivity, with the results using the Heaviside and cubic spline functions being closest to the applied strain.

Strains computed using the Heaviside function follow the applied strain more closely than those computed from the cubic or step function. After a jump from 1% to 2% at the transition from homogeneous to inhomogeneous deformation, the difference from the applied strain decreases to approximately 1.2%. The strains computed using the cubic and step functions, which are more sensitive to motion of the nearest neighbors, showa small but more pronounced deviation from the applied strain. The cubic jumps from 0.7% to 2.0% difference and increases to approximately 3.0%, and the step function jumps from 0.4% to 3.8% and increases to 4.0% over the course of the deformation.

The choice of weight function impacts the strain tensor calculation primarily through increasing or decreasing contribution of remote neighbors. For the Heaviside function, this results in a suppression of the effects of local slip because atoms within the 8Åcutoff radius but outside the 4Å wide slip bands make significant contributions to the local strain calculations, thus, yielding a result closer to the homogeneous applied strain. With the multi-step function, the primary contributions to the formation of the deformation gradient come the first and second nearest neighbors. In this case, the effect of shear along the crystallographic slip planes is locally pronounced and manifest through bulk Green shear strains that are smaller than the applied shear strains.

To directly examine the effect of weight function radius, the Almansi shear strains were computed for atom 2129 using the functions shown in figure 2(b), and the results are plotted in figure 11. By increasing the cutoff radius, additional atoms are considered in determination of the local strain. With a cutoff radius of 4Å, 12 neighbors contribute to the strain. A cutoff of 6Å adds a contribution of six additional neighbors, and so on.

The range of the weight function cutoff radii considered here includes the 5.65Å cutoff radius used for the energy potential. Varying the size of the strain cutoff radius independently from the energy cutoff radius, facilitated examination of the relevant length scale for calculation of strain. In this example, figure 10, the local variations in the computed strain arise due to localized slip. While differences arise in various locations, volume average Green strain shows minimal variation. Volume average Almansi strains may exhibit significant differences for 4 and 6Å cutoff radii and minimal variation for 8 and 10Å radii.

As expected, computed strains were independent of cutoff radius during homogeneous deformation. After 6.4% applied strain, the computed strains follow a similar pattern, but the strain magnitudes vary with cutoff radius size. At 6.4% applied strain, the computed strains all exhibit a jump, with magnitudes ranging from 6.6% for a cutoff of 4Å to 12.9% for a radius of 6Å. At approximately 8.4% applied strain, the computed strains relax. The strain computed considering only nearest neighbors (4 Å) decreases to within 2% of the applied shear curve, while the others indicate slightly higher shear levels. The next jump in the local strain occurs at 14% applied strain. Here the 4Å strain jumps to the same value as the 8Å measure, while the 6Å drops to the applied strain level. The largest jump occurs for the 10Å case which approaches a level of 30%. With this cutoff radius, significant contributions are made from atoms as far as 6Å which includes atoms at the bottom of the sample that are subject to the largest magnitudes of slip. Although the results at atom 2129 are sensitive to cutoff radii in this range, this effect is primarily due to the relative size of the slip region compared with the cutoff radius and the complexity and variation of the slip directions and magnitudes in the region. For regions in which slip is more uniform in magnitude and direction among neighbors, e.g. at atom 764, the variation of computed strain does not exhibit similar jumps.

## Conclusions

A new approach for computing a deformation gradient from discrete atomic data has been presented. The motivation for this algorithm was to develop a means of computing finite deformation strain tensors from discrete atomic simulation data that correspond to continuum strain tensors. The proposed deformation gradient emerges from a weighed least squares optimization of local deformation data and may be used to compute strain or other deformationrelated quantities.

The computations for the proposed formulation utilize position vectors and neighbor lists from atomistic simulations to compute a deformation gradient. The computational complexity of the algorithm, excluding the generation of neighbor lists, scales linearly with the number of atoms, i.e. O(N). Once the deformation gradient is constructed, any desired strain tensor may be computed. In contrast, strain tensors relying on geometric decompositions of interatomic regions, such as a Delaunay triangulation of atomic sites in R3, will require approximately O(N log N)

The simulation of simple shearing of an atomic slab was performed to evaluate the computed strain values of the proposed method. The results demonstrate the consistency between the computed atomic strains and continuum finite deformation strain predictions. For homogeneous deformations, the strain values were shown to be computed precisely, and rotations properly accounted for. Although other methods for computing strain are available, the advantage of the proposed method is that finite strains and rotations can be resolved correctly when based on the deformation gradient. This is not true for methods based on smallstrain definitions. In those methods even relatively small rotations will result in incorrect strain values. Thus, the finite strain model is better suited for establishing cause-effect relations for incorporation into macroscale models.

The computed shear strains in regions of localization are sensitive to the weighting of distant neighbors through either the choice of weight function and/or the size of cutoff radius used. The effect of localization on the computed bulk strains becomes less apparent as the effective size of the cutoff zone is increased. The sensitivity of the local atomic strain to the size of the cutoff radius appears to decrease when the size of the localization increases relative to the local atomic region.

## Acknowledgments

This work was sponsored by Mississippi State University Center for Advanced Vehicular Systems and by the US Department of Energy under Contract No W-7405-ENG-36.