The collective behavior of dislocations can be described thanks to dislocation dynamics codes. In order to reinforce the physical foundation, input data such as mobility laws can be obtained from atomistic calculations ofindividual dislocations. These defects can now be investigated using more accurate ab initio electronics structure methods. We exemplify these studies by focusing in the following section on the properties ofdislocations in bcc metals and especially
iron. In these materials, dislocation properties are known to be closely related to their core structure.

When dealing with dislocations, special care should be taken in the positioning ofthe dislocations and in the boundary conditions of the calculations. For instance, considering (111) screw dislocations, the two cell geometries proposed in the literature — the cluster approach85 and the periodic array of dislocation dipoles86 — have been thoroughly compared.87 The calculations of dislocations are extremely demanding as they can include up to 800 atoms, so studies usually use fast codes such as SIESTA.8 The construction of simulation cells appro­priate for such extended defects should be optimized for cell sizes accessible to DFT calculations, and the cell-size dependence of the energetics evidenced in both the cluster approach and the dipole approach for various cell and dipole vectors should be rationa­lized. The quadrupolar arrangement of dislocation dipoles is most widely used for such calculations87 although the cluster approach with flexible boundary conditions can be considered a reference method when no energies are necessary (i. e., only structures).

DFT calculations in bcc metals such as Mo, Ta, Fe, and W85,87-91 predict a nondegenerate structure for the core, as illustrated in Figure 9 using differential displacement maps as proposed by Vitek.92 The edge component reveals the existence of a significant core dilatation effect in addition to the Volterra field, which can be successfully accounted for by an aniso­tropic elasticity model.93

Thanks to good control of energy, it is also possi­ble to obtain quantitative results on the Peierls poten­tial; namely, the 2D energy landscape seen by a straight screw dislocation as it moves perpendicular to the Burgers vector. This is exemplified in the following Figure 10(a), where a high symmetry direc­tion of the Peierls potential is sampled: the line going between two easy core positions along the glide direc­tion, that is, the Peierls barrier. These calculations

were performed by simultaneously displacing the two dislocations constituting the dislocation dipole in the same direction and by using a constrained relaxation method. In the same work, the behavior of the Ackland-Mendelev potential for iron,45 which gives the correct nondegenerate core structure unlike most other potentials, has been tested against the obtained DFT results. It appears that it compares well with the DFT results for the g-surfaces, but discrepancies exist on the deviation from anisotropic elasticity of both edge and screw components and on the Peierls potential. Indeed, the empirical potential results do not predict any dilatation elastic field exerted by the core. Besides, the Peierls barrier dis­played by the Ackland-Mendelev potential yields a camel hump shape, as illustrated in Figure 10(a), and at the halfway position, the core spreads between

two easy core positions, whereas it exhibits a single hump barrier within DFT and a nearly hard-core structure at halfway position. The effect of the exchange-correlation functional within DFT appears to be significant.87

More insight into the stability of the core structure can be gained by looking at the response of the polarization of the core, as represented in Figure 10(b). In the Ackland-Mendelev and DFT cases, these calculations confirm that the stable core is completely unpolarized, and they prove that there is no metastable polarized core.95

Finally, the methodology exists for calculating the structure and formation and migration energies of single kinks, but using it with DFT96 remains chal­lenging because cells with about 1000 atoms are needed, together with a high accuracy.

In the second-moment approximation to tight bind­ing, the cohesive energy is proportional to the square root of the bandwidth, which can be approximated as a sum of pairwise potentials representing squared hopping integrals. Assuming atomic charge neutral­ity, this argument can be extended to all band occu­pancies and shapes22 (Figure 7).

The computational simplicity of FS and EAM fol­lows from the formal division of the energy into a sum of energies per atom, which can in turn be evaluated locally. Within tight binding, we should consider a local density of states projected onto each atom. The preceding discussion of FS potentials concentrates solely on the d-electron binding, which dominates transition metals. However, good potentials are diffi­cult to make for elements early in d-series (e. g., Sc, Ti) where the s-band plays a bigger role. An extension to the second-moment model, which keeps the idea of

locality and pairwise functions, is to consider two separate bands, for example, s and d.

This was first considered for the alkali and alka­line earth metals, where s-electrons dominate. These appear at first glance to be close-packed metals, forming fcc, hcp, or bcc structures at ambient pres­sures. However, compared with transition metals, they are easily compressible, and at high pressures adopt more complex ‘open’ structures (with smaller interatomic distances). The simple picture of the physics here is of a transfer of electrons from an s — to a d-band, the d-band being more compact but higher in energy. Hence, at the price of increasing their energy (U), atoms can reduce their volumes (V). Since the stable structure at 0 K is determined by minimum enthalpy H = U + PV at high pressure, this s—d transfer becomes energetically favorable. The net result is a metal-metal phase transformation characterized by a large reduction in volume and often also in conductivity, since the s-band is free electron like while the d-band is more localized. Two-bands potentials capture this transition, which is driven by electronic effects, even though the crystal structure itself is not the primary order parameter.

