Atomic scale molecular dynamics simulations of radiation damage have been performed on beryllium. Direct threshold displacement simulations along a geodesic projection of directions were used to investigate the directional dependence with a high spatial resolution. It was found that the directionally averaged probability of displacement increases from 0 at 35 eV, with the energy at which there is a 50% chance of a displacement occurring is 70 eV and asymptotically approaching 1 for higher energies. This is, however, strongly directionally dependent with a 50% probability of displacement varying from 35 to 120 eV, with low energy directions corresponding to the nearest neighbour directions. A new kinetic energy dependent expression for the average maximum displacement of an atom as a function of energy is derived which closely matches the simulated data.

Beryllium (Be) has been chosen as the first wall material for the ITER nuclear fusion reactor,1 as its low atomic number and vapour pressure minimise radiative losses from the plasma. Furthermore, Be or high Be compounds have been proposed as neutron multipliers for tritium breeding in the future DEMO reactor.2 

Owing to the long history of Be in nuclear applications, its behaviour at elevated temperatures during irradiation has been well characterised experimentally for a thermal neutron flux.3–5 More recently, efforts have been undertaken to evaluate its performance under fusion conditions both experimentally6–8 and through modelling, whereby the fundamental defect processes, such as the interaction of interstitials, vacancies, impurities, and radiogenic H and He, have been characterized.9–12 

One issue that modelling has yet to investigate is the threshold displacement energy (ED) of Be, which is a parameter used in estimations of how many point defects will be created during irradiation. This value may be fed into stochastic Monte-Carlo simulations and approaches based on derivatives of the Kinchin-Pease model, to evaluate the long-term irradiation response of materials.13,14

Molecular Dynamic (MD) simulations were carried out using the Lammps code.15 Through MD, equations of motion for each atom are integrated as a function of time,16–18 with forces on each atom calculated from the gradient of the energy, itself calculated from the atom's interaction with its neighbours. This allows the velocities and positions of the atoms to be calculated, with this procedure repeated at discreet time steps, letting the position of the atoms to evolve over time. Canonical (NVT) and isothermal-isobaric (NPT) ensembles are simulated using Berendsen thermostats and barostats.19 

Interactions between atoms are described using Embedded Atom Model (EAM)20,21 potentials of the form described by Agrawal et al.22 with a short-range cutoff range of 5 Å. Note that some constants were misreported in the original paper which has since been corrected with a corrigendum. This potential was chosen as predictions of Be lattice properties corresponding well to experimental data and density functional theory (DFT) simulations, a selection of which are presented in Table I.

The EAM potential takes the form

Ei=12ijϕ(rij)+F(ρi),

where Ei is the energy of atom i, rij is the separation between atoms i and j, ϕ is the pair interaction energy, and F(ρi) is the embedding function, which represents the energy contribution from the interaction of atom i with electron density ρ. The pair potential used is of the Morse form25 

ϕ(rij)=De(e2α(rijre)2eα(rijre)),

where De is the cohesive energy (0.41246 eV), re is the nearest neighbour bond length (0.36324 Å), and α (0.6324 Å−1) is a parameter related to the curvature about the bottom of the potential well. The density, ρi, has an exponential form

ρ=AeB(rijre),

where A (1.597 e/Å3) and B (0.49713 Å−1) are empirically fitted constants. To both the pair potentials and the electron density function, the Voter taper function26 is applied to limit the interaction range to a cut-off of 5 Å and avoid a discontinuity in the 2nd derivative

ϕtapered=ϕ(r)ϕ(rcutoff)+rcutoff10[1(rrcutoff)10]dϕdrrcutoff.

The embedding function, F(ρi), is the Johnson function27 

F(ρi)=F0[1ln((ρiρ0)β)](ρiρ0)βF1(ρiρ0)γ

in which F0 (2.03930 eV), F1 (−12.6178 eV), β (0.18752), ρo (1.0 e/Å3), and γ (−2.28827) are empirically fitted constants.

Threshold displacement simulations were undertaken by first equilibrating 41 × 41 × 28 Be supercells with dimensions of 94.7 × 94.7 × 100.3 Å at a temperature of 300 K for a minimum of 500 ps with timestep 0.001 fs at constant pressure and temperature (NPT) using a Berendsen thermostat and barostat. Periodic boundary conditions were applied. To produce 20 vibrationally different starting configurations, for threshold displacement simulations, equilibration was continued for additional increments of 2 ps to a maximum of 540 ps.

