Interpolation of Microscale Stress and Strain Fields Based on Mechanical Models

Multiscale methods have been attracting more and more attention for the last two decades. One of the main applications of the multiscale methods is to study material with defects [20, 22, 11, 1]. The region containing defects cannot be accurately described by the traditional continuum mechanics and atomic models have been used as a powerful tool to study the defects at small scale. However the high computational cost makes it difficult, if not impossible, to use a fully atomic model for problems of sizes large enough to be of any practical applications, and therefore multiscale methods which can couple the molecular dynamics and the continuum mechanics are called for. Good review papers on such methods can be found at [9, 19, 18, 28, 3]. Among the existing results, we consider the Quasicontinuum method [26, 27, 17] as one of the most systematically developed multiscale framework. It introduces a physically consistent coupling between the continuum mechanics and the molecular dynamics. Though showing promising results in the material defect analysis at nanoscale [26, 25, 16, 7, 5], there are still some problems to fix and many extensions to make, and one of them is to develop a method to interpolate the stress and strain fields for the microscale region of Quasicontinuum models. Stress and strain are only defined in continuum mechanics, therefore, without extra interpolation techniques, one cannot obtain the stress and strain fields for the microscale region which is modeled by molecular dynamics. On the other hand, the microscale stress fields are also found useful for studying the material behaviors at the microscale [8, 13, 21]. Moreover, the microscale stress itself has recently been exploited as part of the coupling variables [10] for obtaining smoother transition at the macro-micro boundary. Therefore, it is desirable to develop techniques that can interpolate the stress and strain fields for the particle methods which are consistent with their counterparts in continuum mechanics.

Various methods for obtaining the microscale stress tensor have already been available in literature. Among them, the so called Virial stress [14] provides simple expressions for the atomics stress and is still widely used in the literature, despite the controversial debate over its validity [30, 24, 6]. An appropriate version was put forward by Lutsko [15] and later improved by Cormier [4], who derived it in a more general form. However, the derivation process is found to be quite complicated. Fourier transformations and inverse transformations have been applied to solve the differential equations [15, 4, 30], without proving their applicability. A non-invertible tensor expansion which is mathematically true but physically weak must be applied during the derivation [4, 30]. The oscillating effect of the interpolated stress [15, 4] could also be merely a mathematical artifact rather than a reflection of the reality. Despite those problems, it offers many helpful insights to capture some essentials of the microscale stress field. First, it uses the local balances of linear momentum as the starting point, consistent with continuum mechanics, suggesting its potential compatibility with the stress determined by continuum mechanics. Regardless of the complicated derivation, the result is a summation of interactions between each interacting pairs, together with characteristic tensors in terms of distance vectors. This suggests a reverse approach: if we assume that along the line of interaction of an interacting pair the stress, when simply considered as force per unit area, over a certain characteristic area is dominated by their interatomic force, then we will be able to relate the stress directly to the atomic interactions by the Cauchy’s relation, skipping the complications brought by the differential equations. On the other hand, for the majority of atoms there are many interacting neighbors, which suggests that we can use the Cauchy’s relation to determine the resultant traction in many directions. Since the distance vectors can be obtained directly from the atomic positions, this then becomes an over-determined problem which can be solved by the least-square method. Of course, the results obtained would in general not be a very accurate one at every position, as the stress is varying from point to point. However, it must be kept in mind that we are not seeking for an accurate analytical expression of the microscale stress at any spatial position, but instead seeking for an approximation for a small region surrounding each individual atoms, i.e. the stress field in a small region around an individual atom is assumed to be constant. The entire stress field can then be interpolated from those point values by shape functions.

