Nucleation of Nano-sized Prismatic Dislocation Loop from Spherical Vacancy Clusters in α-iron: An Atomic-scale Study

Abstract

The concentration of vacancy-type defects in α-iron is increased by plastic deformation and the presence of hydrogen. This leads to the accumulation of monovacancies, either in the form of planar vacancy clusters (VCs) or small voids. The prismatic dislocation loop (PDL) can nucleate from VC as vacancies agglomerate into two-dimensional (2D) VCs and collapse due to attractive force between two interior surfaces. The transition between 2D-VC and PDL is comparatively more straightforward, requiring only a short displacement without the need for atom diffusion to reach stability. However, the most stable VC configuration is three-dimensional (3D) (spherical cluster), which have lower formation energy than 2D-VCs. Despite their stability, the transformation from 3D-VC to PDL is complex, involving the diffusion of multiple atoms. A quantitative energy barrier is established for transitioning from 3D-VC to nano-sized 1/2<111> PDLs using an approach that combines the reaction rate theory and molecular dynamics (MD) simulations. The nucleation of PDL from a spherical cluster composed by 15 vacancies is a rare event at room temperature, even under considerable compressive strain since the activation energy is 1.33 eV. In contrast, 2D-VC with 37 vacancies can be nucleated to PDL with an energy barrier of 0.61 eV.

1. Introduction

Hydrogen (H) holds significant potential as a sustainable energy carrier and source, playing a crucial role in meeting the future energy demand.^1,2) Mainly due to the sensitivity of most metallic materials to hydrogen embrittlement (HE),^3,4,5) H production, storage, and transportation still present significant challenges. HE is characterized by the deterioration of the mechanical properties of susceptible materials, mainly loss of ductility due to H.^6,7) HE has been an essential issue particularly severe in steels because of the high mobility of H in iron (Fe).⁷⁾ α-iron is an ideal choice for studying HE due to the absence of impurities that make empirical results and numerical simulation more complex to interpret. Several theories have been suggested to elucidate the underlying causes of the embrittlement phenomenon. Among these three mechanisms, hydrogen-enhanced decohesion (HEDE),⁸⁾ hydrogen-enhanced localized plasticity (HELP),⁹⁾ and hydrogen-enhanced strain-induced vacancy (HESIV),¹⁰⁾ researchers have shown a particularly strong interest. The HESIV mechanism argues that the observed damage is attributed to abundant vacancies. Experimental evidence supporting this theory includes H outgassing from the specimen, which resulted in incomplete recovery of the mechanical properties to their H-free state. This finding implies that H is not fully responsible for the fracture or damage process. Instead, its role lies in generating excess vacancies, which act as the primary mechanisms of damage. As a result, the fracture process is significantly influenced by vacancies, representing a crucial and decisive factor from this perspective.^11,12)

Moreover, the recent study investigated the relationship between local H concentration and vacancy formation energies, focusing on the observed phenomenon of superabundant vacancy formation.¹³⁾ It was found that regions characterized by high hydrogen concentrations demonstrated significantly reduced vacancy formation energies, providing support for the experimentally observed superabundance of vacancies. Further, the observation indicates that the presence of H atoms enhances the clustering of vacancies. Additionally, observations have revealed that H is released through aging in an H-free environment at room temperature, causing monovacancies to aggregate into clusters or disappear.¹⁴⁾ The recent study discusses the accumulation of high-density vacancies near stress singularities (i.e., triple junction of grain boundaries or crack tip), which results from high-speed vacancy transport facilitated by prismatic dislocation loop (PDL) diffusion.¹⁵⁾ The study mainly investigates the transition behavior between two-dimensional (2D) VCs and PDLs under external stress. Forming a vacancy-type PDL from 2D-VC is relatively explicit, as atoms on the top and bottom sides of the VC displaced into each other due to attractive forces between opposing atoms.¹⁵⁾ However, transforming from a three-dimensional (3D) void to a vacancy-type PDL involves diffusion and movement of many atoms, making it more challenging.¹⁶⁾ The study investigated the defects resulting from the collapse of vacancy clusters in face-centered cubic (FCC) and body-centered cubic (BCC) metals. Fikar et al. investigate a void structure characterized by an almost spherical shape containing a small number of vacancies in their study. The void collapse to form a stacking-fault tetrahedron in FCC crystals. The study suggested that void collapse into vacancy loops with a [001] Burgers vector form on (001) planes in BCC iron. However, the transition was observed by giving hydrostatic stress by reducing the lattice parameter by 5–10%.¹⁷⁾

