Efficient and Reliable Nanoindentation Simulation by Dislocation Loop Erasing Method

Nanoindentation is a useful technique to measure material properties at microscopic level. However, the intrinsically multiscale nature makes it challenging for large-scale simulations to be carried out. It is shown that in molecular statics simulations of nanoindentation, the separated dislocation loops (SDLs) are trapped in simulation box which detrimentally affects the plastic behavior in the plastic zone (PZ); and the long-distance propagation of SDLs consumes much computational cost yet with little contribution to the variation of tip force. To tackle the problem, the dislocation loop erasing (DLE) method is proposed in the work to alleviate the influence of artificial boundary conditions on the SDL–PZ interaction and improve simulation efficiency. Simulation results indicate that the force–depth curves obtained from simulations with and without DLE are consistent with each other, while the method with DLE yields more reasonable results of microstructural evolution and shows better efficiency. The new method provides an alternative approach for large-scale molecular simulation of nanoindentation with reliable results and higher efficiency and also sheds lights on improving existing multiscale methods.


Introduction
Nanoindentation is an effective technique not only to characterize mechanical properties of materials like hardness and elastic response at high spatial resolution in experiments [1], but also to gain insights into plastic mechanisms at atomic level by molecular simulations [2]. In essence, nanoindentation is a typical multiscale problem involving mechanical phenomena at several scales. Firstly, the instability of individual atoms under the indenter tip leads to dislocation nucleation [3]. Afterward, the growing dislocation network constitutes the plastic zone (PZ) which accounts for the hardness, while some dislocations are emitted into the deep of substrate or move along surfaces as separated loops [4]. The expanding PZ reveals the material plasticity, and size effects emerge due to the evolution of dislocation density [5]. Nanoindentation can be directly simulated with empirical-based atomistic methods, such as molecular dynamics (MD) or molecular statics (MS). The computationally intensive nature of nanoindentation, however, limits the simulation scale to a size much smaller than that in experiments. Therefore, developing new methods to simulate nanoindentation with higher reliability and efficiency is of great importance for the research community. Over the past two decades, various simulation methods have been put forward to explore the microscopic mechanism from atomic motion to dislocation evolution, providing in-depth insights into the understanding of experimental results and mechanisms underpinning the plasticity in atomic level.
To make simulations of nanoindentation more realistic, two of the most commonly used methods are atomistic simulation using periodic boundary conditions (PBC) [6] and multiscale modeling [7]. In atomistic simulations, the lateral boundaries of simulation box are usually set as periodic, and the bottom layers with thickness of several nanometers are fixed to support the substrate against translation [6]. PBC can effectively cancel surface effects, but the periodic and fixed boundary conditions will impede the propagation of separated dislocation loops (SDLs) and confine all dislocations in a small volume, and as a result, the residual SDLs will affect the dislocation evolution in PZ. One possible solution is to restrict the maximum depth to 1/10 of the specimen height, but more atoms are needed for deeper indentation, leading to the dramatic increase in computational burden. Alternatively, PBC combined with partially free bottom (PFB) can be used to allow prismatic dislocation loops (PDLs) to escape from the simulation box [2,8,9], such that PDLs have less influence on the evolution of PZ in indentation and retraction stages. However, some half prismatic dislocation loops (HPDLs) can move along the upper surface and enter the simulation box again from the opposite side due to PBC. It can be seen that neither the fixed nor the free bottom boundary can handle this issue. From the viewpoint of efficiency, our previous work shows that long-distance propagation of SDLs will lead to the nonlinear increase in computational cost [10]. Therefore, atomistic simulation using PBC suffers from both accuracy and efficiency issues due to movements of SDLs.
On the other hand, multiscale simulations attract increasing research interest over the past decades based on the idea of reducing degrees of freedom by coupling multiscale representations concurrently or hierarchically [11][12][13][14][15]. Specifically, dislocations and relevant details are captured by atomistic representation, while regions under homogeneous deformation are described by quasi-continuum representation in order to reduce computational cost. The quasi-continuum method (QC), one of the most typical multiscale methods, uses representative atoms to find local stable states by minimizing the coarse-grain potential energy at temperature 0 K [16]. Theoretically, the QC method is expected to simulate nanoindentation of large-scale systems and the efficiency problem could be solved. However, SDLs are emitted frequently from the PZ and move into the substrate or along the surface during indentation, leading to the rapid expansion of atomistic representation regions. Consequently, the efficiency improvement in QC over fully atomistic simulation is effective only when indentation depth is relatively small [17]. Furthermore, the coupled atomistic and discrete dislocation plasticity (CADD) method extends QC by introducing discrete dislocation plasticity. It is shown that CADD has the potential for multiscale modeling of crack initiation and propagation as well as nanoindentation [18]. Recently, the theoretical framework of three-dimensional (3D) version of CADD has been completed [19], which, however, has not been applied to large-scale nanoindentation [20]. Besides, whether the dislocation transition between different representation scales has influence on PZ remains unclear. Therefore, although multiscale methods might be the promising approach to solve the boundary effect and efficiency problems in nanoindentation, there is still a long way to go.
To clarify and solve the aforementioned issues, nanoindentation simulations using MS are carried out to investigate the influence of SDLs on computational accuracy and efficiency. Based on the results and analysis, the dislocation loops erasing (DLE) method is proposed and verified to improve the accuracy and efficiency by controlling the motion of SDLs. MS simulations with DLE are applied to the nanoindentation system of single crystal Cu, showing that reasonable dislocation evolution in deeper indentation can be obtained with higher efficiency. During the simulation, the tip moves toward the substrate with an incremental displacement of 0.1Å, and then the system is relaxed by the conjugate gradient (CG) method. The ratio of the maximum indentation depth h max (3.9 nm) to the thickness in the z-direction, L z (21.7 nm), is greater than 1/6. Small h max /L z (usually less than 1/10) is required to assure the accuracy, which, however, will increase the simulation burden. Here, we intentionally use a larger h max /L z to examine the influence of PBC on simulation results.