Comparing to the attention to the microscale stress, discussions on the microscale strain are rare in literature. The strain is usually considered as a given variable and used together with the microscale stress to obtain local elastic constants [15, 4]. However, in continuum mechanics, strain itself can be used as an indicator of plastic deformation [23] which, when translated into the microscale language, is the onset of dislocations. On the other hand, to estimate the microscale strain is equivalent to estimating the deformation gradient, since all kinds of strain tensors can be obtained from the deformation gradient [29]. In the early date, Born [2] related the deformation gradient to the atomic positions in the lattice by the so called Cauchy-Born rule, under the assumption that the entire lattice is under uniform deformation. This has been later used as the foundation for the well known multiscale method named as the Quasicontinuum method [26, 17]. If we imagine some virtual bonds exist between interacting pairs, then the assumption of homogeneous deformation means all the bonds are under the same deformation. For the microscale model where the particle methods must be applied, this assumption does not hold, i.e. the bond can be under different deformations. However, this does not mean that we have to abandon the Cauchy-Born rule completely. We can still apply it, but only to each bond. If, instead of seeking analytical expressions for the microscale strain, we want to obtain a reasonable engineering approximation for the microscale strain at individual atoms similar to what we hope for the stress tensor, then we can assume an averaged deformation gradient for each individual atom, and then apply the Cauchy-Born rule inversely, i.e. obtaining the deformation gradient from the undeformed and deformed distance vectors of each bond. As shown in the later section, as long as there are more than two linearly independent interacting bonds, such averaged deformation gradient can be estimated by the least-square method, just as the stress tensor. And the obtained deformation gradient can be further used to calculate the strain tensors at the positions of individual atoms.

In this short contribution, we will first derive the microscale strain tensor estimated by the least-square method, followed by the stress tensor. Then we show the results from the numerical verifications, where we compare the estimated microscale stress and strain tensors with the reference values. The numerical test is conducted on lattice model under simple loading conditions: uniform deformation and tensile test. Excellent matches between the least-square estimations and the reference values are obtained for the simple loading cases. After that, we present the results for more complicated loading conditions, including the stress field around a microscale crack, the stress field around an unit edge dislocation and the stress field under an indentation surface. In this contribution, we focus on the discussion of the method itself and the material parameters are of no interest, therefore the results of estimated stress field cannot be numerically compared to existing data. Nevertheless, the distribution pattern of the estimated stress shows good resemblance to the existing results in the literature. Limitations of our method are also discussed at the end.

For atomic systems with only pair potentials, the interaction between any two atoms is independent from the others and it can be treated as virtual bond between these two atoms. Local deformation gradient can then be assigned to each bond, as shown in Fig. 1 and the deformed distance vector $\mathbf{r}^{ij}$ from arbitrary atom $i$ to its neighbor, atom $j$, can be related to their undeformed distance vector $\mathbf{R}^{ij}$ by applying the Cauchy-Born rule [26, 2],

$$\mathbf{r}^{ij}=\mathbf{F}^{ij}\cdot\mathbf{R}^{ij}.$$

(1)

To avoid confusion in their directions, we define the distance vectors here explicitly as

	$\displaystyle\mathbf{r}^{ij}$	$\displaystyle=\mathbf{x}^{j}-\mathbf{x}^{i}$			(2)
	$\displaystyle\mathbf{R}^{ij}$	$\displaystyle=\mathbf{X}^{j}-\mathbf{X}^{i},$			(3)

where $\mathbf{x}^{i}$ is the deformed position of atom $i$ and $\mathbf{X}^{i}$ its original position. All interaction neighbors of atom $i$ are indexed by $j$. We consider the interpolated deformation gradient at atom $i$ as an averaged value, from where we obtain the approximated distance vector after deformation. Since $i$ is arbitrary, without loss of generality we can drop the index $i$ from the notations for simplicity.

$$\tilde{\mathbf{r}}^{j}=\mathbf{F}\cdot\mathbf{R}^{j}$$

(4)

One can rewrite the right hand side of equation (4) as a mapping of the components of the deformation gradient as

$$\underbrace{\left[\begin{array}[]{ccc}{\mathbf{R}^{j}}^{T}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&{\mathbf{R}^{j}}^{T}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&{\mathbf{R}^{j}}^{T}\\ \end{array}\right]}_{\displaystyle\mathbf{D}^{j}}\hat{\mathbf{F}}=\tilde{\mathbf{r}}^{j}.$$

(5)

where $\hat{\mathbf{F}}$ is the Voigt notation of $\mathbf{F}$, defined as