Materials such as cerium have isostructural tran­sitions. It was thought for many years that Cs also had such a transition, but this has recently been shown to be incorrect,23 and the two-band model was origi­nally designed with this misapprehension in mind.24

For systems in which electrons change, from an s-type orbital to a d-type orbital as the sample is pressurized, one considers two rectangular bands of widths W1 and W2 as shown in Figure 8 with widths evaluated using
eqn [3]. The bond energy of an atom may be written as the sum of the bond energies of the two bands on that atom as in eqn [4], and a third term giving the energy of promotion from band 1 to band 2 (see eqn [8]):

Ubond ni1(ni1 N1 )

2Ni

і 1

+ ni2(ni2 — N2) + Eprom [6]

2 N2

where N1 and N2 are the capacities of the bands (2 and 10 for s and d respectively) and ni1 and ni2 are the occupa­tion of each band localized on the ith atom.

For an ion with total charge T, assuming charge neutrality,

Пі1 + na = T [7]

The difference between the energies of the band cen­ters a1 and a2 is assumed to be fixed. The values of a correspond to the appropriate energy levels in the isolated atom. Thus, a2 — a1 is the excitation energy from one level to another. For alkali and alkaline earth metals, the free atom occupies only s-orbitals; the promotion energy term is therefore simply

Eprom = n2(a2 — a1) = n2Eo [8]

where E0 = a2 — a1.

Thus, the band energy can be written as a function of ni1, ni2, and the bandwidths (evaluated at each atom as a sum of pair potentials, within the second — moment approximation). Defining,

Vi = n,1 — n,2 [9]

and using eqn [7], we can write as follows:

D(E)

N/W1 + N/W2

N/W, —

s-band

Figure 8 Schematic picture of density of electronic states in rectangular two-band model. Shaded region shows those energy states actually occupied.

Ubond = — Wa) — T (Wn + W,2)

‘ v2 + T2 /Wn WA

8 N1 + N2

ViT Wn_ Wi2 + 4 N1 N2

+ T—1 Eo [10]

Although this expression looks unwieldy, it is com­putationally efficient, requiring only two sums of pair potentials for Wj and a minimization at each site independently with respect to vi, which can be done analytically.

In addition to the bonding term, a pairwise repul­sion between the ions, which is primarily due to the screened ionic charge and orthogonalization of the valence electrons, is added. In general, this pair potential should be a function of vi and Vj. But to maintain locality, one has to write this pairwise con­tribution to the energy in the intuitive form, as the sum of two terms, one from each ‘band,’ proportional to the number of electrons in that band:

V(r, j ) = (nil + nj 1) Vi (rjj) + (ni2 + nj2) V2 (rjj) [11]

The variational property expressed in eqn [14] can be exploited to derive the force on the ith atom:

r dUtot

j drj

_ dUtot _ 9Utot @V = dr, lv dr,

Hence, the force is simply the derivative of the energy at fixed v. Basically, this is the Hellmann- Feynman theorem25 which arises here because v is essentially a single parameter representation of the electronic structure.

This result means that, like the energy, the force can be evaluated by summing pairwise potentials. Hence, the two-band second-moment model is well suited for large-scale MD. The force derivation itself is somewhat tedious, and the reader is referred to the original papers. There is no Hellman-Feynman type simplification for the second derivative, so analytic expressions for the elastic constants in two-band models are long ranged and complicated. Conseq­uently, elastic constants are best evaluated numerically.

the mechanism may involve shearing rather than isostructural collapse, particularly if a continuous interface between the two phases exists, as in a shock­wave. With the two-band model, the transition is first order, the volume collapse occurring before the bulk modulus becomes negative in the unstable region. Although the shear and tetragonal shear decrease in the unstable region, neither actually goes negative.

The two-band model is applicable to transition metals, but since the d-band is occupied at all pres­sures, electron transfer is continuous and there is no phase transition. This makes the empirical division of the energy into s and d components challenging. However, once appropriately scaled for ionic charge and number of electrons T in principle, the method could be extended to alloys with noninteger T The results of such an extrapolation are extraordinarily good (Figure 10), considering that there is no fitting to any material other than Cs. The extrapolation breaks down at high Z where the amount of sp hybri­dization is not fully captured in the parameterization. As with FS, no information about band shape is included, and so the sequence of crystal structures cannot be reproduced.

