Whatever the system considered and the code used, one needs to provide more inputs than just the atomic positions and types. Most codes suggest some values for these inputs. However, their tuning may still be necessary as default values may very well be suited for some supposedly standard situations and irrele­vant for others. Blind use of ab initio codes may thus lead to disappointing errors. Indeed, not all these choices are trivial, so mistakes can be hard to notice for the beginner. Choices are usually made out of experience, after considering some test cases needing small calculation time.

One can distinguish between choices that should be made only once at the beginning of a study and calculation parameters that can be tuned calculation by calculation. The main unchangeable choices are the exchange and correlation functional and the pseudopotentials or PAW atomic data for the various atomic types in the calculation.

First, one has to choose the flavor of the exchange and correlation functional that will be used to describe the electronic interactions. Most of the time one chooses either an LDA or a GGA func­tional. Trends are known about the behavior of these functionals: LDA calculations tend to overestimate the bonding and underestimate the bond length in bulk materials, the opposite for GGA. However, things can become tricky when one deals with defects as energy differences (between defect-containing and defect-free cells) are involved. For insulating materi­als or materials with correlated electrons, the choice of the exchange and correlation functional is even more difficult (see Section 1.08.5).

The second and more definitive choice is the one of the pseudopotential. We do not mean here the choice of the pseudoization scheme but the choice of the pseudopotential itself. Indeed, calculated ener­gies vary greatly with the chosen pseudopotential, so energy differences that are thermodynamically or kinetically relevant are meaningless if the various calculations are performed with different pseudopo­tentials. The determination ofthe shape ofthe atomic basis set in the case of localized bases is also of importance, and it is close in spirit to the choice of the pseudopotential except that much less basis sets than pseudopotentials are available.

More technical inputs include

• the k-point sampling. The larger the number of k points to sample the Brillouin zone, the more accurate the results but the heavier the calcula­tions will be. This is especially true for metallic systems that need fine sampling of the Brillouin zone, but convergence with respect to the number of k points can be accelerated by the introduction of a smearing of the occupations of electronic levels close to the Fermi energy. The shape and width of this smearing function is then an addi­tional parameter.19

• the number ofplane waves (obviously for plane wave codes but also for some other codes that also use FFT). Once again the larger the number of plane waves, the more accurate and heavier the calculation.

• the convergence criteria. The two major conver­gence criteria are the one for the self-consistent loop of the calculation of the ground-state electronic wave functions and the one to signal the convergence of a relaxation calculation (with some threshold depending on the forces acting on the atoms).

Dislocations in the dominant slip system in bcc metals have (111)/2 Burgers vectors and {110} slip planes. Here, we consider an edge dislocation with Burgers vector b = 1/2 [111] (along x-axis), slip plane normal [110] (along j-axis), and line direction [T 12] (along z-axis). To prepare the atomic configuration, we first create a perfect crystal with dimensions 30[111], 40[H0], 2[ 112] along the x-, y-, z-axes. We then remove one-fourth of the atomic layers normal to the y-axis to create two free surfaces, as shown in Figure 8(a).

We introduce an edge dislocation dipole into the simulation cell by displacing the positions of all atoms according to the linear elasticity solution of the displacement field of a dislocation dipole. To satisfy PBC, the displacement field is the sum of the contributions from not only the dislocation dipole

<№

(a)

Figure 8 (a) Schematics showing the edge dislocation dipole in the simulation cell. b is the dislocation Burgers vector of the upper dislocation. Atoms in shaded regions are removed. (b) Core structure of edge dislocation (at the center) and surface atoms in FS Ta after relaxation visualized by Atomeye.41 Atoms are colored according to their central-symmetry deviation parameter. Adapted from Li, J. In Handbook of Materials Modeling;Yip, S., Ed.; Springer: Dordrecht, 2005; pp 1051-1068; Mistake-free version at http://alum. mit. edu/www/liju99/Papers/05/Li05-2.31.pdf; Kelchner, C. L.; Plimpton, S. J.; Hamilton, J. C. Phys. Rev. B 1998, 58, 11085.