$$\hat{\mathbf{F}}=\left[F_{11},F_{12},F_{13},F_{21},F_{22},F_{23},F_{31},F_{32},F_{33}\right]^{T}.$$

(6)

The error of estimation of $\mathbf{r}^{j}$ is then

	$\displaystyle\Delta\mathbf{r}^{j}$	$\displaystyle=\tilde{\mathbf{r}}^{j}-\mathbf{r}^{j}$			(7)
		$\displaystyle=\mathbf{D}^{j}\hat{\mathbf{F}}-\mathbf{r}^{j}.$			(8)

The best estimation of $\mathbf{F}$ should minimize the error norm summed over $N_{p}$ interaction neighbors of atom $i$.

	$\displaystyle\Pi_{F}$	$\displaystyle=\sum_{j=1}^{N_{p}}\left(\Delta\mathbf{r}^{j}\right)^{T}\Delta\mathbf{r}^{j}$			(9)
		$\displaystyle=\sum_{j=1}^{N_{p}}\left(\mathbf{D}^{j}\hat{\mathbf{F}}-\mathbf{r}^{j}\right)^{T}\left(\mathbf{D}^{j}\hat{\mathbf{F}}-\mathbf{r}^{j}\right)$			(10)

Here we apply the standard linear least-square method, i.e. taking derivatives of equation (10) and equaling it to zero. It yields a linear system,

$$\sum_{j=1}^{N_{p}}\left({\mathbf{D}^{j}}^{T}\mathbf{D}^{j}\right)\hat{\mathbf{F}}-\sum_{j=1}^{N_{p}}\left({\mathbf{D}^{j}}^{T}\mathbf{r}^{j}\right)=\mathbf{0}.$$

(11)

By performing some matrix manipulation, the best estimation of $\hat{\mathbf{F}}$ can be obtained by

$$\hat{\mathbf{F}}={\underbrace{\left(\mathbf{D}^{T}\mathbf{D}\right)}_{\displaystyle\mathbf{M}_{F}}}^{-1}\mathbf{D}^{T}\mathbf{b}_{F},$$

(12)

where

	$\displaystyle\mathbf{D}$	$\displaystyle=\left[\mathbf{D}_{1}^{T},\mathbf{D}_{2}^{T},\ldots,\mathbf{D}_{N_{p}}^{T}\right]^{T}$			(13)
	$\displaystyle\mathbf{b}_{F}$	$\displaystyle=\left[\mathbf{r}_{1}^{T},\mathbf{r}_{2}^{T},\ldots,\mathbf{r}_{N_{p}}^{T}\right]^{T}.$			(14)

Once $\mathbf{F}$ is recovered, the strain tensor can then be easily computed from

$$\mathbf{E}=\frac{1}{2}\left(\mathbf{F}^{T}\mathbf{F}-\mathbf{I}\right),$$

(15)

where $\mathbf{I}$ is the 3-by-3 identity matrix. Notice that the matrix $\mathbf{M}_{F}$ in equation (12) depends only on the original lattice structure and has to be constructed only once, as far as no moving dislocation appears. On the other hand, to ensure that it is invertible, atom $i$ must have at least three interacting neighbors, whose distance vectors are linearly independent.

The stress estimation is a combination of the linear least-square method and the Cauchy’s relation [12],

$$\mathbf{t}=\boldsymbol{\sigma}\cdot\mathbf{n},$$

(16)

where $\mathbf{t}$ is the traction vector associated to the Cauchy stress tensor in a cutting plane with unit normal vector $\mathbf{n}$. In static case, the atoms are always in equilibrium. For each atom $i$, we have

$$\sum_{j=1}^{N_{p}}\mathbf{f}^{j}=\mathbf{0},$$

(17)

where $\mathbf{f}^{j}$ is the force obtained from the interatomic potential between atom $i$ and $j$. Therefore, if we consider a special volume surrounding atom $i$ with its surface cutting through the bond $\mathbf{r}^{j}$, then the total external force at the intersection point can be equivalent to the interatomic force between atoms $i$ and $j$, as shown in Fig. 2, where the traction near by can be approximated by

