As explained above, pseudoization schemes are espe­cially relevant for plane wave codes. All pseudoiza — tion schemes are obtained by calculations on isolated atoms or ions. The real potential experienced by the valence electrons is replaced by a pseudopotential coming from mathematical manipulations. A good pseudopotential should have two apparently contra­dictory qualities. First, it should be soft, meaning that the wave function oscillations should be smoothened as much as possible. For a plane wave basis set, this means that the number of plane waves needed to represent the wave functions is kept minimal. Second, it should be transferable, which means that it should correctly represent the real interactions of valence electrons with the core in any kind of chemical envi­ronment, that is, in any kind of bonding (metallic, covalent, ionic), with all possible ionic charges or covalent configurations conceivable for the element under consideration. The generation of pseudopo­tentials is a rather complicated task, but nowadays libraries of pseudopotentials exist and pseudopoten­tials are freely available for almost any element, though not with all the pseudoization schemes.

One can basically distinguish norm-conserving pseudopotentials, ultrasoft pseudopotentials, and PAW formalism. Norm-conserving pseudopotentials were the first ones designed for ab initio calculations.12 They involve the replacement of the real valence wave function by a smooth wave function of equal norm, hence their name. Such pseudopotentials are rather easy to generate, and several libraries exist with all elements of the periodic table. They are rea­sonably accurate although they are still rather hard, and so they are less and less used in plane wave codes but are still used with atomic-like basis sets. Ultrasoft pseudopotentials13 remove the constraint of norm equality between the real and pseudowave functions. They are thus much softer though less easy to gen­erate than norm-conserving ones. The Projector Augmented Wave14 formalism is a complex pseudoi­zation scheme close in spirit to the ultrasoft scheme but it allows the reconstruction of the real electronic den­sity and the real wave functions with all their oscilla­tions, and for this reason this method can be considered an all-electron method. When correctly generated, PAW atomic data are very soft and quite transferable. Libraries of ultrasoft pseudopotentials or PAW atomic data exist, but they are generally either incomplete or not freely available.

Plane wave codes in use in the nuclear materials community include VASP15 with ultrasoft pseudopo­tentials and PAW formalism, Quantum-Espresso16 with norm-conserving and ultrasoft pseudopotentials and PAW formalism, and ABINIT17 with norm — conserving pseudopotentials and PAW formalism.

Note that for a specific pseudoization scheme many different pseudopotentials can exist for a given element. Even if they were built using the same valence orbitals, pseudopotentials can differ by many numerical choices (e. g., the various matching radii) that enter the pseudoization process.

We present in the following a series of practical choices to be made when one wants to perform ab initio calculations. But the first and certainly most important of these choices is that of the ab initio code itself as different codes have different speeds, accuracies, numerical methods, features, input files, and so on, and so it proves quite difficult to change codes in the middle of a study. Furthermore, one observes that most people are reluctant to change their usual code as the investment required to fully master the use of a code is far from negligible (not to mention the one to master what is in the code).

In this test case, let us consider a body-centered cubic (bcc) crystal of Tantalum (Ta), described by the Finnis-Sinclair (FS) potential.31 The calculations are performed using the MD++ program. The source code and the input files for this and subsequent test cases in this chapter can be downloaded from http://micro. stanford. edu/wiki/Comprehensive_ Nuclear_Materials_MD_Case_Studies.

The cut-off radius of the FS potential for Ta is 4.20 A. To avoid interaction between an atom with its own periodic images, we consider a cubic simulation cell whose size is much larger than the cut-off radius. The cell dimensions are 5[100], 5[010], and 5[001] along x, y, and z directions, and the cell contains N = 250 atoms (because each unit cell of a bcc crystal contains two atoms). PBC are applied in all three directions. The experimental value of the equi­librium lattice constant ofTa is 3.3058 A. Therefore, to compute the equilibrium lattice constant of this potential model, we vary the lattice constant a from 3.296 to 3.316 iA, in steps of 0.001 JA. The potential energy per atom E as a function of a is plotted in Figure 5. The data can be fitted to a parabola. The

ao (A)

location of the minimum is the equilibrium lattice constant, a0 = 3.3058 /A. This exactly matches the experimental data because a0 is one of the fitted parameters of the potential. The energy per atom at a0 is the cohesive energy, Ecoh = —8.100 eV, which is another fitted parameter. The curvature of parabolic curve at a0 gives an estimate of the bulk modulus, B = 197.2 GPa. However, this is not a very accurate estimate of the bulk modulus because the range of a is still too large. For a more accurate determination of the bulk modulus, we need to compute the E(a) curve again in the range of |a — a01 < 10 4 A. The curvature of the E(a) curve at a0 evaluated in the second calcu­lation gives B = 196.1 GPa, which is the fitted bulk modulus value of this potential model.31

When the crystal has several competing phases (such as bcc, face-centered cubic, and hexagonal — closed-packed), plotting the energy versus volume (per atom) curves for all the phases on the same graph allows us to determine the most stable phase at zero temperature and zero pressure. It also allows us to predict whether the crystal will undergo a phase transition under pressure.32