During simulations, primary knock-on energies (PKE) of 5–200 eV in steps of 5 eV are imparted in the desired crystallographic direction to the Primary Knock-on atom (PKA), which is chosen to be the central atom of the cell. Two regions surrounding the PKA, region i and ii, are defined with respect to the application of the thermostat. Region i is a sphere around the PKA with radius optimised to 40 Å such that significant displacements of atoms due to the PKA do not extend outside the sphere. To this region, no thermostat is applied. Region ii constitutes the remainder of the supercell and has a Berendsen thermostat applied to remove the excess heat due to the PKA, without directly scaling the atomic velocities in region i.

Following the impact of the PKA, the simulation is run in for 20 000 timesteps of 0.01 fs, 10 000 timesteps of 0.1 fs, and for a further 10 000 timesteps of 1.0 fs, giving a total simulation time of 11.2 ps. This time period was previously shown as sufficient to capture both displacement events and recombination of unstable defects.28 It is essential to use such a short timestep during the impact phase due to the high velocity of the PKA and other atoms in the initially damaged region.

The crystallographic directions investigated were those that fall within a 60° clockwise arc from the [101¯0] direction (as plotted in spherical polar coordinates) and a 90° arc from the (0001) plane to the [0001] direction. Directions within this range were generated based on the vertices of a geodesic projection of the sphere, with a uniform spacing of 6° in both φ and ρ, ensuring results can be directionally averaged to be representative of the overall material similar to Robinson et al.29 For each direction and energy, 20 replicas were performed. A general schematic of the cell setup is shown in Fig. 1.

The way the threshold displacement energy, Ed, is calculated is by imparting a single atom in a crystal with incrementally increasing PKE in a single direction until it is displaced from its lattice site to another metastable site.28 This must be repeated several times for each PKE from slightly different system starting configurations to account for the effects of thermal vibrations on the exact atomic positions at the time of impact. The results of these simulations are presented in Figures 2 and 4.

Figure 2 shows the probability that an atom is displaced at a given energy averaged across all the directions investigated. The lower bound for the threshold displacement energy is PD0, which is the lowest PKE that has a non-zero probability of resulting in a displacement. The model predicts this threshold to be 35 eV for this material. Above 35 eV, the probability of displacement approaches 1 asymptotically, having reached 0.99 for a PKE of 200 eV. The PKE for which there is an average 50% probability that an atom is displaced (PD50) is 70 eV. Both the PD0 and PD50 have been considered the criteria for threshold displacement energy, with PD50 being suggested as more suitable for Monte-Carlo codes such as SRIM.29 A further criteria considered are the energy with a threshold displacement probability of 0.1 PD0.1, which is 45 eV. These criteria have been shown to compare favourably with experimental estimations of ED by previous MD studies.28 The commonly used Monte-carlo code SRIM uses ED of 25 eV,29 while recent molecular dynamics simulations of damage cascades in Be estimate it to be 21 eV.30 The large difference is likely attributable to the choice of potential used, as demonstrated for Iron by the work of Nordlund et al.31 There is also a large difference between the simulated PD50 value and that used in SRIM which could have serious implications for the reliability of SRIM calculations in Be. There are no experimental values available to the authors knowledge.

The mean displacement, at the end of the total simulation time (11.2 ps) for a given direction, follows a similar pattern, remaining at approximately 0.45 Å (x0) until 35 eV (i.e., PD0). This value is similar to thermal oscillation of the atoms around their lattice sites. Beyond 35 eV, the displacement increases gradually, reaching approximately 4 Å at 200 eV. We can model maximum displacement as a function of PKE by modelling the material as a continuum force field that exerts a drag force on the PKA. We modelled this in two ways: first, with drag as proportional to the momentum of the PKA, and second, as proportional to the kinetic energy of the PKE, giving the two equations describing the drag force on the PKA

FD=αmdxdtandFD=β2m(dxdt)2,

where α is the drag coefficient in the momentum model and β is the drag coefficient in the kinetic energy dependent model. These can be related to the acceleration through Newton's law to form the differential equations

md2xdt2αmdxdt=0andmd2xdt2β2m(dxdt)2=0,

which given the boundary conditions that at t = 0, dxdt=2PKEm, and x = 0, can be solved to give

x=1α2PKEm(eαmt1)andx=2βln(cc+β2t),wherec=12PD0m.

To calculate the maximum displacement as a function of PKE, at the limit where PKE = PD0, the equations become

x=1α2PKEm(PD0PKE1)andx=2βln(PD0PKE).

These equations are valid when PKE > PD0, below which they have no physical significance. In the absence of a collision event, there will be some thermal motion resulting in a non-zero displacement, x0. Thus, the equations become

x=x0x=x0x=1α2PKEm(PD0PKE1)+x0x=2βln(PD0PKE)+x0}0<PKE<PD0,ED<PKE.