inside the cell, but also its periodic images. Care must be taken to remove the spurious term caused by the conditional convergence of the sum.26,40-42 Because the Burgers vector b is perpendicular to the cut — plane connecting the two dislocations in the dipole, atoms separated from the cut-plane by <|b|/2 in the x-direction need to be removed. The resulting struc­ture contains 21414 atoms. The structure is subse­quently relaxed to a local energy minimum with zero average stress. Because one of the two dislocations in the dipole is intentionally introduced into the vacuum region, only one dislocation remains after the relaxa­tion, as shown in Figure 8(b).

The dislocation core is identified by central symmetry analysis,13 which characterizes the degree of inversion-symmetry breaking. In Figure 8(b), only atoms with a central symmetry deviation (CSD) parameter larger than 1.5 A2 are plotted. Atoms with CSD parameter between 0.6 and 6 A2 appear at the center of the cell and are identified with the disloca­tion core. Atoms with a CSD parameter between 10 and 20 A2 appear at the top and bottom of the cell and are identified with the free surfaces.

The edge dislocation thus created will move along the x-direction when the shear stress axy exceeds a critical value. To compute the Peierls stress, we apply shear stress axy by adding external forces on surface atoms. The total force on the top surface
atoms points in the x-direction and has magnitude of Fx = axyLxLz. The total force on the bottom sur­face atoms has the same magnitude but points in the opposite direction. These forces are equally distributed on the top (and bottom) surface atoms. Because we have removed some atoms when creating the edge dislocation, the bottom surface layer has fewer atoms than the top surface layer. As a result, the external force on each atom on the top surface is slightly lower than that on each atom on the bottom surface.

We apply shear stress axy in increments of 1 MPa and relax the structure using the conjugate gradient algorithm at each stress. The dislocation (as identified by the core atoms) does not move for axy < 27 MPa but moves in the x-direction during the relaxation at axy = 28 MPa. Therefore, this simulation predicts that the Peierls stress of edge dislocation in Ta (FS potential) is 28 ± 1 MPa. The Peierls stress computed in this way can depend on the simulation cell size. Therefore, we will need to repeat this calculation for several cell sizes to obtain a more reliable prediction of the Peierls stress. There are other boundary conditions that can be applied to simulate disloca­tions and compute the Peierls stress, such as PBCs in both x — and j-directions,42 and the Green’s function boundary condition.44 Different boundary conditions have different size dependence on the numerical error of the Peierls stress.

The simulation cell in this study contains two free surfaces and one dislocation. This is designed to minimize the effect of image forces from the bound­ary conditions on the computed Peierls stress. If the surfaces were not created, the simulation cell would have to contain at least two dislocations so that the total Burgers vector content was zero. On application of the stress, the two dislocations in the dipole would move in opposite directions, and the total energy would vary as a function of their relative position. This would create forces on the dislocations, in addi­tion to the Peach-Koehler force from the applied stress, and would lead to either overestimation or underestimation of the Peierls stress. On the contrary, the simulation cell described above has only one dislocation, and as it moves to an equivalent lattice site in the x-direction, the energy does not change due to the translational symmetry of the lattice. This means that, by symmetry, the image force on the dislocation from the boundary conditions is iden­tically zero, which leads to more accurate Peierls stress predictions. However, when the simulation cell is too small, the free surfaces in the y-direction
and the periodic images in the x-direction can still introduce (second-order) effects on the critical stress for dislocation motion, even though they do not produce any net force on the dislocation.

Initial work of Bacon and coworkers indicated that the number of stable displacements remaining at the end of a cascade simulation, ND, exhibited a power — law dependence on cascade energy.84 For example, their analysis of iron cascade simulations between 0.5 and 10 keV at 100 K showed that the total number of surviving point defects could be expressed as

Nd = 5-67E MD79 [3]

where EMD is given in kiloelectronvolts. This rela­tionship is not followed below about 0.5 keV because true cascade-like behavior does not occur at these
low energies. Subsequent work by Stoller64-67 indi­cated that Nd also begins to deviate from this energy dependence above 20 keV when extensive subcascade formation occurs. This is illustrated in Figure 9(a) where the values of ND obtained in cascade simula­tions at 100 K is plotted as a function of cascade energy. At each energy, the data point is an average of between 7 and 26 cascades, and the error bars indicate the standard error of the mean. It appears that three well-defined regions with different energy dependencies exist. A power-law fit to the points in each energy region is also shown in Figure 9(a). The best-fit exponent in the absence of true cas­cade conditions below 0.5 keV is 0.485. From 0.5 to