Other elastic constants besides B can be computed using similar approaches, that is, by imposing a strain on the crystal and monitoring the changes in potential energy. In practice, it is more convenient to extract the elastic constant information from the stress-strain relationship. For cubic crystals, such as Ta considered here, there are only three independent elastic constants, Cn, C12, and C44. C11 and C12 can be obtained by elongating the simulation cell in the x-direction, that is, by changing the cell length into L =(1 + exx) • L0, where L0 = 5a0 in this test case. This leads to nonzero stress components sxx, ayy szz, as computed from the Virial stress formula [19], as shown in Figure 6 (the atomic velocities are zero because this calculation is quasistatic). The slope of these curves gives two of the elastic con­stants C11 = 266.0 GPa and C12 = 161.2 GPa. These results can be compared with the bulk modulus obtained from potential energy, due to the relation B = (Cn + 2 C12) / 3 = 196.1GPa.

C44 can be obtained by computing the shear stress axy caused by a shear strain £xy Shear strain exy can be applied by adding an off-diagonal element in matrix H that relates scaled and real coordinates of atoms.

The slope of the shear stress-strain curve gives the elastic constant C44 = 82.4 GPa.

In this test case, all atoms are displaced according to a uniform strain, that is, the scaled coordinates of all atoms remain unchanged. This is correct for simple crystal structures where the basis contains only one atom. For complex crystal structures with more than one basis atom (such as the diamond-cubic structure of silicon), the relative positions of atoms in the basis set will undergo additional adjustments when the crystal is subjected to a macroscopically uniform strain. This effect can be captured by performing energy minimization at each value of the strain before recording the potential energy or the Virial stress values. The resulting ‘relaxed’ elastic constants corre­spond well with the experimentally measured values, whereas the ‘unrelaxed’ elastic constants usually overestimate the experimental values.

Given the short time scale and small volume associated with atomic displacement cascades, it is not currently possible to directly observe their behavior by any avail­able experimental method. Some of their characteristics have been inferred by experimental techniques that can examine the fine microstructural features that form after low doses of irradiation. The experimental work that provides the best estimate of stable Frenkel pair produc­tion involves cryogenic irradiation and subsequent annealing while measuring a parameter such as electri­cal resistivity.26,27 Less direct experimental measure­ments include small angle neutron scattering,28 X-ray scattering, positron annihilation spectroscopy, 0 and field ion microscopy. 1 More broadly, transmission elec­tron microscopy (TEM) has been used to characterize the small point defect clusters such as microvoids, dislocation loops, and stacking fault tetrahedra that are formed as the cascade collapses.32-36

The primary tool for investigating radiation dam­age formation in displacement cascades has been computer simulation using MD, which is a computa­tionally intensive method for modeling atomic sys­tems on the time and length scales appropriate to displacement cascades. The method was pioneered by Vineyard and coworkers at Brookhaven National Laboratory,37 and much of the early work on atomistic simulations is collected in a review by Beeler.38 Other methods, such as those based on the BCA,2 , 1 have also been used to study displacement cascades. The binary collision models are well suited for very high — energy events, which require that the interatomic potential accurately simulate only close encounters between pairs of atoms. This method requires sub­stantially less computer time than MD but provides less detailed information about lower energy colli­sions where many-body effects become important. In addition, in-cascade recombination and clustering can only be treated parametrically in the BCA. When the necessary parameters have been calibrated using the results of an appropriate database of MD cascade results, the BCA codes have been shown to reproduce the results of MD simulations reasonably well.39,40

A detailed description of the MD method is given in Chapter 1.09, Molecular Dynamics, and

will not be repeated here. Briefly, the method relies on obtaining a sufficiently accurate analytical inter­atomic potential function that describes the energy of the atomic system and the forces on each atom as a function of its position relative to the other atoms in the system. This function must account for both attractive and repulsive forces to obtain the appropri­ate stable lattice configuration. Specific values for the adjustable coefficients in the function are obtained by ensuring that the interatomic potential leads to reasonable agreement with measured material para­meters such as the lattice parameter, lattice cohesive energy, single crystal elastic constants, melting tem­perature, and point defect formation energies. The process of developing and fitting interatomic poten­tials is the subject of Chapter 1.10, Interatomic Potential Development. One unique aspect arises when using MD and an empirical potential to investi­gate radiation damage, viz. the distance of closest approach for highly energetic atoms is much smaller than that obtained in any equilibrium condition. Most potentials are developed to describe equilibrium con­ditions and must be modified or ‘stiffened’ to account for these short-range interactions. Chapter 1.10, Interatomic Potential Development, discusses a common approach in which a screened Coulomb potential is joined to the equilibrium potential for this purpose. However, as Malerba points out,41 criti­cal aspects of cascade behavior can be sensitive to the details of this joining process.