$$\mathbf{t}^{j}=\frac{\mathbf{f}^{j}}{A^{j}(\mathbf{x})}.$$

(18)

We define here $A(\mathbf{x})$ as some characteristic area, evaluated at position $\mathbf{x}$, for which equation (18) holds. Together with equation (17), one obtains

$$\boldsymbol{\sigma}(\mathbf{x})\cdot\frac{\mathbf{r}^{j}}{\left|\mathbf{r}^{j}\right|}=\frac{\mathbf{f}^{j}}{A^{j}(\mathbf{x})},$$

(19)

where $\left|\mathbf{r}^{j}\right|$ is the length of $\mathbf{r}^{j}$. For each atom $i$, we only estimate an averaged stress tensor within its surrounding space and therefore the unknown area $A^{j}(\mathbf{x})$ can be replaced by a constant value $A_{c}^{j}$ for each bond $\mathbf{r}^{j}$. The estimated stress tensor results in an approximated traction $\tilde{\mathbf{t}}^{j}$ for each bond with an error of

	$\displaystyle\Delta\mathbf{t}^{j}$	$\displaystyle=\tilde{\mathbf{t}^{j}}-\mathbf{t}^{j}$			(20)
		$\displaystyle=\mathbf{d}^{j}\hat{\boldsymbol{\sigma}}-\frac{\mathbf{f}^{j}}{A_{c}^{j}},$			(21)

where

	$\displaystyle\hat{\boldsymbol{\sigma}}$	$\displaystyle=\left[\sigma_{11},\sigma_{12},\sigma_{13},\sigma_{21},\sigma_{22},\sigma_{23},\sigma_{31},\sigma_{32},\sigma_{33}\right]^{T}$			(22)
	$\displaystyle\mathbf{d}^{j}$	$\displaystyle=\left[\begin{array}[]{ccc}{\mathbf{r}^{j}}^{T}&\mathbf{0}&\mathbf{0}\\ \mathbf{0}&{\mathbf{r}^{j}}^{T}&\mathbf{0}\\ \mathbf{0}&\mathbf{0}&{\mathbf{r}^{j}}^{T}\end{array}\right].$			(26)

Then we define the total error norm for the traction as

	$\displaystyle\Pi^{0}_{\sigma}$	$\displaystyle=\sum_{j=1}^{N_{p}}\left(\Delta\mathbf{t}^{j}\right)^{T}\Delta\mathbf{t}^{j}$			(27)
		$\displaystyle=\sum_{j=1}^{N_{p}}\left(\mathbf{d}^{j}\hat{\boldsymbol{\sigma}}-\frac{\left|\mathbf{r}^{j}\right|}{A_{c}^{j}}\mathbf{f}^{j}\right)^{T}\left(\mathbf{d}^{j}\hat{\boldsymbol{\sigma}}-\frac{\left|\mathbf{r}^{j}\right|}{A_{c}^{j}}\mathbf{f}^{j}\right),$			(28)

whose compact form can be written as

$$\Pi^{0}_{\sigma}=\left(\mathbf{d}\hat{\boldsymbol{\sigma}}-\mathbf{b}_{\sigma}\right)^{T}\left(\mathbf{d}\hat{\boldsymbol{\sigma}}-\mathbf{b}_{\sigma}\right),$$

(29)

with

	$\displaystyle\mathbf{d}$	$\displaystyle=\left[\mathbf{d}_{1}^{T},\mathbf{d}_{2}^{T},\ldots,\mathbf{d}_{N_{p}}^{T}\right]^{T}$			(30)
	$\displaystyle\mathbf{b}_{\sigma}$	$\displaystyle=\left[\frac{\left|\mathbf{r}^{1}\right|}{A_{c}^{1}}{\mathbf{f}^{1}}^{T},\frac{\left|\mathbf{r}^{2}\right|}{A_{c}^{2}}{\mathbf{f}^{2}}^{T},\ldots,\frac{\left|\mathbf{r}^{N_{p}}\right|}{A_{c}^{N_{p}}}{\mathbf{f}^{N_{p}}}^{T}\right]^{T}.$			(31)

