Edge dislocations in Ni monocrystalline structure studied by McChasy 2.0 Monte Carlo code

The main goal of the McChasy code, in general, is to reproduce Rutherford Backscattering Spectrometry experimental spectra recorded in channeling direction (RBS/C) by simulating He-ions travelling inside crystalline structures and to calculate the probability of the backscattering process. The 2.0 version of the code provides the possibility to simulate channeling spectra in large (ca. 108 atoms) arbitrary structures that are created based on crystallographic data or Molecular Dynamic (MD) calculations. In this work, we present the current status of the code and the results of recent investigations of extended structural defects (edge dislocations and dislocation loops) formed inside nickel (Ni) single crystals. Two ways of modelling extended defects are described: one developed using the McChasy code (Peierls-Nabarro approach) and the other one obtained by modification and thermalization of Ni structures by MD (LAMMPS code). The atomic local environment was studied qualitatively and quantitatively by the local projectile-flux density distributions around the defects.


Introduction
Radiation defects in different materials have been extensively studied by many groups over the past several decades. The Rutherford Backscattering Spectrometry in channeling direction (RBS/C) technique has been used as a standard method for the analysis of structural properties for ion-implanted single crystals by many authors [1][2][3][4]. Unfortunately, a lack of the proper tool of RBS/C spectra analysis and oversimplified methodology that has been usually applied could cause misleading results. Hence, it is highly important to develop an appropriate tool that can perform detailed, quantitative analysis separately for various types of defects that are formed in investigated crystals.
McChasy v.1.0 was developed in the late eighties in the National Centre for Nuclear Research [5,6]. The main principle of the first version of the code was to reproduce RBS/C experimental spectra by simulating the movement of He-ions travelling inside 2 a crystalline material. The structure is built on the single unit cell of the considered crystal surrounded by a few/few tens of neighbouring atoms. The code also calculates the probability of the backscattering process. Defects are considered by displacing some atoms "in situ" (during the simulations) based on appropriate models and distributions of defects defined by the user. The new version of McChasy code (v.2.0) is a response to recent interest in materials science and the huge potential of atomistic and molecular simulations, having the possibility to simulate channeling spectra in large (hundreds of millions of atoms) arbitrary structures. Note that the new McChasy 2.0 code can work with different file formats and is compatible with almost every program that provides positions of atoms in a structure of interest as an output. A schematic draft of the main individual procedures of the McChasy 2.0 code is shown in Figure 1. The code uses four repositories of files. TABLE OF TARGETS contains crystallographic data (space groups, lattice parameters, Wyckoff symbols, and initial positions of atoms) for all structures already available for simulations. The table can be upgraded with new structures, if necessary. Using a procedure piece of crystal, the code can save a file containing positions of atoms of a selected crystal in a format readable by MD (MD_INPUT) or read a structure file derived by MD (MD_OUTPUT). No matter how it is made, the structure file is reconstructed using a procedure set of atomic slices to extract atomic slices orthogonal to the direction chosen as the channeling direction. Such prepared samples are saved as a text file (SAMPLES). The simulation procedure is based on calculations of the nuclear encounter probability (NEP) [7] that depends on the impact parameter and energy of He-ions. The procedure is described in detail in Ref. [8]. The possible output of the code is denoted as SIMULATION RESULTS in Fig. 1. An additional feature of the code is the possibility to extract some structural data from the structure files, e.g., positions of atoms along a selected atomic row or within a cross-section at a selected depth (denoted as STRUCTURAL RESULTS in Fig. 1).   [9]. The huge potential of the early stage of the code was already presented elsewhere [8,10]. Already reported simulations performed for magnesium aluminate spinel (MgAl2O4) result in a very good agreement with experimental spectra and also with the previous version of the code (McChasy 1.0). The spinel structure was chosen previously in a manner to simulate only point defects that are formed in the material. The agreement with the previous version of the code proves that the point defects modelled as randomly displaced atoms (RDA) are analyzed similarly in both programs and that the McChasy 2.0 fits the experimental data correctly. The fact that should be highlighted is that the models of defects in version 2.0 of the code are easier to implement and the implementation itself in the large structures is giving much more possibilities. The code can be combined and work with almost every program giving atomic positions (e.g. with analytically added defects). In this paper, we are focused on extended defects (i.e. edge dislocations) and we combine the code with MD-LAMMPS code. The merging gives a possibility to thermalize analytically developed defects (and the structure itself). The use of peripheral add-ons allows also the extraction of structural properties or the reproduction of the spectra used for several structural techniques like TEM or XRD [11,12]. Because of that, we can eventually obtain a full structural picture of a modelled sample. The current status of an early testing phase of the code and recent simulation results are discussed and presented in the paper using the example of Ni monocrystalline structure. Nibased alloys are nowadays one of the most promising materials that can be used in the nuclear industry and in general in the power generation sector due to their radiation resistivity and immune response to harsh environmental conditions [13][14][15].