When this interatomic potential has been derived, the total energy of the system of atoms being simulated can be calculated by summing over all the atoms. The forces on the atoms are obtained from the gradient of the interatomic potential. These forces can be used to calculate the atom’s accelerations according to Newton’s second law, the familiar F = ma (force = mass x acceleration), and the equations of motion for the atoms can be solved by numerical integration using a suitably small time step. At the end of the time step, the forces are recalculated for the new atomic positions and this process is repeated as long as necessary to reach the time or state of interest. For energetic PKA, the initial time step may range from ^1 to 10 x 10~18 s, with the maximum time step limited to ~1—10 x 10~15s to maintain acceptable numerical accuracy in the integration. As a result, MD cascade simulations are typically not run for times longer than 10—100 ps. With periodic boundary conditions, the size of the simulation cell needs to be
large enough to prevent the cascade from interacting with periodic images of itself. Higher energy events therefore require a larger number of atoms in the cell. Typical MD cascade energies and the approxi­mate number of atoms required in the simulation are listed in Table 1. With periodic boundaries, it is important that the cell size be large enough to avoid cascade self-interaction. For a given energy, this size depends on the material and, for a given material, on the interatomic potential used. Different inter­atomic potentials may predict significantly different cascade volumes, even though little variation is even­tually found in the number of stable Frenkel pair.42 Using a modest number of processors on a modern parallel computer, the clock time required to com­plete a high-energy simulation with several million atoms is generally less than 48 h. Longer-term evolu­tion of the cascade-produced defect structure can be carried out using Monte Carlo (MC) methods as discussed in Chapter 1.14, Kinetic Monte Carlo Simulations of Irradiation Effects.

The process of conducting a cascade simulation requires two steps. First, a block of atoms of the desired size is thermally equilibrated. This permits the lattice thermal vibrations (phonon waves) to be established for the simulated temperature and typi­cally requires a simulation time of approximately 10 ps. This equilibrated atom block can be saved and used as the starting point for several subsequent cas­cade simulations. Subsequently, the cascade simula­tions are initiated by giving one ofthe atoms a defined amount of kinetic energy, EMD, in a specified direc­tion. Statistical variability can be introduced by either
further equilibration of the starting block, choosing a different PKA or PKA direction, or some combination of these. The number of simulations required at any one condition to obtain a good statistical description of defect production is not large. Typically, only about 8-10 simulations are required to obtain a small stan­dard error about the mean number of defects pro­duced; the scatter in defect clustering parameters is larger. This topic will be discussed further below when the results are presented. Most of the cascade simulations discussed below were generated using a [135] PKA direction to minimize directional effects such as channeling and directions with particularly low or high displacement thresholds. The objective has been to determine mean behavior, and investigations of the effect of PKA direction generally indicate that mean values obtained from [135] cascades are repre­sentative of the average defect production expected in cascades greater than about 1 keV.43 A stronger influ­ence of PKA direction can be observed at lower ener­gies as discussed in Stoller and coworkers.4 ,

In the course of the simulation, some procedure must be applied to determine which of the atoms should be characterized as being in a defect state for the purpose of visualization and analysis. One approach is to search the volume of a Wigner-Seitz cell, which is centered on one of the original, perfect lattice sites. An empty cell indicates the presence of a vacancy and a cell containing more than one atom indicates an interstitial-type defect. A more simple geometric criterion has been used to identify defects in most of the results presented below. A sphere with a radius equal to 30% of the iron lattice parameter is

centered on the perfect lattice sites, and a search similar to that just described for the Wigner-Seitz cell is carried out. Any atom that is not within such a sphere is identified as part of an interstitial defect and each empty sphere identifies the location of a vacancy. The diameter of the effective sphere is slightly less than the spacing of the two atoms in a dumbbell interstitial (see below). A comparison of the effective sphere and Wigner-Seitz cell approaches found no significant difference in the number of stable point defects identified at the end of cascade simulation, and the effective sphere method is faster computationally. The drawback to this approach is that the number of defects identi­fied by the algorithm must be corrected to account for the nature of the interstitial defect that is formed. In order to minimize the lattice strain energy, most interstitials are found in the dumbbell configuration; the energy is reduced by distributing the distortion over multiple lattice sites. In this case, the single interstitial appears to be composed of two interstitials separated by a vacancy. In other cases, the interstitial configuration is extended further, as in the case of the crowdion in which an interstitial may be visualized as three displaced atoms and two empty lattice sites. These interstitial configurations are illustrated in Figure 4, which uses the convention adopted throughout this chapter, that is, vacancies are displayed as red spheres and interstitials as green spheres. A simple postprocessing code was used to determine the true number of point defects, which are reported below.

Most MD codes describe only the elastic colli­sions between atoms; they do not account for energy

