Identifying Interstitials in MD Simulations — Max Space Clustering Method

Molecular Dynamics (MD) simulations are suitable to study the details of atomistic processes. Identification of interstitials is one of the main difficulties in the analysis of defect dynamics. To identify the interstitials, different methods use different criteria. All the methods, except the Wigner-Seitz method, assume some input parameters which may lead to an erroneous number of interstitials. We describe an unsupervised clustering algorithm, called Max Space clustering method, to identify interstitials without assuming any input parameters. We show that the algorithm is suited for identifying interstitials at high temperatures when amplitude of atomic oscillations become larger, a fact ignored by other models.


Introduction
Safe and successful operation of any nuclear energy system greatly relies on the development of suitable component materials. Changes in the material micro-structure due to high neutron doses will lead to material property changes. This has implications for the material lifetime in reactors. This has led to increased interest in both the experiments and simulations to understand the changes in micro-structure due to radiation damage [1].
Multi-scale methods are necessary to span the length and time scales involved in radiation damage due to neutron irradiation [2,3,4]. Neutrons create energetic primary knock-on atoms (PKAs) that can cause localized collision cascades. Interstitials and vacancies formed in the collision cascades diffuse and recombine to either neutralize each other or to form interstitial and vacancy clusters. This diffusive-recombination can occur throughout the operational lifetime of the reactor material, leading to significant changes in the physical and mechanical properties of materials.
In contrast to the time and length scales required for the observable radiation induced mechanical property changes, the primary radiation damage event that initiates these changes is smaller by many orders of magnitude [3]. Molecular dynamics (MD) is an atomistic simulation method that can provide details of atomistic processes generally upto a few nanoseconds [5,6,7]. Techniques such as Kinetic Monte Carlo (KMC) [8,9] or Mean field rate theory [10] are required to bridge the gap between MD and experimentally observable time and length scales. For example, the diffusion of interstitials, thier migration energy and details of thier diffusion can be obtained from MD and used in KMC simulations at the higherscales [8,9,11].
Interstitialcy diffusion -wherein an interstitial displaces a lattice atom, thereby making the lattice atom an interstitial -plays a prominent role in the diffusive recombination of Frenkel pairs created during the neutron irradiation. The identity of the diffusing atom changes continuously. Moreover, the interstitial can be part of dumbbell and crowdion configurations [3,7]. Therefore, identifying the changing interstitial and its trajectory from MD simulations is a complicated task. We have developed a robust and computationally efficient algorithm called Max-Space clustering method to identify the probable interstitials from MD simulations. Then we use graph optimization technique to identify the interstitial from dumbbell and crowdion configurations thereby creating a diffusion trajectory for the interstitial. In Section 2, we describe existing methods to identify interstitials and thier drawbacks. In Section 3, we describe our MD simulations of Cu,Fe and W interstitialcy diffusion. In Section.4, we describe the Max-Space clustering method and the graph optimization technique. We present our results and carry out a comparative study of existing methods of identifying interstitials with our Max-Space Clustering (MSC) method.

Identification of interstitials-Existing methods
MD has become a powerful tool to study the atomistic processes. Though numerous methods are present for identification of interstitials, we perform a comparative study of few commonly used methods.

Existing Methods
In the course of simulation, usually one has to follow any one of the following available criteria to identify the point defects for visualization and analysis: (i) Potential energy, (ii) Atomic Bonding, (iii) Geometric Methods: (a) Wigner-Seitz cell method (b)Lindemann's sphere (LS) method (c) The nearest neighbour sphere (NNS) method (d) The effective sphere (ES) method.
At very low temperatures, perhaps the simplest way of identifying the interstitials is to consider the atoms with a potential energy significantly higher (typically 2eV) than the equilibrium value as interstitials [12,13]. Though this method is very simple to implement, it can not provide details about the defect structures. When atomic bonding is considered, if an atom does not have correct coordination number of bonds, it is recognized as an interstitial [14].
Using the geometric methods, defect analysis is done by the comparison of atomic positions with respect to a geometric structure such as Wigner-Seitz cell [15], or LS [16,17,18], or NNS [19] or ES [3]. These geometric structures are centered on ideal lattice sites. Wigner-Seitz cell is a Voronoy polyhedron centered on an ideal lattice site. The LS method uses a sphere with a radius equal to atomic mean squared displacement at the melting point of the crystal being considered. The NNS method proposed by M. J. Caturla et al. uses a sphere with radius equal to half of the nearest neighbour distance r N N , whereas the ES method proposed by R. E. Stoller employs a sphere with radius equal to 30% of the inter-atomic separation of the crystal being studied [3]. Stoller, in his recent paper, also used a sphere with radius equal to 0.15 times the lattice parameter before ballistic phase and 0.25 times the lattice parameter after the ballistic phase to study the fragmentation of cascades into sub-cascades [20].
In geometric methods, if the geometric structure contains no atoms, the lattice postion is labelled as a vacancy and either if it contains an additional atom, or if an atom is not associated with any geometric structure, then it is labelled as an interstitial. Nowadays, perhaps ES and NNS methods are the most commonly used methods. Hence we compare our results with these in section 5.