(b)

ratio at low energies was first measured in 1978 and

is well known.57,85 The error bars again reflect the standard error and the dashed line through the points is only a guide to the eye. The MD/NRT ratio is greater than 1.0 at the lowest values of EMD, indicat­ing that the NRT formulation underestimates defect production in this energy range. This is consistent with the low-energy (near threshold) simulations preferentially producing displacements in the ‘easy’ directions.[8] [9] [10] [11] [12] [13] [14] The actual displacement threshold var­ies with crystallographic direction and is as low as ~19 eV in the [100] direction.20,84 Thus, using the recommended average value of 40 eV £d in eqn [2] predicts fewer defects at low energies. The aver­age value is more appropriate for the higher energy events where true cascade-like behavior occurs. In the cascade-dominated regime, the defect density within the cascade increases with energy. Although many more defects are produced, their close proxim­ity leads to a higher probability of in-cascade recom­bination and a lower defect survival fraction.

The surviving defect fraction shows a slight increase as the cascade energy increases above 20, and the indicated standard errors make it arguable that the increase is statistically significant. If signifi­cant, the increase appears to be associated with sub­cascade formation, which becomes prominent above 10-20 keV. In the channeling regions between the high-angle collisions that produce the subcascades shown in Figures 7 and 8, the moving atom loses energy in many low-angle scattering events that pro­duce low-energy recoils. These are essentially like low-energy cascades, which have higher-than-average defect survival fractions (Figure 9). These events could contribute to the incremental increase in defect sur­vival at the highest energies. The average defect sur­vival fraction of ^0.3 NRT shown for cascade energies greater than about 10 keV is consistent with values of Frenkel pair formation obtained from resistivity change measurements following low-temperature

neutron irradiation and ion irradiation.26,27,57,85

The effect of irradiation temperature is shown in Figure 10, which compares the defect survival fractions obtained from simulations at 100, 600, and 900 K. Although it is difficult to discern a consistent effect of temperature between the 600 and 900 K data points, the defect survival fraction at 100 K is always somewhat greater than at either of the two higher temperatures. A similar result for iron was reported in Bacon and coworkers.84 In addition to an interest in radiation temperature itself, the effect of temperature is relevant to the
simulations presented here because no thermostat was applied to the simulation cell to control tem­perature. As mentioned above, the energy intro­duced by the PKA will lead to some heating if the simulation cell temperature is not controlled by a thermostat. For example, in a 1 keV cascade simula­tion with 54 000 atoms, the average temperature rise will be about 140 K when all the kinetic energy of the PKA is distributed in the system. This change in temperature should be more significant at 100 K than at higher temperatures. The fact that defect survival at 600 and 900 K is lower than at 100 K suggests that the 100 K results may be

somewhat biased toward lower survival values by the PKA-induced heating. This is in agreement with the effect of temperature reported by Gao and coworkers77 in their study of 2 and 5 keV cascades with a hybrid MD model that extracted heat from the simu­lation cell. On the other hand, the difference between the 100 and 600 K results is not large, so the effect of ^200 K of cascade-induced heating may be modest.

A simple assessment of this cascade-induced heat­ing was carried out using 10 keV cascades at 100 K. Two independent sets ofsimulations were carried out, seven simulations in a cell of 128 k atoms and eight simulations in a cell of 250 k atoms. A 10 keV cascade will raise the average temperature by 604 and 309 K, respectively, for these two cell sizes. The results of these simulations are summarized in Figure 11, where the parameters plotted are the surviving defect fraction (per NRT), the fraction of interstitials in clusters (per NRT), and the fraction of interstitials in clusters (per surviving MD defect). In each case, the range ofvalues for the two populations are shown, along with their respective mean values with the standard error indicated. The mean and standard error for the combined data sets is also shown. Although the heating differed by a factor of two, it is clear that the defect survival fraction is essentially identical for both populations. There is a slight trend in the interstitial clustering results, which indicates that a higher temperature (due to a smaller number of atoms) promotes interstitial clustering. This is consis­tent with the results that will be discussed below.