loss mechanisms such as electronic excitation and ionization. Thus, the initial kinetic energy, EMD, given to the simulated PKA in MD simulations is more analogous to Td in eqn [2] than it is to the PKA energy, which is the total kinetic energy of the recoil in an actual collision. Using the values of EMD in Table 1 as a basis, the corresponding EPKA and nNRT for iron, and the ratio of the damage energy to the PKA energy, have been calculated using the procedure described in Norgett and coworkers.19 and the recommended 40 eV displacement thresh­old.16 These values are also listed in Table 1, along with the neutron energy that would yield EPKA as the average recoil energy in iron. This is one-half of the maximum energy given by eqn [1]. As mentioned above, the difference between the MD cascade energy, or damage energy, and the PKA energy increases as the PKA energy increases. Discussions of cascade energy in the literature on MD cascade simulations are not consistent with respect to the use of the term PKA energy. The third curve in Figure 3 shows the calculated number of Frenkel pair predicted by the NRT model if the PKA energy is used in eqn [2] rather than the damage energy. The difference between the two sets of NRT values is substantial and is a measure of the ambiguity associated with being vague in the use of terminology. It is recom­mended that the MD cascade energy should not be referred to as the PKA energy. For the purpose of comparing MD results to the NRT model, the MD cascade energy should be considered as approxi­mately equal to the damage energy (Td in eqn [2]).

In reality, energetic atoms lose energy continu­ously by a combination of electronic and nuclear reactions, and the typical MD simulation effectively deletes the electronic component at time zero. The effects of continuous energy loss on defect produc­tion have been investigated in the past using a damp­ing term to slowly remove kinetic energy.46 The related issues of how this extracted energy heats the electron system and the effects of electron-phonon coupling on local temperature have also been exam — ined.47-50 More recently, computational and algorith­mic advances have enabled these phenomena to be investigated with higher fidelity.51 Some of the work just referenced has shown that accounting for the electronic system has a modest quantitative effect on defect formation in displacement cascades. For exam­ple, Gao and coworkers found a systematic increase in defect formation as they increased the effective electron-phonon coupling in 2, 5, and 10 keV cascade simulations in iron,50 and a similar effect was reported
by Finnis and coworkers.47 However, the primary physical mechanisms of defect formation that are the focus of this chapter can be understood in the absence of these effects.

Single vacancies and other vacancy-type defects, such as, SFTs and dislocation loops, have been con­sidered quite extensively since the 1930s because it was recognized that they define many properties of solids under equilibrium conditions. Extensive infor­mation on defect properties was collected before material behavior in irradiation environments became a problem of practical importance. Qualitatively new crystal defects, SIAs and SIA clusters, were required to describe the phenomena in solids under irradiation conditions. This has been studied compre­hensively during the last ^40 years. The properties of these defects and their interaction with other defects are quite different compared to those of the vacancy — type. Correspondingly, the crystal behavior under irradiation is also qualitatively different from that under equilibrium conditions. The basic properties of vacancy — and SIA-type defects are summarized below.

1.13.3.2.1 Point defects

The basic properties of PDs are as follows:

1. Both vacancies and SIAs are highly mobile at tem­peratures of practical interest, and the diffusion coefficient of SIAs, Di, is much higher than that of vacancies, Dv : Д ^ Dv.

2. The relaxation volume of an SIA is much larger than that of a vacancy, resulting in higher interaction energy with edge dislocations and other defects.

3. Vacancies and SIAs are defects of opposite type, and their interaction leads to mutual recombination.

4. SIAs, in contrast to vacancies, may exist in several different configurations providing differ­ent mechanisms of their migration.

5. PDs of both types are eliminated at fixed sinks, such as voids and dislocations.

The first property leads to a specific temperature dependence of the damage accumulation: only limited number of defects can be accumulated at irradiation temperature below the recovery stage III, when vacan­cies are immobile. At higher temperature, when both PDs are mobile, the defect accumulation is practically unlimited. The second property is the origin of the so-called ‘dislocation bias’ (see Section 1.13.5.2) and, as proposed by Greenwood et a/.,47 is the reason for void swelling. A similar mechanism, but induced by external stress, was proposed in the so-called ‘SIPA’ (stress-induced preferential absorption) model of irradiation creep.48-53 The third property provides a decrease of the number of defects accumulated in a crystal under irradiation. The last property, which is quite different compared to that of vacancies leads to a variety of specific phenomena and will be consid­ered in the following sections.

The ‘irradiation’ rate, that is, the rate of impinging particles in the case of neutron and ion irradiation, is usually transformed into a production rate (number per unit time and volume) of randomly distributed displacement cascades of different energies (5, 10, 20, … keV) as well as residual Frenkel pairs (FPs). New cascade debris are then injected randomly into the simulation box at the corresponding rate. The cascade debris can be obtained by MD simulations for different recoil energies T, or introduced on the basis of the number of FP expected from displace­ment damage theory. In the case of KMC simulation of electron irradiation, FPs are introduced randomly in the simulation box according to a certain dose rate, assuming most of the time that each electron is responsible for the formation of only one FP. This assumption is valid for electrons with energies close to 1 MeV (much lower energy electrons may not produce any FP, whereas higher energy ones may produce small displacement cascades with the forma­tion of several vacancies and SIAs).

The dose is updated by adding the incremental dose associated with the scattering event of recoil energy T, using the Norgett-Robinson-Torrens expression8 for the number of displaced atoms. In this model, the accumulated displacement per atom (dpa) is given by:

0 8 T

Displacement per subcascade = — [3]

where T is the damage energy, that is, the fraction of the energy of the particle transmitted to the PKA as kinetic energy and ED is the displacement threshold energy (e. g., 40 eV for Fe and reactor pressure vessel (RPV) steels51).