In this study, the activation energy of the transition from a 3D-VC to a PDL is clarified in comparison with the transition from a 2D-VC, where activation energy is already evaluated by Nudged Elastic Band (NEB) calculation. Section 2 shows the simulation models and calculations methods. The procedure to obtain the 2D and 3D (and spherical) VCs, considering their minimum formation energy is also explained. In Section 3, the study analyzes the formation energies of defects, including 2D and 3D VCs and hexagonal 1/2<111> PDLs. Section 4 explores the activation energy associated with the transitions from planar and spherical VCs to the 1/2<111> PDL.

2. Simulation Model and Methodology

2.1. Simulation Method and Interatomic Potentials

The computational analysis presented in this study employed a rectangular simulation box with periodic boundary conditions (PBC) in all directions. The analyses were performed with the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS), open-source software for molecular dynamics (MD) simulation.¹⁸⁾ To describe the Fe–Fe interactions, an embedded atom method (EAM) potential developed by Wen was used in this study, which was developed using the Fe–Fe model of Ackland-2004 empirical potential.^19,20) The Ackland-2004 potential has been approved to satisfactorily represent the properties of dislocations and dislocation loops in iron.²¹⁾

2.2. Introducing the Planar and Spherical Vacancy Clusters

As shown in Fig. 1, the simulation model was created with α-iron single crystal containing N_Fe=1710 720 Fe atoms with dimensions of 25.18 nm, 24.24 nm, and 32.64 nm along x, y, and z directions, with [211], [011] and [111], respectively. This model represents a good trade-off between avoiding spurious interaction with periodic images and an affordable computational cost.

Fig. 1. Simulation model showing, (a) crystallographic orientation (●, ■, and ▲ indicates atoms in the different layers in [111]) and three-layer V_2D37 (○, □, and △ to express the vacancy in the different layers), (b) schematic of V_3Dn. (Online version in color.)

The planar VCs with n vacancies will be denoted as V_2Dn, while the spherical VCs with n vacancies will be denoted as V_3Dn hereafter. To create the planar VC, atoms are removed from the (111) plane. During the preparation of V_2D37 for the (111) plane, three layers were considered while also considering the symmetrical hexagonal shape, which is the most stable configuration, as illustrated in Fig. 2. The selection of V_2D37 in this study is based on its feature as the smallest symmetrical VC with a tendency for nucleating into a PDL, as established through prior findings. After removing the atoms to create VC, the atomic positions are relaxed using the conjugate gradient (CG) method, and the formation energy of VC then calculated. The formation energy E_f,n of V_2Dn or V_3Dn is computed using

  

$$E_{\text{f},\mathrm{\hspace{0.17em} }n}=E_n-(N_{\text{Fe}}-n)\times E_{\text{Fe}}$$(1)

where E_n is the relaxed total energy of the model with V_2Dn or V_3Dn and E_Fe is the energy of an Fe atom in bulk perfect crystal.

Fig. 2. Structure of V_2D37 (a) ●, ■, and ▲ indicates atoms in the different layers in crystal structure (b) ○, □, and △ to express the vacancy in the different layers.