Simulation Configuration
The interactions between Cu atoms are described by the embedded atom model (EAM) developed by [21], which has been widely applied in atomistic simulations with reliable results [22]. The Morse potential is adopted to describe the Cu-C interaction [23]: where D 0 = 0.087 eV, α = 5.14Å, and r 0 = 2.05Å. The indenter is assumed to be rigid so that the C-C interactions can be ignored. The force evaluation times (n eval ) are used to quantify the computational cost during CG minimization for each loading step [24]: where n k line denotes the force evaluation times during the linear search process and n iter denotes the iteration steps required for convergence. The atomic energy difference is used as the criterion of convergence: where N is the atom number, e i k and e i k+1 are the atomic potential energies at the kth and (k + 1)th iteration steps, respectively. ε = 1 × 10 −5 eV is adopted as the convergence criterion for energy minimization. OVITO is used to visualize atomistic configuration [25]. The centrosymmetry parameter (CSP) [26] and dislocation extraction algorithm (DXA) [27] are used to identify defect atoms and dislocation lines, respectively.  Table 1. It can be seen that although the PZ occupies a small region of the simulation box, SDLs move to  the boundary which breaks the spatial localization of dislocations. Gao et al. carried out systematic molecular simulations of nanoindentation for FCC and BCC metals and found that the ratio of the radius of the PZ to the contact radius lies between 2 and 3, indicating the localization of PZ [4]. The localized PZ and movements of SDLs in our simulations are in agreement with those of previous research [4,28]. When SDLs are blocked by PBC and fixed boundary, dislocation behaviors in PZ and the variation of tip force will be affected. Figure 3a shows the evolution of tip force with respect to indentation depth. Since a conical indenter is used in the simulation, the initial elastic stage of the F -h curve is short (about 0.2 nm), compared with the result using the spherical indenter [9]. The curve increases stably before 1.5 nm, during which dislocations grow steadily without DL emission. After h exceeds 1.5 nm, the curve increases with a lot of drops as a result of inhomogeneous deformation. The whole indentation process can be divided into two stages according to the different responses of tip force to dislocation evolution: the uncrowded stage (h < 3 nm) and the crowded stage (h > 3 nm). In the uncrowded stage, as shown from point 1 to 5 in Fig. 3a, the sharp drops correspond to the emissions of the first five SDLs in Table 1. The other mini drops in the stage are related to the dislocation reaction and annihilation within the PZ. In the crowded stage, three DLs are emitted from the PZ, causing slightly drops of the F -h curve labeled by 6, 7, and 8 in Fig. 3a, whereas the rest of force drops result from SDL-PZ interactions. For example, the fourth PDL in [101] extends the PZ at the maximum depth, and the sixth HPDL and eighth HPDL move back to squeeze the PZ as shown in Fig. 3c. Therefore, formation and emission of DLs are the major mechanism that causes the drops of F -h curve in the uncrowded stage. The F -h curve is dramatically affected by the SDL-PZ interaction in the crowded stage due to the existence of PBC and fixed boundary. Figure 3d shows the dislocation pattern after full retraction. Upon tip retraction, the first PDL, sixth HPDL, and seventh PDL return to the PZ and interact with it. In particular, the first PDL is merged into the dislocation network and extends the size of PZ. Experimentally, however, most emitted DLs will not return to PZ but quickly move deeply into the substrate material [29]; thus, the phenomena in Fig. 3d would be different in a larger system. Figure 3b presents the computational cost of minimization processes for all loading steps. Similarly, the computational cost is significantly affected by dislocation behaviors. According to our previous work [10], long-distance propagation of SDLs is the main cause of the nonlinear increase in computational cost. It can be seen that much more computational cost is required for steps with SDL emissions than for normal steps, and the former could be an order of magnitude higher than the latter. Notably, emissions of the sixth and seventh DLs in the crowded stage are restrained by existing SDLs, and thus less computational cost is needed to move the DLs in a short distance. The abnormal peak labeled A in Fig. 3b has no connection with SDL emission, but is due to the further propagation of the 6th HPDL. When h = 3.5 nm, the sixth HPDL moves back, leading to the extra high computational cost at point 8 in Fig. 3b. The abnormal peaks in the crowded stage in Fig. 3b, therefore, also indicate disordered evolutions of dislocations in deep indentation.