While the extrapolated potentials do not represent the optimal parameterization for specific transition
metals, the recovery of the trends across the group lends weight to the idea that the two-band model correctly reproduces the physics of this series.

The s-d two-band approach has also been applied with considerable success by considering the s-band as an alloying band.26 This has been applied to the FeCr system, which we discuss in more detail later.

Radiation-induced defects, being obstacles to disloca­tion glide, increase yield and flow stress and reduce ductility (see Chapter 1.04, Effect of Radiation on Strength and Ductility of Metals and Alloys for

experimental results). Furthermore, ifthe obstacle den­sity is sufficiently high to block dislocation motion, pre­existing Frank-Reed dislocation sources are unable to operate and plastic deformation requires operation of sources that are not active in the unirradiated state. These new sources operate at much higher stress and give rise to new mechanisms such as yield drop, plastic instability, and formation of localized channels with high dislocation activity and high local plastic deforma­tion. Understanding these phenomena is necessary for predicting material behavior under irradiation and the design and selection ofmaterials for new generations of nuclear devices.

Obstacles induced by irradiation affect moving dislocations in a variety of ways, but can be best categorized as one of two types, namely inclusion­like obstacles and those with dislocation properties.

The first type includes voids, bubbles, and preci­pitates, for example. They usually have relatively short-range strain fields and their properties may not be changed significantly by interaction with disloca­tions. (Copper precipitates in iron are an exception to this — see Sections 1.12.4.1.1-1.12.4.1.2.) Those that are not impenetrable are usually sheared by the dislo­cation and steps defined by the Burgers vector b of the interacting dislocation are created on the obstacle — matrix interface. Unstable precipitates, such as Cu in Fe, may also suffer structural transformation during the interaction, which can change their properties. These obstacles do not usually modify dislocations signifi­cantly, although they may cause climb of edge dis­locations (see Sections 1.12.4.1.1-1.12.4.1.2). Their main effect is to create resistance to dislocation glide. Obstacles such as voids and bubbles are among the strongest, and as a result of their high density, they contribute significantly to radiation-induced harden- ing.7 Materials designed to exploit oxide dispersion strengthening (ODS) are produced with a high con­centration ofrigid, impenetrable oxide particles, which introduce extremely high resistance to dislocation motion.8 These obstacles are also considered here as their scale, typically a few nanometers, is similar to that of obstacles formed under irradiation.

The second obstacle type consists of those with a dislocation character, for example, DLs and SFTs, and so dislocation reactions occur when they are encountered by moving dislocations. Loops have rel­atively long-range strain fields and hence interact with dislocations over distances much greater than their size. SFTs are three-dimensional (3D) struc­tures and have short-range strain fields. Loops with perfect Burgers vectors are glissile, in principle, whereas SFTs and faulted loops, for example, Frank loops in fcc metals, are sessile. In addition to causing hardening, the reaction of these defects with a gliding dislocation can modify both their own structure and that of the interacting dislocation. As will be demon­strated in Section 1.12.4.2, their effect depends very much on the geometry of the interaction, that is, their position and orientation relative to the moving dislo­cation, and the nature of the mutual dislocation seg­ment that may form in the first stage of interaction. The contribution of these obstacles to strengthening can be significant, for their density can be high.

Modification of irradiation-induced microstructure due to plastic deformation is an additional possibly important effect. If mechanical loading occurs during irradiation, it can contribute significantly to the overall microstructure evolution and therefore to change in material properties. Accumulation of internal stress during irradiation is unavoidable in real structural materials and so this effect should not be ignored. The effects of concurrent deformation and irradiation on microstructure are far from clear, for only a few exper­imental studies of in-reactor deformation have been performed.9 This area, therefore, provides a good example of how atomic-scale modeling can help in understanding a little-studied phenomenon.

Many observations contradict the FP3DM. These include the void lattice formation11-14 and higher swelling rates near GBs than in the grain interior in the following cases: high-purity copper and alumi­num irradiated with fission neutrons or 600 MeV protons (see original references in reviews117,118); aluminum irradiated with 225 MeV electrons119 and neutron-irradiated nickel120 and stainless steel.121 Furthermore, the swelling rate at very low dislocation density in copper is higher,122-124 and the depen­dence of the swelling rate on the densities of voids and dislocations is different,125 than predicted by the FP3DM. It gradually became clear that some­thing important was missing in the theory. There was evidence that this missing part could not be the effect of solute and impurity atoms or the crystal structure. Indeed, austenitic steels of significantly different compositions and swelling incubation peri­ods exhibit similar steady-state swelling rates of ~ 1% per NRT dpa.3 , And, although generally the bcc materials show remarkable resistance to swelling, , the alloy V-5% Fe showed the highest swelling rate of ^2% per dpa: 90% at 30dpa.34