These models are compared to the simulated results in Figure 3, with α optimised to 1.06 fs−1 and β to 1.02. Both models reproduce the simulated results; however, the momentum model (linear drag) overestimates displacement at high PKE, while the energy dependent model is more faithful in this energy range. There is prescience for such an energy dependent correction in the Norgett–Robinson–Torrens (NRT) model which predicts the number of displaced atoms.14 The NRT model introduced an energy dependent efficiency term which had the effect of decreasing the number of defects predicted to form with increasing PKA energy.14 

The directional dependence of PD0 and PD50 is examined in Figure 4. By both measures, ED shows strong directional anisotropy, with PD0 ranging from 30 to 75 eV and PD50 ranging from 35 to 120 eV. Minima in ED are apparent in {101¯0} directions, corresponding to the direction of the nearest-neighbour Be in the basal plane. Further minima are observed in {21¯1¯1}, which correspond to the second nearest Be neighbours out of the plane. Finally, three further minima are apparent surrounding the {0001} directions, which are themselves shallow minima. These relationships hold true for both measures of ED. They also support the previous findings of Thomas et al.,32 who investigated ED in rutile using similar methodology and found that primary knock-on events in the nearest neighbour directions caused a collision sequence resulting in a larger separation of the interstitial—vacancy pair, which was thus more likely to remain stable. Conversely, it was found that glancing-angle collisions resulted in the highest ED as there was little separation between the PKA and vacancy, as well as a larger disordered region which encouraged recombination. Although these results are expected, previous studies have not investigated the directional dependence in such spatial resolution. Such information may be a useful input for binary collision models that take into account the directional dependence of threshold displacement.

We have investigated the threshold displacement energy and its directional dependency using MD in conjunction with classical EAM potentials. It was found that the probability of an atom being permanently displaced increases with primary knock-on energy following a roughly s shaped curve, with a non-zero probability of displacement occurring at 35 eV, and a 50% probability of displacement at 70 eV. The threshold displacement energy was determined to be strongly directionally dependent, with the nearest neighbour directions having a low threshold displacement energy due to the initiation of a collision sequence. Finally, a momentum dependent and an energy dependent model for the maximum displacement of the primary knock-on atom have been developed. The energy dependent model better reproduces the simulated results, especially at the highest displacement energies.

M.L.J. acknowledges Culham Centre for Fusion Energy for sponsorship. Computing resources were provided by the Imperial College London High Performance Computing Service. Michael Rushton is acknowledged for productive discussions.

1.
M. F.
Smith
and
A. W.
Mullendore
, “
An evaluation of beryllium for magnetic fusion components
,”
J. Nucl. Mater.
122
,
855
857
(
1984
).
2.
H.
Yamada
 et al, “
Preliminary neutronic estimation for demo blanket with beryllide
,”
Fusion Eng. Des.
69
,
269
273
(
2003
).
3.
R. K.
Klueh
,
R. E.
Stoller
,
D. S.
Gelles
, and
H. L.
Heinisch
, “
Neutron damage in beryllium
,”
J. Nucl. Mater.
191-194
,
194
198
(
1992
).
4.
M.
Dalle Donne
,
F.
Scaffidi-Argentina
,
C.
Ferrero
, and
C.
Ronchi
, “
Modelling of swelling and tritium release in irradiated beryllium
,”
J. Nucl. Mater.
212–215
,
954
960
(
1994
).
5.
V. I.
Dubinko
and
V. R.
Barabash
, “
Theoretical assessment of irradiation swelling in beryllium
,”
J. Nucl. Mater.
233–237
,
832
836
(
1996
).
6.
J.
Reimann
and
H.
Harsch
, “
Thermal creep of beryllium pebble beds
,”
Fusion Eng. Des.
75–79
,
1043
1047
(
2005
).
7.
Y.
Mishima
 et al, “