The interatomic potential for pair-interactions can written be in terms of the distance between the two atoms, i.e. $\Phi^{j}=\Phi^{j}(\left|\mathbf{r}^{j}\right|)$, and therefore the force can be derived as

$$\mathbf{f}^{j}=\frac{{\Phi^{j}}^{\prime}}{\left|\mathbf{r}^{j}\right|}\mathbf{r}^{j}.$$

(32)

By substituting equation (32) into equation (31), $\mathbf{b}_{\sigma}$ can be simplified as

$$\mathbf{b}_{\sigma}=\left[\frac{{\Phi^{1}}^{\prime}}{A_{c}^{1}}{\mathbf{r}^{1}}^{T},\frac{{\Phi^{2}}^{\prime}}{A_{c}^{2}}{\mathbf{r}^{2}}^{T},\ldots,\frac{{\Phi^{N_{p}}}^{\prime}}{A_{c}^{N_{p}}}{\mathbf{r}^{N_{p}}}^{T}\right]^{T}.$$

(33)

Meanwhile, the conservation of angular momentum for non-polar continua requires $\boldsymbol{\sigma}$ to be symmetric, i. e.

$$\boldsymbol{\sigma}={\boldsymbol{\sigma}}^{T}.$$

(34)

When written in terms of elements it yields

$$\sigma_{12}=\sigma_{21},\sigma_{13}=\sigma_{31},\sigma_{23}=\sigma_{32},$$

(35)

which can be converted to matrix form,

$$\underbrace{\left[\begin{array}[]{ccccccccc}0&1&0&-1&0&0&0&0&0\\ 0&0&1&0&0&0&-1&0&0\\ 0&0&0&0&0&1&0&-1&0\\ \end{array}\right]}_{\displaystyle\mathbf{A}}\hat{\boldsymbol{\sigma}}=\mathbf{0}$$

(36)

The symmetric constraint defined by equation (36) can be imposed on the estimation by adding a penalty term to equation (29), as

$$\Pi_{\sigma}=\left(\mathbf{d}\hat{\boldsymbol{\sigma}}-\mathbf{b}_{\sigma}\right)^{T}\left(\mathbf{d}\hat{\boldsymbol{\sigma}}-\mathbf{b}_{\sigma}\right)+\lambda\left(\mathbf{A}\hat{\boldsymbol{\sigma}}\right)^{T}\left(\mathbf{A}\hat{\boldsymbol{\sigma}}\right),$$

(37)

where $\lambda$ is the penalty factor. By minimizing equation (37) one can obtain the best estimation of the stress tensor $\hat{\boldsymbol{\sigma}}$ for atom $i$.

$$\hat{\boldsymbol{\sigma}}={\underbrace{\left({\mathbf{d}}^{T}\mathbf{d}+\lambda\mathbf{A}^{T}\mathbf{A}\right)}_{\displaystyle\mathbf{M}_{\sigma}}}^{-1}\mathbf{d}^{T}\mathbf{b}_{\sigma}.$$

(38)

The matrix form of $\hat{\boldsymbol{\sigma}}$ can be recovered by equation (26). To guarantee $\mathbf{M}_{\sigma}$ is invertible, at minimum 3 neighbors with linearly independent distance vectors are required for each atom, the same as the condition for the strain estimation.

So far we have not been able to derive an analytical form for the characteristic area $A_{c}^{j}$ used in equation (29). By computational studies we found that the formula

$$A_{c}^{j}=\frac{4\pi\left(\frac{1}{2}\left|\mathbf{R}^{j}\right|\right)^{2}}{N_{p}}$$

(39)

provides the estimated $\mathbf{\sigma}$ with averaged values that are closest to the macroscale values. $\left|\mathbf{R}^{j}\right|$ is length of the undeformed distance vector between atom $i$ and $j$. It can be interpolated as a fraction of area of a spherical surface which is centered at atom $i$ and crosses the middle of the undeformed bond $\mathbf{R}^{j}$. And the fraction is determined by the number of interaction neighbors of atom $i$.

An atomic lattice model is constructed for conducting the numerical verifications. The interpolated microscale strain and stress fields are then compared with reference values.