Nanoindentation Simulation Using PBC
To explore more details of the artificial boundary effects and the abnormal computational cost caused by SDLs, Fig. 4 presents the complete evolution of tip force in the energy minimization process when h = 1.53 nm, during which a DPL and a HPDL are formed and propagate along respective slip directions. The whole minimization process can be divided into three stages according to the reduction magnitude of tip force as shown in Fig. 4a. In stage I, atoms near the tip adjust positions to reduce the high repulsive force caused by tip penetration. Meanwhile, an HPDL moves along the [110] direction and a new PDL is formed and moves along the [011] direction. Dislocations in the PZ keep evolving in this stage. It is worth noting that 99% of the reduction in tip force happens in this stage with only 10% of the total computational cost. As the iteration continues in stage II, as shown in Fig. 4b, the PDL and HPDL continue moving and the tip force decreases slowly, while the PZ keeps unchanged afterward. The PDL stops moving at the end of stage II, and the tip force decreases by only 0.6%. In stage III, dislocations hardly move and the tip force stays nearly constant. The HPDL reaches its equilibrium position when it arrives at the corner of the box. Therefore, the atomic rearrangement in PZ and dislocation evolution in stage I make the major contribution to the decrease in tip force, while 90% of the total computational cost is consumed in stages II and III.
To investigate the abnormal computational cost due to the long-distance propagation of SDLs, the work of Gerberich et al. [30] is used to estimate the theoretical equilibrium position of SDLs after emission. The local stress acting on an SDL is given by where p 0 is the maximum pressure, d i is the distance from the DL to the tip, a is the contact radius, and v is the Poisson's ratio (0.33 for Cu). The maximum pressure for a spherical contact is 3F/2πa 2 , where F is the tip force. Although the conical tip is used in our simulation, Eq. (4) can provide an instructive estimation for general nanoindentation. Here, we ignore the image force from the free surface and consider the propagation of a single SDL; thus, only external driving force caused by indentation takes effect. Figure 5 shows the local stress acting on a single SDL when it moves away from the PZ under different tip forces. The horizontal line represents the Peierls stress of Cu [31] which is the barrier that needs to be overcome for SDL movement. Therefore, the intersection points of the horizontal line with different curves define the equilibrium positions where driving forces are equal to the lattice friction. It can be found that, for F = 200 nN (the first DL emission), the theoretical moving distance of the SDL reaches 562.5 nm, and for the tip force F = 1000 nN (at the end of indentation in Fig. 4a, an SDL can move 1284 nm away from the tip. Therefore, the theoretical equilibrium position of DLS is up to micrometers from the tip due to the low Peierls stress of Cu, which goes far beyond the capability of atomistic simulations. As a result, SDLs will always reach the boundaries of a simulation box with sizes in hundreds of nanometers. From the viewpoint of accuracy, SDLs are formed at h = 1.53 nm and move a long distance of 22 nm and 28 nm, respectively, until they are stuck by PBC in the current simulation, leading to unphysically strengthening of the tip force. From the viewpoint of efficiency, one can see that the larger the simulation box, the more computational cost will be consumed in stages II and III. A trade-off between accuracy and efficiency needs to be addressed for the current nanoindentation simulation. In nanoindentation simulation, if the computational efficiency is preferential, a small configuration can be used with deeper indentation. In this case, SDLs will get stuck in a small volume and the corresponding high computational cost due to long-distance propagation of SDL can be avoided. However, the trapped SDLs will interact with PZ, leading to the loss of accuracy. On the other hand, a large configuration can mitigate the SDL-PZ interaction, but increase computational cost in two aspects: (i) the larger system needs more CPU resources to calculate interatomic forces and update positions for extra atoms, and (ii) more SDLs will nucleate and propagate with longer distances accompanied with nonlinear increase in computational cost [10]. Although some multiscale methods such as QC can handle the first aspect well by using representative atoms to reduce degrees of freedom, they still suffer from the second aspect. CADD seems to be hopeful to solve the issue, but its application in 3D simulations is still ongoing [32]. The efficiency problem originated from SDLs, therefore, is the common challenge for both atomistic and multiscale simulations in nanoindentation and becomes the motive for the development of a new simulation method.