Methods
Having already experience regarding the RBS/C simulations of single-phase concentrated solid solution alloys (SP-CSAs), the current investigations were done using mainly the representation of extended, dislocation-like defects. The model of them was built based on the Peierls -Nabarro dislocation approach [16]. In both MC programs (McChasy 1.0 and 2.0), the edges of dislocations, as well as their Burgers vectors, are oriented perpendicular to the analyzing ion beam. RBS/C is highly sensitive to such orientated dislocations, which contribute to intensified dechanneling of the beam resulting in increased backscattering yield measured in the experiment [17]. While in version 1.0 dislocation edges are tossed randomly in four (for cubic structures) or six (for hexagonal ones) crystallographic directions [18,19], in the 2.0 version the user may choose between different edge orientations. For the study, we used a/2<100> oriented dislocations. The other parameters and simulation conditions were chosen carefully to match the experimental ones. A detailed description of the simulation procedures can be found elsewhere [6,18]. A simulated structure (called a sample) is about 30 nm x 30 nm wide and 200 nm long. The cubature, of course, can be adjusted to meet the appropriate conditions. Such samples typically contain from millions up to a few hundred million atoms. The MD simulations provide models of defected Ni crystal structures at T=300 K and p=10 atm. For each case well-thermalized atom positions were obtained after 20,000 steps, which correspond to 20 ps of real-time, the step interval was set equal to 1 fs; as a gauge of a final thermalization, the fluctuation of total energy less than 1% was assumed. The temperature  [20], which emulates the NPT ensemble, where N corresponds to the constant number of atoms, P -constant pressure and Tconstant temperature. To avoid spurious surface effects, calculations were done under periodic boundary conditions, which create virtual infinite crystals. The final 250 steps of the simulation were recorded and next averaged to produce atomic positions free of thermal motion. Interactions between atoms for MD simulations are approximated by potential energy functions. In this case, the Embedded Atom Model (EAM) potential with parameters given by Foiles was used [21]. For graphical interpretation and 3D visualization of atomistic data obtained from molecular dynamics or Monte Carlo simulations, the custom applications combined with Origin or Gnuplot and Open Visualization Tool (OVITO) software were used [22].

Discussion
In the McChasy code, the interaction between an ion and an oriented (aligned) monocrystalline structure is considered at every atomic layer by calculating the probability of the projectileatom collision with the impact parameter small enough that the projectile is backscattered [7]. Due to the atomic thermal vibrations, this probability has a Gaussian-like distribution around the atomic location [23]. Sampling plane by plane and tracking the changes that occurred due to the projectile-atoms interactions, the depth-sensitive channeling spectra as well as projectileflux density distributions may be obtained. In this paper, we are focused on structural changes observed near a dislocation edge. For that purpose, we investigated a single edge dislocation inside the Ni structure after its implementation in two ways:

Single edge dislocation inside Ni structure -MD-LAMMPS implementation
The first approach was performed by sampling based on the distances between projectiles and equilibrium positions determined by MD. For this purpose, we developed a model with two lattice half-planes cut off from the sample starting from the middle of it (100 nm) along the Z direction, until the end of the sample (see Figure 2 (a) and (b)). Although Ni is a cubic (fcc) structure, the atoms lay on two distinguished planes (cube base and the centre of the cube). Simply getting rid of only one of the planes in <001> direction will not preserve the long-range crystal symmetry and will create unsymmetrical defects and clustering. It seems that the results sufficiently agree with the Peierls-Nabarro model of an edge dislocation and can be fitted using the arctan function as follows: where: 0, 1, 2, 3 are fitting parameters. Among them, p0 corresponds to the distance between the asymptotes (denoted H in Figure 2(d)), while p2 and p3 correspond to the (x,y) coordinates of the inflexion point of the arctan function. The fitting procedure was based on a derivativefree Nelder-Mead method [24] implemented in Octave v.4.2.2 software package [25]. Moreover, we can extract interesting parameters regarding, i.e., how the amplitude changes as a function of the distance from the dislocation core (D = ⋅ 0 , see Eq. 1). What we obtained is a parabolic decrease of the arctan amplitude D (Figure 2 (d)). Such parameters can be used to improve or refine the ones used inside the program or to improve the way of describing extended defects inside the investigated structure.   [18,19,26]. To assess distortion created by the defect, two ion-flux density maps were produced; at the depth before and after the dislocation edge. The first set of maps (Figure 3) was produced using a custom crystallization procedure that is implemented into the McChasy 2.0 code. Two orientations of the dislocation are considered, namely inward and outward (positive and negative regarding the slip plane position), as shown in Figure 3 (a) and (b). Both have their dislocation lines (edges) as well as Burger vectors laying in the same plane parallel to the surface. In both cases, the edge was chosen to be placed at 100 nm from the sample surface. However, in the first case, the sample's surface is located in the region, at which two additional half-planes were added (inward or positive dislocation), and opposite in the other case (outward or negative dislocation).   Figure 3 Cross-sectional view along Z direction of two rows inserted along a/2<100> inside the Ni monocrystalline structure (a),(b). Cross-sections along X direction giving atomic positions taken shortly before ((c),(e)) and after dislocation core ((d),(f)). Corresponding ion flux densities maps near edge dislocation ((g)-(i)).
Because a single channeling projectile trajectory is usually defined by a sequence of simultaneous projectile-atom (far) collisions, the maps can be produced precisely (with 1 nm depth resolution) only with a large number of projectiles (ca. 10 7 ). The projectiles momentum is determined based on the impact parameter deviations.

Conclusions
The current state of the McChasy 2.0 code testing phase and the recent results regarding edge dislocations formed using manual, McChasy internal procedures, and using MD simulations are presented. An influence on the ion flux densities and atomic positions close to the dislocation edges are investigated. The presented results will be investigated further using different edge dislocation orientations like a/2 <110> {111} that are typically observed in the fcc systems [27]. Moreover, the merge of MD and MC simulations should be also performed for modelled e.g. dumbbells reported elsewhere [15,28] and/or for different types