Nucleation of small clusters in supersaturated solu­tions has been of significant interest to several genera­tions ofscientists. The kinetic model for cluster growth and the rate of formation of stable droplets in vapor and second phase precipitation in alloys during aging was studied extensively. The similarity to the con­densation process in supersaturated solutions allows the results obtained to be used in RDT to describe the formation of defect clusters under irradiation.

The initial motivation for work in this area was to derive the nucleation rate of liquid drops. Farkas63 was first to develop a quantitative theory for the so-called homogeneous cluster nucleation. Then, a great number of publications were devoted to the kinetic nucleation theory, of which the works by Becker and Doring,64 Zeldovich,65 and Frenkel66 are most important. Although these publications by no means improved the result of Farkas, their treat­ment is mathematically more elegant and provided a proper background for subsequent works in for­mulating ME and revealing properties of the clus­ter evolution. A quite comprehensive description of the nucleation phenomenon was published by Goodrich.67,68 Detailed discussions of cluster nucle — ation can also be found in several comprehensive reviews.69,70 Generalizations of homogeneous cluster nucleation for the case of irradiation were developed by Katz and Wiedersich71 and Russell.72 Here we only give a short introduction to the theory.

For small cluster sizes at high enough tempera­ture, when the thermal stability of clusters is rela­tively low, the diffusion of clusters in the size-space governs the cluster evolution, which is nucleation of stable clusters. In cases where only FPs are produced by irradiation, the first term on the right-hand side of eqn [18] is equal to zero and cluster nucleation, for example, voids, proceeds via interaction between mobile vacancies to form divacancies, then between vacancies and divacancies to form trivacancies, and so on. By summing eqn [18] from x = 2 to 1, one finds

^ = J (x)k=i = Jnucl [24]

1

where Nc = ^2 f (x) is the total number of clusters.

x=2

The nucleation rate in this case, Jcnucl, is equal to the rate of production of the smallest cluster (divacancies in the case considered); hence the flux J(x)|x=1 is the main concern.

When calculating Jcnucl, one can obtain two limit­ing SDFs that correspond to two different steady — state solutions of eqn [18]: (1) when the flux J(x, t) = 0, for which the corresponding SDF is »(x), and, (2) when it is a constant: J(x, t) = Jc, with the SDF denoted as g(x). Let us first find »(x). Using equation P(x)n(x) — Q(x + 1)»(x + 1, t) = 0 and the condition »(1) = C, one finds that

P(y)

Q (y + 1)

Using function »(x), the flux J(x, t) can be derived as follows

J(x, t) = P(x)n(x)

The SDF g(x) corresponding to the constant flux, J(x, t) = Jc, can be found from eqn [26]:

1

PG0«M

Using the boundary conditions g(1) = »(1) = C one finds that Jcnucl is fully defined by »(x):

nucl _ ______ 1_______

Jc 1

[P(x)n(x)]—1

x=1

Generally, »(x) has a pronounced minimum at some critical size, x = xcr, and the main contribution to the denominator of eqn [28] comes from the clusters with size around xcr. Expanding »(x) in the vicinity of xcr up to the second derivative and replacing the sum­mation by the integration, one finds an equation for Jcnucl, which is equivalent to that for nucleation of second phase precipitate particles.64,65 Note that eqn [28] describes the cluster nucleation rate quite accu­rately even in cases where the nucleation stage coexists with the growth which leads to a decrease of the concentration of mobile defects, C. This can be seen from Figure 2, in which the results of numerical integration of ME for void nucleation are compared with that given by eqn [28].73

In the case of low temperature irradiation, when all vacancy clusters are thermally stable (C = Cv in the case) and only FPs are produced by irradiation,

the void nucleation rate, eqn [21], can be calculated analytically. Indeed, in the case where the binding energy of a vacancy with voids of all sizes is infinite, Ev(x) = 1 (see eqn [75]), it follows from eqn [25] that the function n(x) is equal to

C1 DyCy x 1

x1/3 DiCi

Substituting eqn [29] in eqn [28], one can easily find that the nucleation rate, Jcnucl, takes the form

rnucl

where w = (48я2/02)1/3 is a geometrical factor of the order of 1020m-2 (see Section 1.13.5). The sum in the dominant eqn [30] is a simple geometrical progression and therefore it is equal to