As outlined in Section 1.13.3.1, the primary dam­age production under neutron and ion irradiations is more complicated; in addition to PDs, both vacancy

and SIA clusters are produced in the displacement cascades. This is the reason the FP3DM predictions fail to explain microstructure evolution in solids under cascade damage conditions. In fact, it has been shown that it is the clustering of SIAs rather than vacancies that dominates the damage accumu­lation behavior under such conditions. The PBM proposed in the early 1990s and developed during the next 10years (see Section 1.13.1) essentially resolved many of the issues; the phenomena men­tioned have been properly understood and described. This model is described in the next section.

For an anisotropic, linearly elastic crystal, Mura29 derived a line integral expression for the elastic dis­tortion of a dislocation loop, as

uij(x) E/»kCpqmnb?

where nk is the unit tangent vector of the dislocation loop line L, d/ is the dislocation line element, jh is the permutation tensor, Cijk/ is the fourth-order elas­tic constants tensor, Gy,/(x — x0) = dGy (x — x’)/dxi, and Gij(x — x0) are the Green’s tensor functions, which correspond to displacement components along the xi direction at point x due to a unit point force in the xj direction applied at point x0 in an infinite medium.

It can be seen that the elastic distortion eqn [10] involves derivatives of the Green’s functions, which need special consideration. For general anisotropic solids, analytical expressions for Gij, k are not available. However, these functions can be expressed in an integral form (see e. g., Barnett,3 Willis,3 Bacon

et a/.,32 and Mura33), as

Gij, k(x x ) [—?kNi^- (k)D—1 (k)

Ck

+ kk Clpmq (rpkq + kprq )Nil (k)Njm(k)D 2(k)]df [11]

where r = x — x0, r = r/r, k is the unit vector on the plane normal to r, the integral is taken around the unit circle Ck on the plane normal to r, and Nj(k) and D(k) are the adjoint matrix and the determinant of the second-order tensor Cikjlkkkl, respectively.

At a point P on a dislocation line L, the external force per unit length is obtained by the Peach — Koehler formula FA = (b-oA) x t, where sA is the sum of the stress from an applied load and that arising from internal sources at P, b is the Burgers vector of the dislocation, and t is the unit tangent vector of the element dl. For the special, yet important case of a