The spherical clusters were produced by selecting a specific region using a defined radius (R), as shown in Figs. 1(b) and 3. The spherical region comprised N atoms, and the center of the spherical region was set on an atom to maintain symmetry. The vacancies were introduced by removing atoms one by one, starting with the most stable VC with vacancy site (i), determined by selecting the atom with the lowest formation energy (see Eq. (1)) from (N−i) sample atoms. The stable V_3D15 were introduced using our in-house python-based algorithm. The detail of stable structure until V_3D15 is illustrated in Fig. 4. The inclusion of V_3D15 in this study is motivated by its representation as the smallest 3D-VC exhibiting spherical cluster characteristics and its ability to provide enhanced symmetry due to its minimal structural configuration. As illustrated in Fig. 4(a), there are three layers in the [011] direction where vacancies are introduced and labeled as I, II, and III. In the [111] direction, the layers are denoted by the characters A–G (see Fig. 4(b)). Also, Fig. 4(c) presents a comprehensive 2D representation of the vacancy sequence, employing the layer labels for a clear correlation with the 3D perspective utilizing the [011] and [111] direction layers.

Fig. 3. The process of creating a larger spherical cluster, (a) the selection volume used to obtain suitable atomic coordinates, and (b) establishing a stable spherical cluster by defining the number of vacancies based on the formation energy. (Online version in color.)

Fig. 4. Structure of V_3D15 (a) Top view; ●, ■, and ▲ indicates atoms in the different layers in crystal structure, (b) Side view; ○, □, and △ to indicates the vacancies in the different layers, (c) Side view; express the sequence of the introduced vacancies to form V_3D15. (Online version in color.)

2.3. Calculation Conditions for Determining Activation Energy for VC to PDL Transition

Two distinct transitions were observed in this study. Firstly, the transition between V_2D37 and PDL was investigated. Secondly, the transition from V_3D15 was examined. The activation energies for the transition between 2D-VC and PDL using NEB were presented in our previous study. However, this study considers the transition from V_3D15. In the case of V_3D15, they involve atomic diffusion, making finding the saddle points complex. Therefore, in this study, we used an approach combining the reaction rate theory and MD simulations to calculate the transition activation energy. Multiple MD simulations are performed with various initial atomic velocities sampled from the Maxwell-Boltzmann distribution for the same initial atomic structures. Notably, each case was subjected to ~50 independent MD simulations (replicas) for a given temperature and strain along [111](z axis). The objective is to gather statistical distribution data on the probability of successive transitions occurring and their survival probabilities. The simulations were performed with the canonical (NVT) ensemble and the Berendsen barostat, allowing the system to equilibrate quickly for the given temperature with constant pressure. Consequently, after equilibration, MD steps were considered to calculate the PDL nucleation time. The simulations were performed using time step of 1 fs.

The activation energy for the transition from V_2D37 to PDL was determined without external stress at temperatures 315 K and 330 K. In contrast, for the V_3D15 to PDL transition, a strain equivalent to 4% and 5% at temperatures 1100 K and 1150 K, respectively, was applied to achieve the transition within the simulation time. The transition from V_2D37 to PDL was determined by analyzing the displacement of middle atoms on the top and bottom faces of VC (see Fig. 5). Pairwise inter-distances between these top and bottom faces were considered. When the interspacing between these pairs becomes shorter than ~1.3b, where b is the Burgers vector, it was considered as PDL. These conditions were applied based on the transition state obtained from the NEB calculations.¹⁵⁾ Also, the PDL structure is confirmed by DXA modifier.²²⁾ However, in the case of V_3D15, using the same approach used for V_2D37 to observe the transition is challenging. This discrepancy arises because while the 2D-VC structure maintains symmetry, the 3D-VC structure does not exhibit the same level of symmetry. Consequently, observing nucleation, similar to the 2D-VC case, becomes complex due to challenges in tracking the pairwise inter-distance within the 3D structure. However, the V_3D15 to the V_2D15 (approximately 2D) structure was decided based on a similar approach as explained for V_2D37 to PDL (to reduce the calculation cost), and later, DXA modifier was used to judge the PDL.