DyCy

DvCv — Di Ci

Substituting eqn [31] to eqn [30], one can finally obtain the following equation

Jcnucl = wCv(DvCv — DiCi) [32]

Note that the function g(x) in this case takes a very simple form, g(x) = Cc/x1/3, and hence decreases with increasing cluster size. In contrast, in R-space,
g(R) (see eqn [16]) increases with increasing cluster size: g(R) = (36я/0)1/3 CcR (see also eqns [43] and [44] in Feder et a/.69).

The real time-dependent SDF builds up around the function g(x) with the steadily increasing size range (see, e. g., Figure 2 in Feder et a/.69). Also note that homogeneous nucleation is the only case where an analytic equation for the nucleation rate exists. In more realistic scenarios, the nucleation is affected by the presence ofimpurities and other crystal imper­fections, and numerical calculations are the only means of investigation. Such calculations are not trivial because for practical purposes it is necessary to consider clusters containing very large numbers of defects and, hence, a large number of equations. This can make the direct numerical solution of ME impractical. As a result several methods have been developed to obtain an approximate numerical solu­tion of ME (see Section 1.13.4.4 for details).

The equations formulated in this section govern the evolution of mobile and immobile defects in solids under irradiation or aging and provide a frame­work, which has been used for about 50 years. Appli­cation of this framework to the models developed to date is presented in Sections 1.13.5 and 1.13.6.

The second example demonstrates mesoscale KMC model simulations of the diffusion of silver (Ag) through the pyrolytic carbon and silicon carbide containment layers of a TRISO nuclear fuel particle. The model atomically resolves Ag, but provides a nonatomic, mesoscale medium of carbon and silicon carbide that includes a variety of defect features including grain boundaries, the carbon-silicon car­bide interfaces, cracks, precipitates, and nanocavities. These defect features can serve as either fast diffu — sional pathways or traps for the migrating silver. The model consists of a 2D slab geometry incorporating the pyrolytic carbon and silicon carbide layers, with incident silver atoms placed at the innermost pyro­lytic carbon layer, as described in more detail in Meric de Bellefon and Wirth.88

The key input parameters to the model (diffusion coefficients, trap binding energies, interface charac­teristics) are determined from available experimental data, or parametrically varied, until more precise values become available from lower length scale modeling. The predicted results, in terms of the time/temperature dependence of silver release dur­ing postirradiation annealing and the variability of silver release from particle to particle have been compared to available experimental data from the German High-Temperature Reactor (HTR) Fuel Program,89 as shown below in Figure 7, and studies performed by the Japan Atomic Energy Research Institute (JAERI).90

Figure 6 presents KMC simulation results, which shows the effect of different grain geometries in SiC on silver release during postirradiation annealing. In this model, the grains are considered to have a rect­angular geometry. The smaller dimension is parallel to the interfaces and has a fixed length of 1 pm. The longer dimension, parallel to the radial direction in a TRISO fuel particle, has a variable length that is uniformly distributed among grains over a range from 1 to 40 pm, as shown in the upper plot of Figure 6. Such a grain distribution mimics a highly columnar structure, as is often observed experimentally.91

The diffusion coefficient for silver transport within the grain boundaries has been assumed to be three orders of magnitude higher than in bulk SiC. As expected, within this model, the presence of grains provides fast diffusion paths for silver trans­port and accelerates the released fraction. Adding a

1

Decreasing grain width range

0.1

OPyC

0.001

columnar structure in those grains further increases the release of silver; the released fraction increases from 50 to 80% when the grain distribution shifts from an isotropic structure (grains are 1 mm squares) to a highly columnar structure (length of grains uni­formly distributed from 1 to 40 mm) after 270 h of heating at 1700 °C.

Figure 7 presents KMC simulation results of the transport of silver through a given PyC/SiC/PyC microstructure during postirradiation thermal an­nealing at 1600, 1700, and 1800 °C, as well as results from three German experimental release measure­ments performed at annealing temperatures of 1600, 1700, and 1800 °C. The simulated microstructures include reflective interfaces, trapping cavities, and a grain boundary structure in the SiC layer. The microstructures for the 1600 and 1700 °C simulations (particle 23 and 24) are identical and contain an isotropic grain geometry in SiC that consists of 1-mm long square grains with a grain boundary diffu — sivity 100x higher than in bulk. In the 1800 °C simu­lation (particle 25), the SiC grain characteristics are varied to match the measured release at 1800 °C. Faster transport through grain boundaries is required to match the experimental results, which is obtained by implementing a highly columnar structure in which the grains are 0.5 mm-wide and have a radial length between 10 and 40 mm into SiC, as well as a much higher grain boundary diffusivity that is 2000 times higher than in bulk.