Drawbacks of existing Methods
Methods using potential energy, atomic bonding and LS overestimate number of interstitials even at lower temperatures. They do not provide any details of defect structure [6]. Wigner-Seitz cell method is computationally expensive despite of its advantage that it does not involve any arbitrary cut-off size as an input parameter. Though LS criteria overestimates the number of interstitials even at low temperatures, it could be the motivation behind the widely used NNS and ES methods. Both NNS and ES methods give a rough number of interstitials. Though ES method is the most successful method among the others and computationally faster than Wigner-Seitz cell method, it overestimates the number of interstitials at higher temperatures and also it can not correct the number of interstitials for dumbbell and crowdion configurations. While using the any one of the above mentioned criteria, results are affected by the thermal vibration of atoms in the simulation cell. To eliminate the thermal vibrations several techniques can be used, [21,22,23]. But they are very computer intensive and cannot be employed to each configuration obtained from the simulation.
To overcome the above mentioned difficulties, we present a rigorous unsupervised clustering algorithm, called Max-Space Clustering (MSC) [24,25]. Apart from the MSC algorithm, we also use graph optimization technique to correct the number of interstitials for dumbbell and crowdion configurations and to get a smooth trajectory. These are described in Section 4.

MD Simulations
We use Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) code [26] for the classical MD study of interstitialcy diffusion in Cu, Fe and W single crystals. A cubic simulation box of 10 × 10 × 10 unit cells is used for each element. A single atom is introduced close to the center of the crystal. It causes an abrupt increase in pressure and leads to the deformation of the crystal. To overcome this,energy minimization of the system is performed till the potential energy of the crystal becomes minimum. This is done by iteratively adjusting the atom coordinates. The EAM potential by Foiles et al. [27] was used to simulate Cu. The EAM/alloy potential by Zhou et al. was used to simulate W [28]. For Fe, the concentration dependent EAM potential by Stukowski et al. [29] was used.
An NPT ensemble was used to carry out MD simulations at 0 bar pressure with periodic boundary conditions in all the three directions till the system equilibrates. After the equilibration, an NVE ensemble is used to study interstitialcy diffusion for 10ns. MD simulations were carried out at intervals of 100 K in the temperature range 300K-1100K and 300K-1200K for Cu and Fe respectively and at the intervals of 200 K in the temperature range 300K-2900K for W. The chosen temperature range is so as to have a valid potential well before the melting point. A time step of 1 femto-second was used for both, NPT and NVE ensembles. The atomic co-ordinates were output for every 100 time steps during the NVE simulations to study the interstitial transport. Hitherto, each output for every 100 time-steps is referred as a frame and the displacement of an atom from its nearest lattice position is referred to as offset.

Max-Space Clustering -Algorithm
Identification of interstitials in the analysis of defect dynamics is an issue of critical importance [6,22,23]. In this section, we describe Max-Space clustering method, which can identify interstitials without assuming any input cut-off values. The only required input is the positions of all atoms in a given time-step. At any given time-step, the interstitial exists as a single interstitial atom (SIA) or as part of a dumbbell or crowdion configuration. It has to be identified in each successive MD frame to produce a smooth trajectory. Dumbbell and crowdion configurations are shown in Fig.1 and Fig.2. In dumbbell configuration, two atoms share one lattice position. In crowdion configuration, three or more atoms share two or more lattice positions such that n atoms share n − 1 lattice positions. The MSC algorithm helps us to seperate these probable interstitials from the vibrating lattice atoms from low temperatures to the temperatures close to the melting point.
The MSC algorithm to identify the probable interstitials in a given MD time-step is as follows: • Find out the offsets of each atom in the frame. • Sort the atoms by their offsets in descending order.
• Find the difference in offsets between successive members of the sorted list of atoms.
• Find mid-point of successive offsets with the maximum difference. This is hitherto referred to as "threshold offset". • Group the atoms above and below the mid point of successive offsets with the maximum difference.
Two groups of atoms are created by the above algorithm as shown in Fig3. One group having one or more number of atoms with large offsets (probable interstitials) and the other group of atoms having smaller offsets (lattice atoms). The second group of atoms with smaller amplitudes also includes another sub-group of atoms with slightly higher amplitudes. This subgroup of atoms are the strongly displaced atoms nearby the interstitials and can be seen in Figs.1 & 2 and they are indicated by squares in Fig.3. The atoms in the group having larger offsets are all probable candidate interstitials.

Identifying the true interstitial
The graph optimization algorithm to identify the interstitialcy diffusion path is as follows: • The positions of probable interstitials obtaned from MSC algorithm, in each MD frame, are considered to be nodes of a graph structure [30]. • All the nodes across successive MD frames are connected thereby creating possible diffusion pathways. • A cost function is constructed for the different diffusion pathways such that it minimizes the distance between interstitials across consecutive frames and it maximizes the offset of interstitial. • The diffusion path having the lowest value for the cost function is chosen.
This technique is described more in detail in [24]   Figs.4-6 show that NNS method underestimates the number of interstitials and in some frames, it can not detect any interstitial even though it exists. ES method is ambiguous in choosing the threshold offset such as "15% or 25% or 30% of lattice parameter" [3,20] and does not provide any justification for the selected threshold offset. In this article we use "30% of lattice parameter" criterion of ES method for comparision. With the other criteria, number of interstitials predicted by the ES method further increases. Note that the ES method is the best amongst the "sphere methods". It however overestimates the number of interstitials as the temperature increases and it can not correct for dumbbell and crowdion configurations.
Though it is not a good idea to have a single threshold offset value for all the frames even at a given temperature, for a simple and most efficient way of finding interstitials, threshold offset separation values at 300K and the slope of the separating lines for Cu, Fe and W are given in table 1. Using these two parameters, threshold offsets at any temperature can be calculated.   Fig.7 shows the interstitialcy diffusion trajectory of Cu using the MSC and the graph optimization algorithm. These methods are being used to obtain interstitialcy diffusion parameters from MD simulations and are being used in KMC simulations [11]. We conclude by stating that these methods can be used for multiple interstitial diffusion tracking in collision cascades.