Fig. 5. Method used to judge the transition from VC to corresponding PDL. (Online version in color.)

3. Formation Energy of V_2Dn and V_3Dn

Figure 6 illustrates the formation energy for different V_2Dn, PDL, and V_3Dn sizes. The formation energy of the V_2Dn approximately linearly increases with the number of vacancies. The V_3Dn have a substantially lower formation energy than the V_2Dn. Furthermore, the difference increases as the number of vacancies increases. The lowest-energy configurations of V_2D37 and V_3D15 are displayed in Figs. 2 and 4, respectively. From its symmetricity, V_3D15 is quite stable compared to different sizes.

Fig. 6. The figure illustrates a comparison of the formation energies of V_2Dn, V_3Dn, and PDL corresponding to V_2Dn. (Online version in color.)

For larger VCs (n≥37), PDLs are stable than 2D-VCs. However, 2D-VC is more stable than PDL for n=19. The extrapolation from this trend gives more energy difference for n<15. 3D-VCs are more stable than 2D-VCs except very small VCs. These results give the prediction that the transition to PDL is quite difficult for V_3D15. However, when PDL becomes too small, maintaining the hexagonal symmetry is difficult and possibly have different structure and different energy trend.

4. Energy Barrier for the Transition from V_2D37 and V_3D15 to PDL

The rate of PDL nucleation, according to the reaction rate theory, can be expressed as a first-order rate process:²³⁾

  

$$\dot{N}=(N_0-N)v\text{exp}\left[\raisebox{1ex}{$-G_{\text{ac}}(\epsilon ,\mathrm{\hspace{0.17em} }T)$}\!\left/ \!\raisebox{-1ex}{$k_{\text{B}}T$}\right.\right]$$(2)

Where (N₀−N) represent the available nucleation site at given time t, v is attempt frequency, G_ac(ε, T) is activation free energy, ε is applied strain, k_B is Boltzmann constant and T is absolute temperature.^24,25)

Let F be the survival probability for PDL nucleation. Let N₀ represent the total number of MD simulations, while N refers to the successive MD simulations that initiate PDL nucleation at time t, i.e., $F=\left(1-\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$N_0$}\right.\right)$. Thus, F is one at t = 0 and F approaches zero after long sampling time. The nucleation rate $\dot{N}$ for PDL can be defined as (i.e., $\text{d}N/\text{d}t$ dF/dt = $-\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$N_0$}\right.\text{d}N/\text{d}t$)

  

$$\raisebox{1ex}{$\dot{N}$}\!\left/ \!\raisebox{-1ex}{$(N_0-N)$}\right.=-\frac{\text{d}F/\text{d}t}{F}$$(3)

Where Eq. (3) can be written as

  

$$ln[F(t)]=-\int \raisebox{1ex}{$\dot{N}$}\!\left/ \!\raisebox{-1ex}{$(N_0-N)$}\right.dt$$(4)

Substituting Eqs. (2), (3), (4) it can be written as

  

$$ln[F(t)]=-\int v\text{exp}\left[\raisebox{1ex}{$-G_{\text{ac}}(\epsilon ,\mathrm{\hspace{0.17em} }T)$}\!\left/ \!\raisebox{-1ex}{$k_{\text{B}}T$}\right.\right]dt$$(5)

Solving Eq. (5), fraction of surviving entities can be expressed as

  

$$\raisebox{1ex}{$N$}\!\left/ \!\raisebox{-1ex}{$N_0$}\right.=1-\text{exp(}-At\text{)}$$(6)