1.14.4.1 Transmutation Rate

The rate of producing transmutations can also be included in KMC models, as deduced from the reac­tion rate density determined from the product of the neutron cross-section and neutron flux. Like the irra­diation rate, the volumetric production rate is used to introduce an appropriate number of transmutants, such as helium that is produced by (n, a) reactions in the fusion neutron environment, where the species are introduced at random locations within the material.

1.14.4.2 Diffusion Rate

Usually the rates of diffusion can be obtained from the knowledge of the migration barriers which have to be known for all the diffusing ‘objects’; that is, for the point defects in AKMC, OKMC, and EKMC or the clusters in OKMC or EKMC. For isolated point

defects, the migration barriers can be from experi­mental data, that is, from diffusion coefficients, or theoretically, using either ab initio calculations as described in Caturla et a/.49 and Becquart and Domain50 or MD simulations as described in Soneda and Diaz de la Rubia.22 Since the migration energy depends on the local environment of the jumping species, it is generally not possible to calculate all of the possible activation barriers using ab initio or even MD simulations. Simpler schemes such as broken bond models, as described in Soisson et a/.,52 Le Bouar and Soisson,53 and Schmauder and Binkele,54 are then used. Another kind of simpler model is based on the calculation of the system configurational energies before and after the defect jump. In this model, the activation energy is obtained from the final Ef and the initial E as follows:

DE

Ea0 + ~2~ where Ea0 is the energy of the moving species at the saddle point. The modification of the jump activation energy by DE represents an attempt to model the effect of the local environment on the jump frequen­cies. Indeed, detailed molecular statics calculations suggest that this represents an upper-bound influence of the effect,55 and although this is a very simplified model, the advantage is that this assumption maintains the detailed balance of jumps to neighboring positions.

The system configurational energies E and Ef, as well as the energy of the moving species at the saddle point Ea0 can be determined using interatomic poten­tials as described in Becquart et a/.,26 Bonny et a/.,44 Wirth and Odette,55 and Djurabekova et a/.56 when they exist. However, at present, this situation is only available for simple binary or ternary alloys. This approach allows one to implicitly take into account relaxation effects as the energy at the saddle point which is used in the KMC and is obtained after relax­ation of all the atoms. The challenge in that case is the total number of barriers to be calculated, which is determined by the number of nearest neighbor sites included in the definition ofthe local atomic environ­ment. Without considering symmetries, this number is sN, where s is the number of species in the system. In spite of using the fast techniques that were devel­oped to find saddle points on the f/y such as the dimer method,57 the nudged elastic band (NEB) method,58 or eigen-vector following methods,59 this number quickly becomes unmanageable. Ideally, the alterna­tive should be to find patterns in the dependence of the energy barriers on the configuration. This is the
approach chosen by Djurabekova and coworkers,56 using artificial intelligence systems. For more complex alloys, for which no interatomic potentials exist, Ei and Ef can be estimated using neighbor pair interactions60- 63 A recent example of the fitting procedure of a neigh­bor pair interactions model can be found in Ngayam Happy et a/63 A discussion of the two approaches applied to the Fe-Cu system has been published by Vincent et a/64 Also note that in the last 10 years, methods in which the possible transitions are found in some systematic way from the atomic forces rather than by simply assuming the transition mechanism a priori (e. g., activation-relaxation technique (ART) or dimer methods)65-68 have been devised. The accuracy of the simulations is thus improved as fewer assumptions are made within the model. However, interatomic poten­tials or a corresponding method to obtain the forces acting between atoms for all possible configurations is necessary and this limits the range of materials that can be modeled with these clever schemes.

The attempt frequency (nX in eqn [1]) can be calculated on the basis of the Vineyard theory69 or can be adjusted so as to reproduce model experiments.

Regular or subregular solution and variable stoichi­ometry representations, while relatively successful, lack a sense of reproducing physical processes. Specifically, they are constrained with regard to accu­rately dealing with entropy contributions because of the defect structures in nonstoichiometric phases and substitutional solutions. A practical advance has
been the sublattice approach, which has been further refined for crystalline systems in the compound energy formalism (CEF). As typical for cation — anion systems, the structure of a phase can be repre­sented by a formula, for example, (A, B)k(D, E)/ where A and B mix on one sublattice and D and E mix on a second sublattice. The constitution of the phase is made up of occupied site fractions, and allowing one of the constituents to be a vacancy permits treatment of nonstoichiometric systems.

Even with a sublattice approach such as CEF, the relationship of eqn [10] is still applicable, but with an interpretation related to a sublattice model. The sum
of the standard Gibbs free energies in this case is the sum of the values for the paired sublattice constitu­ents, which for the example above might be AkD/. Each is a unique set with the Gibbs free energies for the constituents derived from the end-member standard Gibbs free energies, typically through sim­ple geometric additions with any necessary additional configurational entropy contributions. The entropy contribution from mixing on the sublattice sites is defined as

Gid = R^zsys ln(ys) [14]