The lattice is a face-centered Bravais lattice with a lattice constant of 4, as illustrated in Fig. 3. Lennard-Jones potential is applied, as

$$\Phi^{ij}=\frac{B}{\left|\mathbf{r}^{ij}\right|^{12}}-\frac{A}{\left|\mathbf{r}^{ij}\right|^{6}},$$

(40)

where $A$ is chosen to be $1.0239\times 10^{8}$ and $B$ to be $2.6211\times 10^{10}$. The equilibrium distance is $2.8284$ and the cutoff radius is $3$ so that only immediate neighbors are interacting with each other. Since here the material properties are of no interest in the test, all parameters are therefore dimensionless.

In our numerical tests, uniform deformation is applied to the boundary atoms by

$$\mathbf{u}_{bc}=\underbrace{\left(\mathbf{F}-\mathbf{I}\right)}_{\displaystyle\mathbf{H}}\mathbf{X}_{bc}.$$

(41)

When the static equilibrium is reached the entire lattice will be under the same deformation gradient [2], and the reference strain field can then be obtained by equation (15). On the other hand, the reference stress field is obtained by

$$\boldsymbol{\sigma}_{ij}^{ref}=\frac{1}{A_{i}}\sum_{k}\mathbf{f}_{j}^{i,k},$$

(42)

where $A_{i}$ is the area of the boundary surface $i$ defined in Fig. 3, and $\mathbf{f}_{j}^{i,k}$ is the $jth$ component of constraint force applied on atom $k$ on the boundary surface $i$. The surface change due to the deformation is ignored here, due to the small deformation applied.

On the other hand, in order to investigate the possibility of applying our interpolation to multiscale method, we also compared our stress interpolation with the Quasicontinuum stress (QC) field [26]. To obtain the QC stress field, a simple reference lattice can be constructed as shown in Fig. 4 whose radius is the same as the cutoff radius. The atoms at the center will be considered as the representative atom. Under uniform deformation, the 1st Piola-Kirchhoff stress tensor can be derived as [26]

$$\mathbf{P}=\frac{1}{2V}\sum_{j=1}^{N_{p}}{\Phi^{j}}^{\prime}\frac{1}{\left|\mathbf{r}^{j}\right|}\mathbf{r}^{j}\otimes\mathbf{R}^{j},$$

(43)

where $j$ is the index of the interaction neighbors of the representative atom and $V$ is the volume of the sphere surrounding the representative atom, at half of the cutoff radius. And the Cauchy stress tensor [29] can then be obtained by

$$\boldsymbol{\sigma}=J^{-1}\mathbf{P}\mathbf{F},$$

(44)

where $J$ is the determinant of the deformation gradient $\mathbf{F}$. The QC stress field represents the stress field for a crystal lattice of infinite size under uniform deformation.

For the first test, combined uniform strain field is applied to the lattice,

$$\mathbf{E}=\left[\begin{array}[]{ccc}0.010&0.010&-0.008\\ 0.010&0.010&0.010\\ -0.008&0.010&0.020\\ \end{array}\right].$$

(45)

The interpolation results for selected lattice planes are plotted in Fig. 5a and Fig. 5b, for the strain field and stress field accordingly, and the corresponding statistics are listed in Table. 1. Surface atoms are excluded from the final results. The recovered microscale stress and strain fields are uniform, same as the reference values. The strain estimation shows almost perfect match with the reference values. The stress estimation on the other hand shows about $20\%$ error compared to the reference values, as listed in Table. 1. This shift is caused by the surface atoms which are under unsymmetric interactions. We figured out that by combining the stress at each atom $i$ with its interacting neighbors $j$ by weighted summation, as

$$\bar{\mathbf{\sigma}}^{i}=\frac{1}{2}\mathbf{\sigma}^{i}+\frac{1}{2N_{p}}\sum_{j=1}^{N_{p}}\mathbf{\sigma}^{j},$$

(46)