Dislocation Loop Erasing Method
As discussed in Sect. 3, the destination of SDLs in a large system should be far away from the tip and PZ for the materials with low Peierls stress. Since we primarily concern the tip force and dislocation evolution in the PZ after minimization, the results will not change whether DLs move slowly or instantly to the destination. Therefore, DLs can be erased on appropriate occasions to keep dislocation events locally, such that the SDL-PZ interaction can be alleviated with higher accuracy and the long-distance propagation of SDL can be avoided with higher efficiency. Here, a special method called dislocation loop erasing (DLE) is proposed to achieve the goal. Figure 6a shows the schematic of the implementation of the DLE with PBC. DLs will be erased directly from the simulation box once they move to the distance of R DLE . The determination of R DLE will be discussed below. Lattice steps are supposed to form at the outer boundaries like PFB to maintain mass conservation. Considering that the atoms near the lattice steps will not affect the nanoindentation process, they can be erased directly to restore the original boundaries. Using the DLE operation, only small atomic region with radius of R DLE is needed. Moreover, there is less SDL-PZ interaction once SDLs are erased, and much computational cost can be saved without long-distance propagation of SDLs. Notably, a special boundary proposed by Li et al. [2], van Vliet et al. [8], and Lee et al. [9] in MD simulations allows PDLs to move out of the box from PFB with parallel lattice steps being left at the bottom of substrate. The original intention of this method is to avoid PDL-PZ interactions. Alhafez et al. compared the difference of PZ evolution using different boundary conditions at the bottom and found that PFB indeed can effectively alleviate the influence of PDL on PZ evolution [28]. Although PFB fails on HPDLs due to different slip planes, it is shown that erasing SDLs directly is an efficient and reliable way to transport mass and plasticity.
One concern about DLE is its influence on tip force. Since it is the lattice friction that impedes the movement of SDLs, if SDLs are erased, contribution of lattice friction on tip force will be cancelled. However, this contribution can be ignored due to the very low Peierls stress of Cu. Gilman listed several reasons to demonstrate that Peierls stress in pure metals is negligible [33]. In the current simulation,  Fig. 4b also supports this point for Cu. The HPDL travels about 10 nm in this stage, while the tip force keeps almost the same level, indicating that propagation of the HPDL has little influence on tip force, and the contribution of lattice friction can therefore be ignored.