where z is the stoichiometric coefficient, s defines the lattice, and y is the site fraction for species j. Excess terms represent the interaction energetics between each set of sublattice constituents, for exam­ple, AkD/: BkD/. Again, a Redlich-Kister-Muggianu formulation that includes expansion terms for inter­actions between the constituents can be used:

Gxs = EXE j :*

t j t

+ EEE*jv [15]

t j к

where the sums are associated with components on each sublattice 1 and 2 and the L values are terms for the interaction energies between cations i and j on one sublattice when the other sublattice is occupied only by cation к, and vice versa for the second term.

The PuO2_x phase has been successfully repre­sented by a CEF approach by Gueneau et al. The phase can be described by two sublattices with vacan­cies only on the anion sites

(Pu4+,Pu3+)'(O2~,Va)2.

Including the end members, the constituent species are then

(Pu4+MO2-)2,

(Pu4+)'(Va)2,

(Pu3+MO2-)2, and (Pu3+MVa)2.

A schematic of the relationship between the constituents is seen in Figure 5 where the corners represent each of the constituents listed above. The charged constituents must sum to neutrality, and the line designating neutrality is seen in Figure 5. Gibbs free energy expressions for each of the units can be determined from standard state values. Optimizations using all available thermochemical information, for example, oxygen potentials and phase equilibria, can thus yield the necessary corrections to the Gibbs free energies for the nonstandard constituents together with obtained interaction parameters (L values). The results are shown in Figure 6 where oxygen potential isotherms overlay the phase diagram and which shows mO2 results of models for other phases in the system.

The CEF approach has recently begun to be more widely applied to nuclear fuels. Besides the PuO2_x system noted above, Gueneau eta/.31 also applied the model to accurately describe solid solution phases in the U-O system, as has Chevalier et a/.32 who also addressed the U-O-Zr ternary system.33 Kinoshita et a/.34 used a sublattice approach to model fluorite — structure oxides including ThO2_x and NpO2_x although they did not include charged ionic cations and anions on the sublattices. Zinkevich et a/.35 suc­cessfully modeled the CeO2_x phase using the CEF approach in their comprehensive assessment of the Ce-O phase diagram.

The first mean-field lattice rate models included two thermally activated jump frequencies, one for the vacancy and the other for the interstitial. A direct interstitial diffusion mechanism14 and later a dumb­bell diffusion mechanism1 have been modeled in detail. The vacancy jump frequency parameters are fitted to available thermodynamic and tracer dif­fusion data, and the interstitial parameters are fitted to effective migration energies derived from resistiv­ity recovery measurements.121 The resulting local — concentration-dependent jump frequencies describe both the kinetics of thermal alloys toward equilib­rium and irradiation-induced surface segregation in concentrated alloys. The surface and its vicinity are modeled by the stacking of N parallel atomic planes perpendicular to the diffusion axis, which is taken to be a [100] direction of an fcc alloy. A mirror bound­ary condition is used at one end, and a free surface, which can act as both a source and sink for point defects at the other end. Above the surface, a buffer plane almost full of vacancies is added. Fluxes between the buffer plane and vacuum are forbidden. The resulting equilibrium segregation profiles are controlled by the nominal composition, temperature, and two interaction contributions, the first one expressed in terms of the surface tensions and the second in terms of the ordering energies. Note, that the predicted equilibrium vacancy concentration at the surface is much higher than in the bulk.

Time dependence of mean occupations in an atomic plane of point defects and atoms results from a balance between averaged incoming and out­going fluxes. Fluxes are written within a mean-field approximation, decoupling the statistical averages into a product of mean occupations and mean jump frequencies for which the occupation numbers in the exponential argument are replaced by the corresponding mean occupations. The resulting first order differential kinetic equations are integrated using a predictor corrector variable time step algo­rithm because of the high jump frequency disparities between vacancies and interstitials.

It is observed that interstitial contribution to RIS is of the same order as that of the vacancy. The predicted formation of a ‘W-shaped’ profile as a tran­sient state from the preirradiated enrichment to the strong depletion of Cr is shown to be governed by both thermodynamic properties and the relative values of the transport coefficients between Fe, Cr, and Ni (Figure 9). Thermodynamics not only plays a part in the transport coefficients but also arises in the establishment of a local equilibrium between the surface and the adjacent plane, explaining the oscil­latory behavior of the Cr profile: an equilibrium tendency toward an enrichment of Cr at the grain boundary plane, which competes with a Cr depletion tendency under irradiation. However, the predicted profile is not as wide as the experimental one.

What needs to be improved is the interstitial dif­fusion model. The lack of experimental and ab initio data leads to approximate interstitial jump frequen­cies. Coupling between fluxes is described partially as correlation effects are not accounted for. The recent mean-field developments98 should be integrated in this type of simulation.

The actual number of defects that survive the dis­placement cascade and their spatial distribution in solids will determine the effect on the irradiated microstructure. Figure 7 summarizes the effect of