Where $A=v\text{exp}\left[\raisebox{1ex}{$-G_{\text{ac}}(\epsilon ,\mathrm{\hspace{0.17em} }T)$}\!\left/ \!\raisebox{-1ex}{$k_{\text{B}}T$}\right.\right]$, the fitting parameter obtained using a statistical process of PDL nucleation at a given temperature. This is achieved by plotting fraction of surviving entities as a function of time (see Fig. 7). The activation free energy G_ac(ε, T) can be expressed as G_ac(ε, T)$=E_n^a-\epsilon \stackrel{\prime }{\mathrm{\Omega }}$, where $E_n^a$ is the activation energy at zero applied strain, and Ω is the activation volume, i.e., $\stackrel{\prime }{\mathrm{\Omega }}$=EΩ, where E is Young’s modulus. We assumed $E_n^a$ is temperature independent. By evaluating A for different temperature and applied strain, we can evaluate v, $E_n^a$, and Ω. Then we can predict the activation free energy G_ac at realistic small strain, and finally we can know the transition rate of rate event at room temperature.

Fig. 7. The fraction of surviving entities as a function of time. (Online version in color.)

In the case of V_2D37, the transition occurs in the absence of strain with the help of thermal activation (see Table 1). A snapshot of the nucleation is shown in Fig. 7. However, transitioning from V_3D15 to PDL is difficult at room temperature without strain (see Table 1). At 4% strain, the activation free energy is 1.34 eV and the attempt frequency is 8.025×10¹⁵. Activation free energy decreased to 1.33 eV, while attempt frequency was 13.72×10¹⁵ at 5% strain. The influence of the applied strain on the activation free energy is subtle and closely connected to the relatively low activation volume, Ω=8×10⁻³ $(\stackrel{\circ }{\mathrm{A}}{)}^3$. Accordingly, the activation energy under no strain conditions is obtained as 1.38 eV. This value is close the sum of activation energy of vacancy diffusion (~0.6 eV) and the vacancy binding energy to vacancy cluster (0.6–1.14 eV).^26,27) The small activation volume is also reasonable if it is related to vacancy diffusion because it is almost stress independent phenomenon. It is considered that the transition is controlled by vacancy diffusion for the appropriate VC shape to cause transition to PDL. Figure 8 illustrates the structure and associated formation energies of V_3D14, V_3D15 and V_3D16 with symmetrical nature (Fig. 8(a)) considered in this study, using our in-house algorithm, and 3D-VC with same n vacancies which deviate from the symmetrical nature (Fig. 9(b)). The formation energies increase by ~0.25 eV as the VC loses its symmetric nature. Also, ~0.6 eV is necessary to remove one vacancy from these clusters, and it agrees with the above mentioned binding energy. The nucleation process of 1/2<111> PDL from V_3D15 is shown in Fig. 9. The vacancy diffusion during the transition was observed using the Wigner-Seitz defect analysis modifier.²⁸⁾ This modifier can identify point defects within crystalline structures using the Wigner-Seitz cell method. It can count vacancies and trace their movement throughout the lattice (see Fig. 9(c)). Furthermore, the application of compressive strain controls the mode of diffusion mechanism, such as the diffusion direction. Compressive stress retards diffusivity in the stress direction; consequently, the applied compressive strain facilitates diffusion normal to the [111] direction, hence V_3D15 nucleated into 1/2<111> PDL.²⁹⁾ Previous research showed that in the presence of H, 2D-VC becomes more stable compared to 3D-VC. Therefore, transition from 3D to 2D VC is enhanced.³⁰⁾ Moreover, the transition between 2D to PDL is already established, and it is possible at room temperature, which implies that transition from 3D to PDL is possible at room temperature in the presence of H, from the viewpoint of formation energies (2D-VC more stable than 3D-VC in the presence of H). However, in this study, we did not treat the transition with H, which is left for future work.

Table 1. Energy barriers of VC–PDL transitions for V_2D37 and V_3D15.

Defect type	Activation energy (NEB)14)	Activation energy (This study, zero strain)
V2D37	0.68 eV	0.61 eV
V3D15	–	1.38 eV

Fig. 8. Structure of V_3D14, V_3D15, and V_3D16 and corresponding formation energies (a) symmetrical structure considering lowest formation energies using the in-house algorithm (b) structure deviating from symmetrical nature. (Online version in color.)