the estimation for the internal stress can be significantly improved. The recovered stress values after applying the averaging technique equation (46) are given in Table. 2. One can observe that the internal stress field remains to be uniform while its accuracy is significantly improved, with the error reduced to the level of $2\%$. Moreover, on can observe that the stress field becomes much more closer to the QC stress field, which implies a promising application of the technique to the multiscale methods based on the QC framework.

Moreover, to investigate the performance of our interpolation changing with the severity of the deformation a numerical tensile test is conducted on the lattice model. An uniform strain in the z-direction is applied,

$$\mathbf{E}=\left[\begin{array}[]{ccc}0&0&0\\ 0&0&0\\ 0&0&E_{33}\\ \end{array}\right]$$

(47)

where $E_{33}$ varies from $-0.1$ to $0.1$ for $20$ loading steps. The recovered values of the normal stress in the loading direction are plotted in Fig. 6, together with reference values. The averaging technique defined by equation (46) is applied during the recovery. From Fig. 6 one can observe that the recovered stress values matches very well with the reference value and the QC stress field during the entire test. It can be also noticed that the QC stress field shows better consistency with the reference values than the previous test where a complicated strain field was applied. The unsymmetric response between tension and compression comes from the strong repulsion of the interatomic potentials when atoms are pressed together.

After the verification for simple loading cases, we then applied the interpolation method to study some more practical problems. In this contribution we show the interpolated stress fields around a micro crack, around an unit edge dislocation in the FCC lattice and below an indentation surface. Figure. 7a shows the Von Mises stress around a micro crack, generated by removing a single layer of atoms in the middle of lattice. The lattice is modeled in the shape of a thin plate with a thickness of 10, and is uniformly pulled in the $[001]$ direction. The figure shows that the highest stresses interpolated by our method occur at the corners of the crack and stretches in the diagonal directions, whereas the atoms right above and below the crack surface experience almost no stress. The places with high stress concentrations can be considered as potential positions for dislocations.

In Figure. 7b, an unit edge dislocation with a burger vector of $\frac{1}{2}[110]$ is created in the middle and the resultant stress field $\sigma_{zz}$ after relaxation is illustrated. Due to the limitation of computational power, only immediate neighbors are interacting with each other. The interpolated result indicates a symmetric distribution below and above the core, resembling the pattern predicted by the continuum theories [13]. However, this symmetry is distorted by the boundary effects. On the other hand, the dislocation introduces a pair of symmetric concentrations far away from the dislocation, indicating a tendency of splitting during the relaxation. This could be an intrinsic behavior of the dislocation or could merely be an influence of the boundary effects.

At the end, examples of stress distribution under an indentation surface are shown in Figure. 7c and Figure. 7d, where a square indentor is pressed into the center of the surface. The figure shows the distribution of the Von Mises stress. In Figure. 7c one can observe that when there is no defect, stress is concentrated at the corners and the bottom of the indentor. When there are defects (point vacancies) closely beneath the indentation surface, the stress distribution is notably changed. Stress is severely concentrated around the defects and the maximum value is significantly increased. On the other hand, under the combined effects of the indentation and the lattice defects, the stress in certain region is observed to be lowered.

In this contribution we introduced a straightforward approach to interpolate the microscale strain and stress for particle systems. The numerical test shows very good consistency between the recovered strain field and the applied strain field for simple loading conditions. The stress estimation shows relatively large errors compared to the reference values, which can be effectively reduced by applying averaging techniques. On the other hand, we observe good consistency between the interpolated microscale stress field and the stress field obtained by the Quasicontinuum method, which indicates potential applications of the method in multiscale methods. The tensile test shows that the error between the interpolated stress field and the reference stress field stays at low level during the entire loading process, further justifying the applicability of our method. For complicated loading conditions, such as dislocations and indentations, we observed that the patterns of the estimated stress resembles the results from the literature. There are also some foreseen limitations associated with our method. First, there is no dynamic effects considered for the stress estimation, or in other words the influence of mass flow, if any, to the stress is not considered. Second, it is possible that the results given in this contribution only valid for particles systems under pair potentials, nevertheless, we believe similar procedures can be applied to systems with many-body interactions.

The authors gratefully acknowledge support from the Graduate Research Group 615 (GRK615), funded by the German Research Foundation (DFG).