In contrast to using associates for liquid solutions is a sublattice approach in which cations and anions are mixed on respective lattice sites. With anions and cations assigned to specific sublattices, it is possible
to capture interactions and short-range order with species occupying the sites and additional energetic terms. The components can essentially be allowed to independently mix on each sublattice within the energetic constraints and the system free energy minimized.46 The approach has been successfully used by Gueneau et a/.47 to model the liquid in the U-O and Pu-O systems where ionic metal species reside on one lattice and oxygen anions, neutral UO2 or PuO2, charged vacancies, and O species on the other.

An improvement to the simple sublattice approach is the quasichemical method introduced by Fowler and Guggenheim48 and later further devel­oped by Pelton and coworkers.49-52 It approaches short-range order in liquids through the formation of nearest-neighbor pairs on a quasilattice. It thus differs significantly from the modified associate species approach such that in the quasichemical method short-range order is accommodated by com­ponents pairing and the energetically described extent of like and unlike components pairing. The technique thus avoids the paradox where a high

degree of short-range order in the associate species approach causes minimal end-member species con­tent and therefore fails in the limit to be a Raoult’s law solution.52

In the modified quasichemical approach for a sim­ple binary A—B system, the components are treated as distributed on a quasilattice and that an energetic term governs exchange among the pairs.

(A — A) + (B — B)=2(A — B) [19]

A parameter, Z, represents the nearest-neighbor coordination number such that each component forms Z pairs. For one mole of solution,

ZXA = 2nAA + 2nAB [20]

ZXb = 2»bb + 2»ab [21]

where the moles of each pair are nAA, nBB, and nAB. The relative proportions of each pair is Xj where

Xij = n, j / (»AA + «BB + »ab) [22]

The configurational entropy contribution is captured from the random distribution of the pairs over the quasilattice positions. The result is a heat of mixing of

Hmix = (Xab/2)Eqm [23]

where Eqm is the quasichemical model energetic parameter, and a configurational entropy

Sconfig = — R(Xa ln Xa + Xb ln Xb)

— R(Z/2)(XAAln(XAA/xA ) + Хвв1п(Хвв/Хв2)

+ Хав1п(Хав/2ХаХв)) [24]

Utilizing Gibbs free energy functions for the compo­nents and expanding EQm as a polynomial in a scheme for minimizing the system free energy pro­vide a system for optimization of the liquid using known thermochemical values and phase equilibria. Issues such as the displacement of the composition of maximum short-range order from 50% composition are dealt with by assuming different coordination numbers for each component. Greater accuracy is obtained by inclusion of the Bragg-Williams model, thus incorporating lattice interactions beyond nearest neighbors. The modification to the quasi­chemical model yields

Hmix = (Eqm/2)Xab + EbwXaXb [25]

The extension of the modified quasichemical model to ternary systems is directly possible using only binary model parameters.

An issue for the modified quasichemical model is that it fails at high deviations from random ordering, although that is generally not a problem because immiscibility will occur before the deviations grow too large. The model can also predict a large amount
of ordering that can result in a negative configu­rational entropy, a physical impossibility.53

The unique feature ofelectron irradiations in compari­son to ions and neutrons is that they create defects in very low-energy recoil events. As a consequence, nearly all FPs are produced in isolation. This has been of foremost importance in developing our understanding of radiation damage, as it made studies of defect crea­tion mechanisms as well as the fundamental properties of FPs possible. Recall that the properties of vacancies and vacancy clusters, for example, formation and migra­tion energies, stacking fault energies, etc., could be determined from quenching studies. It is not possible, however, to quench in interstitials in metals. Very little was therefore known about this intrinsic defect prior to about 1955 when irradiation experiments became widely employed. In this section, we highlight some of the key findings derived from these past studies.