Fig. 9. Snapshots of the V_3D15 transition to PDL (a) Side view, (b) Top view, (c) Side view; vacancy diffusion during transition. Grey color represents the non-BCC atoms, while blue color indicates BCC atoms and red color represents vacancies. (Online version in color.)

5. Conclusion

This study used molecular statics simulation to investigate the formation energy of 3D (spherical-shaped) VCs and compared it with regular 2D (hexagonal-shaped) VCs, vacancy-type PDLs with Burgers vector 1/2<111> in α-iron. Furthermore, the reaction rate theory combined with MD simulation was used to calculate the energy barrier for VC and PDL transitions without/with applied strain. Accordingly, the main findings are summarized as follows:

(1) The transition between V_2D37 and PDL under zero external compressive strain conditions is possible at room temperature, which aligns with previous calculations.¹⁵⁾

(2) The transition between V_3D15 and PDL is difficult at room temperature since the activation energy is 1.38 eV. However, slight temperature increase or long-term aging enables the transition.

Acknowledgements

This research was supported by ISIJ Research Group II “Extraction of essential factor technology for hydrogen embrittlement evaluation” and JSPS KAKENHI Grant Number 19H02025.

References

• 1)   B.  Johnston,  M. C.  Mayo and  A.  Khare: Technovation, 25 (2005), 569. https://doi.org/10.1016/j.technovation.2003.11.005
• 2)   R. L.  Amaro,  N.  Rustagi,  K. O.  Findley,  E. S.  Drexler and  A. J.  Slifka: Int. J. Fatigue., 59 (2014), 262. https://doi.org/10.1016/j.ijfatigue.2013.08.010
• 3)   A. M.  Abdalla,  S.  Hossain,  O. B.  Nisfindy,  A. T.  Azad,  M.  Dawood and  A. K.  Azad: Energy Convers. Manag., 165 (2018), 602. https://doi.org/10.1016/j.enconman.2018.03.088
• 4)   K.  Mazloomi and  C.  Gomes: Renew. Sustain. Energy Rev., 16 (2012), 3024. https://doi.org/10.1016/j.rser.2012.02.028
• 5)   S.  Sharma and  S. K.  Ghoshal: Renew. Sustain. Energy Rev., 43 (2015), 1151. https://doi.org/10.1016/j.rser.2014.11.093
• 6)   M.  Nagumo: Fundamentals of Hydrogen Embrittlement, Vol. 921, Springer Singapore, (2016), 137. https://doi.org/10.1007/978-981-10-0161-1
• 7)   H. K. D. H.  Bhadeshia: ISIJ Int., 56 (2016), 24. https://doi.org/10.2355/isijinternational.ISIJINT-2015-430
• 8)   W. W.  Gerberich,  R. A.  Oriani,  M. J.  Lji,  X.  Chen and  T.  Foecke: Philos. Mag. A, 63 (1991), 363. https://doi.org/10.1080/01418619108204854
• 9)   H. K.  Birnbaum and  P.  Sofronis: Mater. Sci. Eng. A, 176 (1994), 191. https://doi.org/10.1016/0921-5093(94)90975-X
• 10)   M.  Nagumo: Mater. Sci. Technol., 20 (2004), 940. https://doi.org/10.1179/026708304225019687
• 11)   F. F.  Dear and  G. C.  Skinner: Philos. Trans. Royal Soc. A, 375 (2017), 20170032. https://doi.org/10.1098/rsta.2017.0032
• 12)   E.  Hayward and  C. C.  Fu: Phys. Rev. B, 87 (2013), 174103. https://doi.org/10.1103/PhysRevB.87.174103
• 13)   K.  Arakawa,  A.  Kageyama,  H.  Hiroshima,  H.  Yasuda and  S.  Ogata: ISIJ Int., 61 (2021), 2305. https://doi.org/10.2355/isijinternational.ISIJINT-2021-140
• 14)   L.  Chiari and  M.  Fujinami: ISIJ Int., 62 (2022), 832. https://doi.org/10.2355/isijinternational.isijint-2021-422
• 15)   M.  Vijendran and  R.  Matsumoto: Comput. Mater. Sci., 225 (2023), 112195. https://doi.org/10.1016/j.commatsci.2023.112195
• 16)   J.  Fikar,  R.  Schäublin,  D. R.  Mason and  D.  Nguyen-manh: Nucl. Mater. Energy, 16 (2018), 60. https://doi.org/10.1016/j.nme.2018.06.011
• 17)   C. C.  Matthai and  D. J.  Bacon: J. Nucl. Mater., 135 (1985), 173. https://doi.org/10.1016/0022-3115(85)90075-3
• 18)   A. P.  Thompson,  H. M.  Aktulga,  R.  Berger,  D. S.  Bolintineanu,  W. M.  Brown,  P. S.  Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott and S. J. Plimpton: Comput. Phys. Commun., 271 (2022), 108171. https://doi.org/10.1016/j.cpc.2021.108171
• 19)   M.  Wen: Comput. Mater. Sci., 197 (2021), 110640. https://doi.org/10.1016/j.commatsci.2021.110640
• 20)   G. J.  Ackland,  M. I.  Mendelev,  D. J.  Srolovitz,  S.  Han and  A. V.  Barashev: J. Phys. Condens. Matter, 16 (2004), S2629. https://doi.org/10.1088/0953-8984/16/27/003
• 21)   D.  Terentyev,  D. J.  Bacon and  Y. N.  Osetsky: J. Phys. Condens. Matter, 20 (2008), 445007. https://doi.org/10.1088/0953-8984/20/44/445007
• 22)   A.  Stukowski and  K.  Albe: Model. Simul. Mater. Sci. Eng., 18 (2010), 085001. https://doi.org/10.1088/0965-0393/18/8/085001
• 23)   K.  Krūkle-Bērziņa,  A.  Actiņš and  A.  Bērziņš: J. Math. Chem., 50 (2012), 2120. https://doi.org/10.1007/s10910-012-9999-4
• 24)   K. S.  Ng and  A. H. W.  Nga: Modelling Simul. Mater. Sci. Eng., 16 (2008), 055004. https://doi.org/10.1088/0965-0393/16/5/055004
• 25)   M.  Wen,  B.  An,  S.  Fukuyama,  K.  Yokogawa and  A. H. W.  Ngan: Carbon, 47 (2009), 2070. https://doi.org/10.1016/j.carbon.2009.03.065
• 26)   S. E.  Restrepo,  H.  Lambert and  A. T.  Paxton: Phys. Rev. Mater, 4 (2020), 113601. https://doi.org/10.1103/PhysRevMaterials.4.113601
• 27)   D.  Kulikov,  L.  Malerba and  M.  Hou: Philos. Mag., 86 (2006), 141. https://doi.org/10.1080/14786430500156930
• 28)   K.  Nordlund,  M.  Ghaly,  R. S.  Averback,  M.  Caturla,  T.  Diaz de la Rubia and  J.  Tarus: Phys. Rev. B, 57 (1998), 7556. https://doi.org/10.1103/PhysRevB.57.7556
• 29)   J. W.  Jang,  J.  Kwon and  B. J.  Lee: Scr. Mater., 63 (2010), 39. https://doi.org/10.1016/j.scriptamat.2010.02.045
• 30)   X. D.  Pan,  T.  Lu,  Y. M.  Lyu,  Y. P.  Xu,  H. S.  Zhou,  Z. S.  Yang,  G. J.  Niu,  X. C.  Li,  F.  Gao and  G. N.  Luo: J. Nucl. Mater, 542 (2020), 152500. https://doi.org/10.1016/j.jnucmat.2020.152500