damage morphology from the viewpoint of the grain boundary and how the defect flow affects radiation — induced grain boundary segregation. Of the total defects produced by the energetic particle, a fraction appears as isolated, or freely migrating defects, and the balance is part of the cascade. The fraction of the ‘ballistically’ produced FPs that survive the cascade quench and are available for long-range migration is an extremely important quantity and is called the migration efficiency, e. These ‘freely migrating’ or ‘avail­able migrating’ defects10 are the only defects that will affect the amount of grain boundary segregation, which is one measure of radiation effects. The migra­tion efficiency can be very small, approaching a few percent at high temperatures. The migration effi­ciency, e, comprises three components:

gi v: the isolated point defect fraction,

8iv: clustered fraction including mobile defect clusters such as di-interstitials, and Z: fraction initially in isolated or clustered form after the cascade quench that is annihilated during subsequent short-term (>10~ns) intracascade thermal diffusion.

They are related as follows:

e = di + gi + Ci = dv + gv + Cv I11]

Figure 8 shows the history of defects born as vacan­cies and interstitials as described by the NRT model.

_______ у______

Clustered point defect fraction (CDF) (div)

Due to significant recombination in the cascade, only a fraction (^30%) is free to migrate from the displacement zone. These defects can recombine out­side of the cascade region, be absorbed at sinks in the
matrix (voids, loops), or be absorbed at the grain boundaries, providing for the possibility of radiation — induced segregation.

The fraction of defects that will be annihilated after the cascade quench by recombination events among defect clusters and point defects within the same cascade (intracascade recombination), Z, is about 0.07, for a migration efficiency of 0.3 (see below for additional detail).10 The clustered fraction, d includes large, sessile clusters and small defect clusters that may be mobile at a given irradiation temperature and will be different for vacancies and interstitials. For a 5 keV cascade, di is about 0.06 and dv is closer to 0.18.10 Some of these defects may be able to ‘evaporate’ or escape the cluster and become ‘available’ defects (Figure 8).

This leaves g, the isolated point defect fraction that are available to migrate to sinks, to form clus­ters, to interact with existing clusters, and to partic­ipate in the defect flow to grain boundaries that gives rise to radiation-induced segregation. Owing to their potential to so strongly influence the irra­diated microstructure, defects in this category, along with defects freed from clusters, make up the freely migrating defect (FMD) fraction. Recall that electrons and light ions produce a large fraction of their defects as isolated FPs, thus increasing the likeli­hood of their remaining as isolated rather than clus­tered defects. Despite the equivalence in energy among the four particle types described in Figure 5, the average energy transferred and the defect pro­duction efficiencies vary by more than an order of magnitude. This is explained by the differences in the cascade morphology among the different parti­cle types. Neutrons and heavy ions produce dense cascades that result in substantial recombination during the cooling or quenching phase. However, electrons are just capable of producing a few widely spaced FPs that have a low probability of recombi­nation. Protons produce small widely spaced cas­cades and many isolated FPs due to the Coulomb interaction and therefore, fall between the extremes in displacement efficiency defined by electrons and neutrons.

The value of g has been estimated to range from 0.01 to 0.10 depending on PKA energy and irradia­tion temperature, with higher temperatures resulting in the lower values. Naundorf12 estimated the freely migrating defect fraction using an analytical treat­ment based on two factors: (1) energy transfer to atoms is only sufficient to create a single FP, and (2) the FP lies outside a recombination (interaction)

radius so that the nearby FPs neither recombine nor cluster. The model follows each generation of the collision and calculates the fraction of all defects produced that remain free. Results of calculation using the Naundorf model are shown in Table 1 for several ions of varying mass and energy. Values of Z range between 24% for proton irradiation to 3% for heavy ion (krypton) irradiation. Recent results,13 however, have shown that the low values of FMD efficiency for heavy ion or neutron irradiation cannot be explained by defect annihilation within the parent cascade (intracascade annihilation). In fact, cascade damage generates vacancy and interstitial clusters that act as annihilation sites for FMD, reducing the efficiency of FMD production. Thus, the cascade remnants result in an increase in the sink strength for point defects and along with recombination in the original cascade, account for the low FMD efficiency measured by experiment.

Ferritic steels are an important class ofnuclear mate­rials, which include reactor pressure vessel (RPV) steels and high chromium steels for elevated temper­ature structural and cladding materials in fast reactors and fusion reactors, see Chapter 4.03, Ferritic Steels and Advanced Ferritic-Martensitic Steels. From a basic science point of view, the modeling of these materials starts with that of pure iron, in the ferro­magnetic bcc structure. Iron presents several difficul­ties for DFT calculations. First, being a three dimensional (3D) metal, it requires rather large basis sets in plane wave calculations. Second, the calculations need to be spin polarized, to account for magnetism, and this at least doubles the calculation time. But most of all, it is a case where the choice of the exchange — correlation functional has a dramatic effect on bulk properties. The standard LDA incorrectly predicts the paramagnetic face-centered cubic (fcc) structure to be more stable than the ferromagnetic bcc structure. The correct ground state is recovered using gradient corrected functionals,39 as illustrated in Figure 1. Finally, it was pointed out that pseudopotentials tend to overestimate the magnetic energy in iron,40 and therefore, some pseudopotentials suffer from a lack of transferability for some properties. In practice, how­ever, in the large set of the results obtained over the last decade for defect calculations in iron, a quite remark­able agreement is obtained between the various computational approaches. With a few exceptions, they are indeed quite independent on the form of the GGA functional, the basis set (plane wave or localized), and the pseudopotential or the use ofPAW approaches.