Procedure of DLE
There are two tasks in DLE to handle SDLs reasonably and efficiently: (1) SDL detection To detect atoms associated with PZ and SDLs, CSP between 0.5 and 3 is used to identify "defect atoms." Cluster analysis is utilized to divide atoms into different groups according to interatomic distances. The atom group of PZ can be easily found beneath the tip, while SDLs can be distinguished from PZ according to their relative positions.
(2) SDL erasing Once an SDL is detected, a minimum box containing the SDL is marked, as shown in Fig. 7a. The outer layer of the box should have no defect atoms. Basis layers on the upper and lower surfaces of the covering box are shown in Fig. 7c-d. All atoms within the box except the two basis layers are removed. Finally, the removed region is reconstructed by replicating the basis layers periodically according to the stacking order of lattice structure. Figure 7b shows the atoms with CSP > 0.5 after reconstruction. The reconstructed region looks disordered since there is a little mismatch between the basis layers and the original layers. However, it is shown that such disorder will disappear after several minimization iterations, and after that, SDLs are erased completely. The same operation is applicable for either PDL or HPDL.

Determination of R DLE
DLE takes effect whenever an SDL moves to R DLE . To illustrate the influence of R DLE , Fig. 8a plots the evolution of distance between SDL and the tip during minimization when h = 1.53 nm. Similar to Fig. 4, three stages remark different motion behaviors of the PDL and HPDL. Five R DLE located in different stages are selected to examine the influence of R DLE on the response of tip force. As shown in Fig. 8b, the original curve represents the tip force evolution without DLE, and more than 40,000 force evaluations are required to satisfy a given convergence criterion. The other curves represent the results using DLE with different R DLE , in which the minimization process finishes earlier to some extent. It is also worth noting that the final tip forces with DLE are lower than that without DLE after convergence.
The difference between the tip forces with and without DLE (δF ) originates from the strengthening force caused by SDL-PZ interactions. When the DLE is invoked with R DLE in stage I, δF equals to 2.7 nN. Snapshots of atomistic configuration show that applying DLE in stage I will affect the dislocation evolution in the PZ. However, when DLE is invoked in stages II and III, δF is obtained as the same value of 1.4 nN. Since the dislocation evolution of PZ has almost completed after stage From the aspect of efficiency, R DLE should be as small as possible to reduce computational cost. For example, the computational costs with R DLE = 11.68 nm and 21.75 nm are, respectively, 16.8% and 55.4% of those without the DLE operation. Obviously, DLE simulations with R DLE of 11.68 nm can save more computational cost. Since the specific value of R DLE depends on the indented materials, indentation depth, and size of PZ, it is a challenge to predict R DLE before simulation. Physically, the analysis in Fig. 8 reveals that R DLE should be the least distance between the DL and PZ once the SDL enters stage II, which can be directly reflected from the evolution of tip force in Fig. 4a. Therefore, it is proposed that the reduction magnitude of tip force can be used as a simple indictor to determine the minimum R DLE as shown in Fig. 4a, which is a relative quantity and easy to monitor regardless of the indented materials, indentation depth, and size of PZ. The SDL detection in DLE is similar to the