Recent results on beryllium and beryllides in Japan
,”
J. Nucl. Mater.
367–370
,
1382
1386
(
2007
).
8.
V.
Chakin
,
J.
Reimann
,
A.
Moeslang
,
R.
Latypov
, and
A.
Obukhov
, “
Thermal conductivity of highly neutron-irradiated beryllium in nuclear fusion reactors
,”
Prog. Nucl. Energy
57
,
2
7
(
2012
).
9.
P. A.
Burr
,
S. C.
Middleburgh
, and
R. W.
Grimes
, “
Crystal structure, thermodynamics, magnetics and disorder properties of Be–Fe–Al intermetallics
,”
J. Alloys Compd.
639
,
111
122
(
2015
).
10.
S. C.
Middleburgh
and
R. W.
Grimes
, “
Defects and transport processes in beryllium
,”
Acta Mater.
59
,
7095
7103
(
2011
).
11.
P.
Zhang
,
J.
Zhao
, and
B.
Wen
, “
Retention and diffusion of H, He, O, C impurities in Be
,”
J. Nucl. Mater.
423
,
164
169
(
2012
).
12.
D.
Duffy
, “
Modelling materials for fusion power
,”
Int. Mater. Rev.
56
,
324
340
(
2011
).
13.
G. H.
Kinchin
and
R. S.
Pease
, “
The displacement of atoms in solids by radiation. Reports
,”
Prog. Phys.
18
,
1
51
(
1955
).
14.
M. J.
Norgett
,
M. T.
Robinson
, and
I. M.
Torrens
, “
A proposed method of calculating displacement dose rates
,”
Nucl. Eng. Des.
33
,
50
54
(
1975
).
15.
S.
Plimpton
, “
Fast parallel algorithms for short–range molecular dynamics
,”
J. Comput. Phys.
117
,
1
19
(
1995
).
16.
D. C.
Rapaport
,
The Art of Molecular Dynamics Simulation
(
Cambridge University Press
,
2004
).
17.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Clarendon Press
,
1989
).
18.
D.
Frenkel
and
B.
Smit
, “
Understanding molecular simulation: From algorithms to applications
,”
Comput. Sci. Ser.
1
,
1
638
(
2050
).
19.
H. J. C.
Berendsen
,
J. P. M.
Postma
,
W. F.
van Gunsteren
,
A.
DiNola
, and
J. R.
Haak
, “
Molecular dynamics with coupling to an external bath
,”
J. Chem. Phys.
81
,
3684
(
1984
).
20.
M. S.
Daw
and
M. I.
Baskes
, “
Embedded-atom method: Derivation and application to impurities, surfaces, and other defects in metals
,”
Phys. Rev. B
29
,
6443
6453
(
1984
).
21.
S. M.
Foiles
,
M. I.
Baskes
, and
M. S.
Daw
, “
Embedded-atom-method functions for the fcc metals Cu, Ag, Au, Ni, Pd, Pt, and their alloys
,”
Phys. Rev. B
33
,
7983
7991
(
1986
).
22.
A.
Agrawal
,
R.
Mishra
,
L.
Ward
,
K. M.
Flores
, and
W.
Windl
, “
An embedded atom method potential of beryllium
,”
Model. Simul. Mater. Sci. Eng.
21
,
085001
(
2013
).
23.
M.
Chou
,
P.
Lam
, and
M.
Cohen
, “
Ab initio study of structural and electronic properties of beryllium
,”
Phys. Rev. B
28
,
4179
4185
(
1983
).
24.
A.
Migliori
, “
Beryllium's monocrystal and polycrystal elastic constants
,”
J. Appl. Phys.
95
,
2436
(
2004
).
25.
P. M.
Morse
, “
Diatomic molecules according to the wave mechanics. II. Vibrational levels
,”
Phys. Rev.
34
,
57
64
(
1929
).
26.
A. F.
Voter
,
Los Alamos Unclassified Technical Report No. 93-3901
,
1993
.
27.
R. A.
Johnson
, “
Analytic nearest-neighbor model for fcc metals
,”
Phys. Rev. B
37
,
3924
3931
(
1988
).
28.
M.
Robinson
,
N. A.
Marks
,
K. R.
Whittle
, and
G. R.
Lumpkin
, “
Systematic calculation of threshold displacement energies: Case study in rutile
,”
Phys. Rev. B
85
,
104105
(
2012
).
29.
J.
Ziegler
and
J.
Biersack
, “
The stopping and range of ions in matter
,” in
Treatise on Heavy-Ion Science
(
Springer
,
1985
), pp. 93–129.
30.
V. A.
Borodin
and
P. V.
Vladimirov
, “
Molecular dynamics simulation of atomic displacement cascades in beryllium
,” in
Proceedings of the 11th IEA International Workshop on Beryllium Technology
(
KIT Scientific Publishing
,
2013
), p.
236
.
31.
K.
Nordlund
,
J.
Wallenius
, and
L.
Malerba
, “
Molecular dynamics simulations of threshold displacement energies in Fe
,”
Nucl. Instrum. Methods Phys. Res. Sect. B
246
,
322
332
(
2006
).
32.
B. S.
Thomas
,
N. A.
Marks
,
L. R.
Corrales
, and
R.
Devanathan
, “
Threshold displacement energies in rutile TiO2: A molecular dynamics simulation study
,”
Nucl. Instrum. Methods Phys. Res. Sect. B
239
,
191
201
(
2005
).