1.08.4.1.1 Self-interstitials and self­interstitial clusters in Fe and other bcc metals

The structure and migration mechanism of self­interstitials in iron is a very good illustrative example of the impact of DFT calculations on radiation defect

studies. Progress in methods, codes, and computer performance made this archetype of radiation defects accessible to DFT calculations in the early 2000s, since total energy differences between simu­lation cells of 128+1 atoms could then be obtained with a sufficient accuracy. In 2001, Domain and Becquart reported that, in agreement with the experiment, the (110) dumbbell was the most stable structure.41 Quite unexpectedly, the (111) dumbbell was predicted to be ^0.7 eV higher in energy, at variance with empirical potential results that pre­dicted a much smaller energy difference. DFT cal­culations performed in other bcc metals revealed that this is a peculiarity of Fe,42 as illustrated in Figure 2, and magnetism was proposed to be the origin of the energy increase in the (111) dumb­bell in Fe. The important consequence of this result in Fe, which has been confirmed repeatedly since

Figure 2 Formation energies of several basic SIA configurations calculated for bcc transition metals of group 5B (left) and group 6B (right), taken from Nguyen-Manh et a/.42 Data for bcc Fe are taken from Fu et a/.43 Reproduced from Nguyen-Manh, D.; Horsfiels, A. P.; Dudarev, S. L. Phys. Rev. B 2006, 73, 020101.

then, is that it excludes the SIA migration to occur by long 1D glides of the (111) dumbbell followed by on-site rotations of the (110) dumbbell, as predicted previously from empirical potential MD simulations. Moreover, DFT investigation of the migration mechanism yielded a quantitative agreement with the experiment for the energy of the Johnson translation-rotation mechanism (see Figure 3), namely ^0.3 eV.43

These DFT calculations were followed by a very successful example of synergy between DFT and empirical potentials. The DFT values of interstitial formation energies in various configurations and interatomic forces in a liquid model have indeed been included in the database for a fit of EAM type potentials by Mendelev et a/.45 This approach has resulted in a new generation of improved empirical potentials, albeit still with some limitations. When considering SIA clusters made of parallel dumbbells, the Mendelev potential agrees with DFT for predict­ing a crossover as a function of cluster size from the (110) to the (111) orientation between 4 and 6 SIA
clusters.44 However, discrepancies are found when considering nonparallel configurations.46 More pre­cisely, new configurations of small SIA clusters were observed in MD simulations performed at high tem­perature with the Mendelev potential. The energy ofthe new di-interstitial cluster, made ofa triangle of atoms sharing one site (see Figure 4), is even lower than that of the parallel configuration within DFT but higher by 0.3 eV with the Mendelev potential (see also Section 1.08.4.3 on dislocations). The new tri — and quadri-interstitial clusters, with a ring structure (see Figure 4), are one ofthe few examples in which a significant discrepancy is found between various DFT approaches. Calculations with the most accu­rate description of the ionic cores predict that the new tri-interstitial configuration is slightly more sta­ble than the parallel configuration, whereas more approximate ones predict that it is 0.7 eV higher. The first category includes calculations in the PAW approach, performed using either the VASP code or the PWSCF code and also ultrasoft pseudopotential calculations. The second one includes calculations

with less transferable ultrasoft pseudopotentials with VASP and norm-conserving pseudopotentials with SIESTA.46 Such a discrepancy is not common in defect calculations in metals. Further investigations are required to understand more precisely its origin, in particular the possible role of magnetism.

The structures of the most stable SIA clusters in Fe, and more generally of their energy landscape, remain an open question. One would ideally need to combine DFT calculations with methods for ex­ploring the energy surface, such as the Dimer47 or ART4 methods. Such a combination is possible in principle, and it has indeed been used for defects in semiconductors,49 but due to computer limitations this is not the case yet in Fe. The alternative is to develop new empirical potentials in better agreement with

DFT energies in particular for these new structures, to perform the Dimer or ART calculations with these potentials, and to validate the main features of the energy landscape thus obtained by DFT calculations.

To summarize, the energy landscape of interstitial type defects has been revisited in the last decade driven by DFT calculations, in synergy with empiri­cal potential calculations.

For a free electron gas with Fermi wavevector kF, the energy U of volume O is5 h2kF

10я2 me

This contribution to the energy of the condensed phase generates no interatomic force since U is independent of the atomic positions. However, its contribution is significant: metallic cohesive energy and bulk moduli are correct to within an order of magnitude. Consider­ation of this term gives some justification for ignoring the cohesive energy and bulk modulus in fitting a potential, and fitting shear moduli, vacancy, or surface energies instead. The discrepancy is absorbed by a putative free electron contribution which does not contribute to the interatomic atomic forces in a con­stant volume ensemble calculation.