Nanoindentation Simulation with DLE
To investigate the flexibility and efficiency of DLE, nanoindentation simulations with DLE and PBC are carried out. Figure 9 shows the flowchart of the new method, in which the DLE block is inserted into the standard MS procedure. Since the CSP calculations of all atoms for SDL detection consume much computational cost, it should be invoked conditionally. Considering the propagation speed of SDLs, detecting SDLs every N SDL ≥ 500 iterations is suggested. The optimal N SDL can be determined by trial simulation tests.
The same configuration in Fig. 2 is used in the new simulations with the maximum indentation depth of 3.90 nm. Figure 10a presents the nine erased SDLs with information of their formation depth and dislocation length calculated by DXA [34]. The δF associated with each HPDL and PDL is plotted in Fig. 10b. For both HPDL and PDL, δF increases with indentation depth because SDL-PZ interactions become stronger as SDL and PZ get closer. It is seen that the extra strengthening effect also depends on the SDL type. As shown in Fig. 10b, δF of PDL is much higher than that of HPDL, which can be attributed to two reasons: (i) the size of HPDLs is usually half of that of PDLs at the same depth, and thus PDLs exert larger force on the tip; (ii) the tip force is the sum of atomic forces in the vertical direction, while the motion of HPDLs is horizontal. Figure 11a shows the F -h curves with and without DLE. The depth of SDL emission is marked in the two curves. In general, the two curves are close to each other. Once DLE is invoked at h = 1.53 nm, the two curves start to diverge slightly due to the different numbers and appearance sequences of SDLs. Figure 11d shows the dislocation pattern at the end of indentation when all SDLs have been erased. In comparison with the pattern from the original simulation in Fig. 3c, there are no residual SDLs in the space between PZ and the boundaries. Moreover, no SDL returns to PZ after full retraction, as shown in Fig. 3d. Therefore, the unphysical strengthening effect due to boundary effects has been well relieved after DLE.
Another significance of DLE in nanoindentation is efficiency improvement. In Fig. 11b, it is seen that the efficiency is improved whenever DLE occurs, such that the nonlinear increase in computational cost found in the original simulations [10] can be mitigated. Figure 11c shows the comparison of computational cost with and without DLE when SDLs are emitted. For the first DLE, the HPDL should have moved from the vicinity of tip to the corner of the box in Fig. 4, but DLE saves the long-distance propagation. The computational cost with DLE is only 16.8% of that in the original simulation. The remaining DLE operations have non-obvious efficiency improvement because the SDLs are trapped in the small simulation box and no long-distance propagation is required. The efficiency improvement by DLE, however, will be more significant for larger-scale simulations in which more DLs will propagate in long distance. On the other hand, the new method is more reasonable considering the dislocation pattern in PZ and the strengthening effect due to boundary conditions. To achieve the comparable accuracy, a much larger configuration has to be used in the traditional MS method. In this sense, the implicit efficiency improvement in the new method will be more attractive.

Applicability of DLE to Other Materials
It is worth noting that the proposed method is applicable when some dislocations exist in local region, while some others move far away in the form of SDL. Fortunately, such scenario is quite common for metals such as FCC Al, Cu, BCC Fe, Ta [4] and HCP Mg, Ti, Zr [28], and for some zincblende structure ceramics [35,36]. Although DLE removes the long-range elastic interaction between SDLs and the rest of atoms, it has less influence on simulation accuracy for materials with low Peierls stress. Recent studies have shown that the Peierls stress is very small for most FCC and HCP metals (< 1 MPa), moderate for BCC metals (390-960 MPa), and relatively high for zinc-blende materials (1300-4000 MPa) [31]. The motion of SDLs can be effectively constrained in materials with high Peierls stress, such that DLE is unnecessary for zinc-blende materials. Nonetheless, considering the magnitude of hardness in MD simulations (∼GPa), DLE provides a more efficient and reliable simulation approach for large-scale nanoindentation of FCC and HCP metals.

Conclusion
The accuracy and efficiency issues caused by SDLs in nanoindentation simulations are investigated systematically and a DLE method is proposed in this paper to tackle the problem. In nanoindentation simulations with PBC, it is found that SDLs are strongly related to the breakage of dislocation localization, strengthening of tip force, and deterioration of computational efficiency. The long-range propagation of SDLs consumes much of the computational cost, but has little contribution to the tip force. The DLE algorithm is designed to release the restriction of SDL motion and alleviate the SDL-PZ interaction due to the existence of artificial boundary conditions. Simulation results indicate that the force-depth curves obtained from simulations with and without DLE are consistent with each other, and the simulation with DLE shows more reasonable results of microstructural evolution and better efficiency.
Moreover, the idea of DLE can be applied to the development of multiscale methods. Since DLE gets rid of the long-range propagation of SDLs, only small regions around the PZ need fully atomistic description and simulation with deeper indentation could be achieved. The analysis also shows that it is unnecessary to model SDLs with atomistic representation when it is far away from the tip, because they have little influence on the dislocation evolution of PZ or the variation of tip force.