-o— 3-Nodes ■o…. 4-Nodes -o-— 5-Nodes
L/4 L/4 L/4 L/4 Figure 2 Nodal displacements for the first time step of an initially straight segment.
with only two linear and equal segments connected at point A, as shown in Figure 2. We will compute the shape of the line, advancing it from its initial straight configuration to a curved position. Under these sim­plifications, the variation in Gibbs free energy 8G for any one of the two segments is given by
1 ft dr |ds| 0
V 8r |ds
B	[15]
8G
Now, we expand the virtual displacement and veloc- ity in only two shape functions: C1 = u, C2 = 1 — u. Thus
drk = dqtkC,	[16]
Vk = q, k,t C, [17]
Since we allow the displacement to be only in a direction normal to the dislocation line (y direction), we drop the subscript k as well. For arbitrary varia- tions of 8qik, the following equation is applicable to any of the two segments.
At x (Pk + fs)C, |ds| = — B
DqmCmCi | ds|	[18:
Equation [18] can be explicitly integrated over a short time interval At. The resistivity matrix elements are defined by g, m = BjC, Cm|ds|, and the force vector 0 1 elements by f = J C, x (fpK + fs) | ds | . With these definitions, we have the following (2 x 2) algebraic system for each of the two elements: Aqmgim = f x At [19]

For any one linear element, the line equation can be determined by

And the resistivity matrix can be simplified as [21]

Furthermore, as a result of the shear stress t and the absence of self-forces during the first time step only, the distributed applied force vector reads

tbi 1 "2 1

Since the dislocation line is divided into two equal segments, we can now assemble the force vector, stiff­ness matrix, and displacement vector in the global coordinates, and arrive at the following equation for the global nodal displacements AQ:

F1

+ At 0

F3

An important point to note here is that, at the two fixed ends, we know the boundary conditions but the reaction forces needed to satisfy overall equilibrium are unknown. These reactions act on the fixed obsta­cles at L and R, and are important in determining the overall stability of the configuration (e. g., if they exceed a critical value, the obstacle is destroyed, and the line is released). If AQ = AQ = 0 at both fixed ends, we can easily solve for the nodal displacement AQ = 3 lbAr and for the unknown reaction forces at the two ends: F1 = F3 = 1|rbl. If we divide the disloca­tion line into more equal segments, the size ofthe matrix equation expands, but the nodal displacements and reaction forces can be calculated similarly. Results of analytical solutions for successively larger number of nodes on the dislocation segment are shown in Figure 2.

In the following section, we present several applica­tions of the computational method to irradiated mate­rials, where modeling defect behavior has contributed to a greater understanding of radiation effects.

In an interdiffusion experiment, a sample A (mostly composed of A atoms) is welded to a sample B (mostly composed of B atoms) and annealed at a temperature high enough to observe an evolution of the concentration profile. According to eqn [12], the flux of component i in the reference crystal lattice is proportional to its concentration gradient:

J = — D, ГС, [18]

where the so-called intrinsic diffusion coefficient Di is a function of the phenomenological coefficients and the thermodynamic factor:

LV

D = <f — C F [19]

An interdiffusion experiment consists of measure­ment of the intrinsic diffusion coefficients as a func­tion of local concentration. The resulting intrinsic diffusion coefficients are observed to be dependent on the local concentration. Within the TIP, while the driving forces are locally defined, the L-coefficients are considered as equilibrium constants. It is not easy to ensure that the experimental procedure satisfies these TIP hypotheses, especially when concentration gradients are large, and the system is far from equilib­rium. When measuring diffusion coefficients, one implicitly assumes that a flux can be locally expanded to first order in chemical potential gradients around an averaged solid solution defined by the local con­centration. Starting from atomic jump frequencies and applying a coarse-grained procedure, a local expansion of the flux has been proved to be correct in the particular case of a direct exchange diffusion mechanism.62

An interdiffusion experiment is not sufficient to characterize all the diffusion properties. For example, in a binary alloy with vacancies, in addition to the two intrinsic diffusion coefficients, another diffusion coefficient is necessary to determine the three inde­pendent coefficients LAA, LAB, and LBB.

Calculations of defect production, eqn [2], require knowledge of the damage function, n(T). While it is not possible to measure this function directly, as no irradiation creates monoenergetic recoils except near the surface, it can be obtained by measuring defect production for a wide range of ion irradiations and subsequently deconvoluting eqn [3]. Low-energy light ions, for example, weight the recoil spectrum near the threshold energy, «25-100 eV, while more energetic heavy ions weight it at high energies. Results are shown for Cu in Figure 13. Here, electri­cal resistivity measurements are employed to monitor the absolute number of FPs produced per unit dose of irradiation. Included in this figure are the damage efficiency function, X(T1/2), deduced from the experi­ments and X(T) calculated using molecular dynamics computer simulation. The damage efficiency function is defined as

v(T ) = £(T )vNRT (T), [20]

where nNRT(T) is the NRT damage function defined by eqn [1]. The good agreement between experiment and simulations illustrates that the damage function in Cu is now well understood. This is now true for many other pure metals as well.24 In alloys and ceramic

materials, however, the damage function remains poorly known.

Our simulation system is typically a parallelepiped supercell in which particles are placed either in a very regular manner, as in modeling a crystal lattice, or in some random manner, as in modeling a gas or liquid. For the simulation of perfect crystals, the number of particles in the simulation cell can be quite small, and only certain discrete values, such as 256, 500, and 864, should be specified. These numbers pertain to a face — centered-cubic crystal that has four atoms in each unit cell. If our simulation cell has l unit cells along each side, then the number of particles in the cube will be 4l3. The above numbers then correspond to cubes with 4, 5, and 6 cells along each side, respectively.

Once we have chosen the number of particles we want to simulate, the next step is to choose the system density we want to study. Choosing the density is equivalent to choosing the system volume since density p = N/O, where N is the number of particles and O is the supercell volume. An advantage of the Lennard-Jones potential is that one can work in dimensionless reduced units. The reduced density pa3 has typical values of about 0.9—1.2 for solids and 0.6-0.85 for liquids. For reduced temperature kB T/e, the values are 0.4-0.8 for solids and 0.8-1.3 for liquids. Notice that assigning particle velocities according to the Maxwellian velocity distribution probability = (m/2nkB T)3/2exp[-m( v2 + vj + ^)/2кв T]dvx dvy dvz is tantamount to setting the system temperature T.

For simulation of bulk properties (system with no free surfaces), it is conventional to use the periodic boundary condition (PBC). This means that the cubical simulation cell is surrounded by 26 identical image cells. For every particle in the simulation cell, there is a corresponding image particle in each image cell. The 26 image particles move in exactly the same manner as the actual particle, so if the actual particle should happen to move out of the simulation cell, the image particle in the image cell opposite to the exit side will move in and become the actual particle in the simulation cell. The net effect is that par­ticles cannot be lost or created. It follows then that the particle number is conserved, and if the simulation cell volume is not allowed to change, the system density remains constant.

Since in the pair potential approximation, the particles interact two at a time, a procedure is needed to decide which pair to consider among the pairs between actual particles and between actual and image particles. The minimum image convention is a procedure in which one takes the nearest neighbor to an actual particle as the interaction partner, regardless of whether this neighbor is an actual parti­cle or an image particle. Another approximation that is useful in keeping the computations to a manage­able level is the introduction of a force cutoff distance beyond which particle pairs simply do not see each other (indicated as rc in Figure 4). In order to avoid a particle interacting with its own image, it is necessary to set the cutoff distance to be less than half of the simulation cell dimension.

Another book-keeping device often used in MD simulation is a neighbor list to keep track of who are the nearest, second nearest, … neighbors of each particle. This is to save time from checking every particle in the system every time a force calculation is made. The list can be used for several time steps before updating. In low-temperature solids where the particles do not move very much, it is possible to do an entire simulation without, or with only a few, updating, whereas in simu­lation ofliquids, updation every 5 or 10 steps is common.

If one uses a naive approach in updating the neighbor list (an indiscriminate double loop over all particles), then it will get expensive for more than a few thousand particles because it involves N x N operations for an N-particle system. For short-range interactions, where the interatomic potential can be safely taken to be zero outside of a cutoff rc, acceler­ated approaches exist that can reduce the number of operations from order-N2 to order-N. For example, in the so-called ‘cell lists’ approach,18 one partitions the supercell into many smaller cells, and each cell maintains a registry of the atoms inside (order-N operation). The cell dimension is chosen to be greater than rc, so an atom cannot possibly interact with more than one neighbor atom. This will reduce the number of operations in updating the neighbor list to order-N.

With the so-called Parrinello-Rahman method,1 the supercell size and shape can change dynamically during a MD simulation to equilibrate the internal stress with the externally applied constant stress. In these simulations, the supercell is generally non­orthogonal, and it becomes much easier to use the so-called scaled coordinates sy to represent particle positions. The scaled coordinates Sj are related to the real coordinates ry through the relation, ry = H • sy, when both rj and sj are written as column vectors. H is a 3 x 3 matrix whose columns are the three repeat vectors of the simulation cell. Regardless of the shape of the simulation cell, the scaled coordi­nates of atoms can always be mapped into a unit cube, [0, 1) x [0, 1) x [0, 1). The shape change of the simulation cell with time can be accounted for by including the matrix H into the equation of motion. A ‘cell lists’ algorithm can still be worked out for a dynamically changing H, which minimizes the number of updates.13

For modeling ionic crystals, the long-range elec­trostatic interactions must be treated differently from short-ranged interactions (covalent, metallic, van der Waals, etc.). This is because a brute-force evaluation of the electrostatic interaction energies involves computation between all ionic pairs, which is of the order N2, and becomes very time­consuming for large N. The so-called Ewald sum — mation20,21 decomposes the electrostatic interaction into a short-ranged component, plus a long-ranged component, which, however, can be efficiently summed in the reciprocal space. It reduces the computational time to order N3/2. The particle mesh Ewald22-24 method further reduces the computational time to order N log N.

with t taken to be as long as possible. In terms of discrete time steps, eqn [10] becomes

1 1

(A )=yE A(tk) [11]

L k=1

where L is the number of time steps in the trajectory. The second is a time-dependent quantity of the form

1L0

(A(0)B(t))= A(tk)B(tk + t) [12]

L k=1

where B is in general another dynamical variable, and L’ is the number of time origins. Equation [12] is called a correlation function oftwo-dynamical vari­ables; since it is manifestly time dependent, it is able to represent dynamical information of the system.

We give examples of both types of averages by considering the properties commonly calculated in MD simulation.

N

U = V(ry) potential energy [13]

i<j

N

m Vi • Vi

i=1

@V (rj’) h @rj ij,

1N

g (r )=rwN(|=: j=:d(r

radial distribution function

1N

MSD(t) =Njri(t) ^ri (0)l2

N i = 1

mean squared displacement

In eqn [19], va is the average volume of one atom, va is the a-component of vector Vi, and rja is the a-component of vector ri — ry. The interest in writing the stress tensor in the present form is to suggest that the macroscopic tensor can be decomposed into individual atomic contributions, and thus is known as the atomic level stress25 at atom i. Although this interpretation is quite appealing, one should be aware that such a decomposition makes sense only in a nearly homogeneous system where every atom ‘owns’ almost the same volume as every other atom. In an inhomogeneous system, such as in the vicinity of a surface, it is not appropriate to consider such decomposition. Both eqns [15] and [19] are written for pair potential models only. A slightly different expression is required for potentials that contain many-body terms.26

1.10.13.1 Useful Concepts Without True Physical Meaning

For very large simulations, imaging is a problem, since showing all atoms in a massive simulation is likely to obscure the important regions. There are a number of heuristic quantities arising from the interatomic potential which can be used to pick out the atoms associated with atypical configurations of interest.

• Most empirical potentials define the energy per atom.

• The atomic level stress55

“b 1 a b a b

s = 2O ~ m, v’vP’

1 j

where f is the force on atom i due to atom j which determines the glass transition.56

• The concept of ‘local crystal structure’ can be used to locate twin boundaries, phase transitions, etc. This may be done by common neighbor analysis, or bespoke investigation of pair and angular distribution functions to search for particular configurations.57

• The balance between pair and many-body ener­gies in glue-type models, or more significantly the various angular density functions of MEAM.

• The magnetization in ‘magnetic’ potentials.

Although uniquely defined for a given potential, many of these are not well-defined concepts in electronic structure. Yet, they can be extremely useful in identi­fying the atoms in far-from-equilibrium environments.

1.10.9 Summary

Interatomic potential development is a continuing challenge for materials modeling. They represent the only way to perform MD, which in turn is crucial for the nonequilibrium and off-lattice processes, which dominate radiation damage. Despite best efforts, few potentials can be reliably employed to predict quanti­tative energies beyond where they are fitted. Their most useful role is to reveal processes and topologies that might be of importance in real materials.

1.13.1 Introduction

The study of radiation effects on the structure and properties of materials started more than a century ago,1 but gained momentum from the development of fission reactors in the 1940s. In 1946, Wigner2 pointed out the possibility of a deleterious effect on material properties at high neutron fluxes, which was then confirmed experimentally.3 A decade later, Konobeevsky et al4 discovered irradiation creep in fissile metallic uranium, which was then observed in stainless steel.5 The discovery of void swelling in neutron-irradiated stainless steels in 1966 by Cawthorne and Fulton6 demonstrated that radiation effects severely restrict the lifetime ofreactor materials and that they had to be systematically studied.

The 1950s and early 1960s were very productive in studying crystalline defects. It was recognized that atoms in solids migrate via vacancies under thermal — equilibrium conditions and via vacancies and self­interstitial atoms (SIAs) under irradiation; also that the bombardment with energetic particles generates high concentrations of defects compared to equilib­rium values, giving rise to radiation-enhanced diffu­sion. Numerous studies revealed the properties of point defects (PDs) in various crystals. In particular, extensive studies of annealing of irradiated samples resulted in categorizing the so-called ‘recovery stages’ (e. g., Seeger7), which comprised a solid basis for understanding microstructure evolution under irradiation.

Already by this time, which was well before the discovery of void swelling in 1966, the process of interaction of various energetic particles with solid

targets had been understood rather well (e. g., Kinchin and Pease8 for a review). However, the primary dam­age produced was wrongly believed to consist of Frenkel pairs (FPs) only. In addition, it was com­monly believed that this damage would not have serious long-term consequences in irradiated materi­als. The reasoning was correct to a certain extent; as they are mobile at temperatures of practical interest, the irradiation-produced vacancies and SIAs should move and recombine, thus restoring the original crys­tal structure. Experiments largely confirmed this sce­nario, most defects did recombine, while only about 1% or an even smaller fraction survived and formed vacancy and SIA-type loops and other defects. How­ever small, this fraction had a dramatic impact on the microstructure of materials, as demonstrated by Cawthorne and Fulton.6 This discovery initiated extensive experimental and theoretical studies of radiation effects in reactor materials which are still in progress today.

After the discovery of swelling in stainless steels, it was found to be a general phenomenon in both pure metals and alloys. It was also found that the damage accumulation takes place under irradiation with any particle, provided that the recoil energy is higher than some displacement threshold value, £d, (^30-40 eV in metallic crystals). In addition, the microstructure of different materials after irradiation was found to be quite similar, consisting of voids and dislocation loops. Most surprisingly, it was found that the microstructure developed under irradiation with ~1 MeV electrons, which produces FPs only, is similar to that formed under irradiation with fast neutrons or heavy-ions, which produce more compli­cated primary damage (see Singh et a/.1). All this created an illusion that three-dimensional migrating (3D) PDs are the main mobile defects under any type of irradiation, an assumption that is the foundation of the initial kinetic models based on reaction rate the­ory (RT). Such models are based on a mean-field approximation (MFA) of reaction kinetics with the production of only 3D migrating FPs. For conve­nience, we will refer to these models as FP produc­tion 3D diffusion model (FP3DM) and henceforth this abbreviation will be used. This model was devel­oped in an attempt to explain the variety of phenom­ena observed: radiation-induced hardening, creep, swelling, radiation-induced segregation (RIS), and sec­ond phase precipitation. A good introduction to this theory can be found, for example, in the paper by Sizmann,9 while a comprehensive overview was pro­duced by Mansur,10 when its development was already completed. The theory is rather simple, but its general methodology can be useful in the further development of radiation damage theory (RDT). It is valid for ^1MeV electron irradiation and is also a good introduction to the modern RDT, see Section

1.13.5.

Soon after the discovery of void swelling, a number of important observations were made, for example, the void super-lattice formation11-14 and the microm­eter-scale regions of the enhanced swelling near grain boundaries (GBs).15 These demonstrated that under neutron or heavy-ion irradiation, the material micro­structure evolves differently from that predicted by the FP3DM. First, the spatial arrangement of irradia­tion defects voids, dislocations, second phase particles, etc. is not random. Second, the existence of the micrometer-scale heterogeneities in the microstruc­ture does not correlate with the length scales accounted for in the FP3DM, which are an order of magnitude smaller. Already, Cawthorne and Fulton6 in their first publication on the void swelling had reported a nonrandomness of spatial arrangement of voids that were associated with second phase precipi­tate particles. All this indicated that the mechanisms operating under cascade damage conditions (fast neu­tron and heavy-ion irradiations) are different from those assumed in the FP3DM. This evidence was ignored until the beginning of the 1990s, when the production bias model (PBM) was put forward by Woo and Singh.16,17 The initial model has been changed and developed significantly since then18-2 and explained successfully such phenomena as high swelling rates at low dislocation density (Section

1.13.6.2.2) , grain boundary and grain-size effects in void swelling, and void lattice formation (Section

1.13.6.2.3) . An essential advantage of the PBM over the FP3DM is the two features of the cascade dam­age: (1) the production of PD clusters, in addition to single PDs, directly in displacement cascades, and (2) the 1D diffusion of the SIA clusters, in addition to the 3D diffusion of PDs (Section 1.13.3). The PBM is, thus, a generalization of the FP3DM (and the idea of intracascade defect clustering intro­duced in the model by Bullough et a/. (BEK29)). A short overview of the PBM was published about 10 years ago.1 Here, it will be described somewhat differently, as a result of better understanding of what is crucial and what is not, see Section 1.13.6.

From a critical point of view, it should be noted that successful applications of the PBM have been limited to low irradiation doses (< 1 dpa) and pure metals (e. g., copper). There are two problems that prevent it from being used at higher doses. First, the PBM in its present form1 predicts a saturation of void size (see, e. g., Trinkaus et at19 and Barashev and Golubov30 and Section 1.13.6.3.1). This originates from the mixture of 1D and 3D diffusion-reaction kinetics under cascade damage conditions, hence from the assumption lying at the heart of the model. In contrast, experiments demonstrate unlimited void growth at high doses in the majority of materials and conditions (see, e. g., Singh et at.,31 Garner,32 Garner et at.,33 and Matsui et at.34). An attempt to resolve this contradiction was undertaken23,25,27 by including thermally activated rotations of the SIA-cluster Burgers vector; but it has been shown25 that this does not solve the problem. Thus, the PBM in its present form fails to account for the important and common observation: the indefinite void growth under cascade irradiation. The second problem of the PBM is that it fails to explain the swelling satura­tion observed in void lattices (see, e. g., Kulchinski et at.13). In contrast, it predicts even higher swelling rates in void lattices than in random void arrange — ments.25 This is because of free channels between voids along close-packed directions, which are formed during void ordering and provide escape routes for 1D migrating SIA clusters to dislocations and GBs, thus allowing 3D migrating vacancies to be stored in voids.

Resolving these two problems would make PBM self-consistent and complete its development. A solu­tion to the first problem has recently been proposed by Barashev and Golubov35,36 (see Section 1.13.7). It has been suggested that one of the basic assumptions of all current models, including the PBM, that a random arrangement of immobile defects exists in the material, is correct at low and incorrect at high doses. The analysis includes discussion of the role of RIS and provides a solution to the problem, making the PBM capable of describing swelling in both pure metals and alloys at high irradiation doses. The solu­tion for the second problem of the PBM mentioned above is the main focus of a forthcoming publication by Golubov et at.37

Because of limitations of space, we only give a short guide to the main concepts of both old and more recent models and the framework within which radiation effects, such as void swelling, and hardening and creep, can be rationalized. For the same reason, the impact of radiation on reactor fuel materials is not considered here, despite a large body of relevant experimental data and theoretical results collected in this area.