Energetics and reactivity of small beryllium deuterides

Enthalpies and free energies of reaction for small neutral and charged beryllium deuterides BeD, BeD2, and BeD3 that have been calculated are reported for a temperature range of 0 K to 1000 K. We discuss probable dissociation channels and possible ways of producing BeD by localizing the relevant transition states and by calculating corresponding rate constants. BeD and BeD+ are found to be the most stable ones among the considered compounds. BeD2 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathrm{BeD}}_2^{+} $$\end{document}BeD2+ are more likely to decompose into Be0,+ + D2 than into BeD0,+ + D. The metastable BeD3 and \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$ {\mathrm{BeD}}_3^{+} $$\end{document}BeD3+ predominantly decompose into BeD0,+ + D2. In light of our results on the reaction energetics, we can interpret the pathways for production of BeD via BeD2 and BeD3 intermediates observed in molecular dynamics simulations.


Introduction
The development of new technologies for controlled fusion caused beryllium compounds, especially hydrides that can be created by D/T bombardment from plasma, to become one focus of materials research. In the ITER reactor, beryllium is planned to be the first-wall material, and hence it will be directly exposed to particles, predominantly deuterons, that escape the confinement as is already observed in the JET tokamak with ITER-like walls [1,2]. For many years, plasma-wall interactions (PWI) have been extensively studied experimentally as well as theoretically. The main source of experimental data concerning beryllium-deuterium interactions are linear devices like PISCES-B [3,4] or tokamaks [5][6][7]. However, the underlying processes like sputtering, transport and deposition are hard to reproduce and quantify experimentally. Modeling and theoretical approaches to obtain data for codes like Wall-DYN [8], ERO [9] or SDTrimSP [10,11] are therefore developed to get insight into such processes. Theoretical studies of plasma wall interactions involve the description of the interaction of surfaces with the fusion plasma [12], the characterization of elementary processes [13] as well as the validation of experimental results [14,15]. Concerning beryllium experiments, there is evidence of the formation of BeD molecules and of a linear drop of the BeD:Be ratio with increasing temperature in the temperature range of 500-700 K while no larger molecules like BeD 2 or BeD 3 were observed [7]. In contrast, molecular dynamics (MD) simulations employing analytical bond-order potentials (ABOP) [16,17], as well as multiscale modeling extrapolated from them [18] predict BeD 2 and BeD 3 as the main eroded species for the same This paper belongs to Topical Collection 7th Conference on Modeling & Design of Molecular Materials in Trzebnica (MDMM 2016) temperature range [19]. Since at lower temperatures (<500 K) MD simulations were in agreement with experiment, a more complete description of the possible fragmentation processes is needed from both the energetic and kinetic points of view. Dissociation and reactivity of beryllium hydrides and their isotopes were briefly discussed by Safi et al. [19] and Virot et al. [20] on the basis of standard thermodynamic data and possible reaction channels were given for the dissociation of BeD 2 and BeD 3 at various temperatures. The dissociation and ionization rates for primary reactions of BeD due to electron collisions were reported by Björkas et al. [18]. In the present work, we employ quantum chemical methods including comparisons between various levels of theory and we study the reaction kinetics to determine the reactivity of different channels. Highly accurate data have been published before on the beryllium hydrogen systems: BeH as well as BeH 2 have received a significant amount of interest as a test system for quantum chemical methods including non-standard ones [21][22][23][24][25][26][27][28][29]. The multi-reference averaged coupled-pair functional method (MR ACPF) [21] was used to calculate the accurate ground state potential energy functions, vibration-rotation energy levels for BeH, BeD, and BeT and their ions, which agree excellently with spectroscopic experimental data, i.e., the equilibrium bond length R e = 1.341 Å [22,23]. Non-Born-Oppenheimer variational calculations employing explicitly correlated Gaussian basis functions were performed in order to determine the ionization energy of BeH and the dissociation energies of BeH and BeH + [24]. Penotti's [26] non-orthogonal single and multi-configurational calculations with a highly optimized even-tempered STO basis set yielded a value of R e = 1.329 Å for the D ∞h geometry of the BeH 2 molecule. A value of 2053.0 cm −1 was obtained for the harmonic symmetric-stretch frequency. Very precise results for the BeH system also included non-adiabatic effects and extrapolation of the basis set up to the spdfgh level as well as extrapolation of correlation effects to the full configuration interaction (FCI) limit [28]. An equilibrium distance of R e = 1.341 Å and a ground state frequency of ω e = 2062.1 cm −1 were reported. Hinze et al. [29] published potential energy surfaces (PESs) for BeH 2 and BeH þ 2 obtained with the multireference configuration interaction method (MRCI) and documented the insertion reaction of Be into H 2 . Koput and Peterson [30] obtained vibrational and rotational energy levels of beryllium dihydride and of its isotopes from an accurate potential energy surface using CCSD (T) and extrapolation to the full basis set limit. The IR emission spectra for BeH and BeH 2 were measured by Bernath and coworkers [31][32][33][34]. They obtained R e = 1.342 Å, ω e = 2061.4 cm − 1 for BeH, and R e = 1.326 Å, ω e = 2255.2 cm −1 for the asymmetric stretch for BeH 2 .
Reaction enthalpies for the dissociation channels of the BeD 3 molecule were calculated in ref. [19] and the thermodynamic stability of neutral and anionic BeH 3 was analyzed in ref. [35]. The knowledge of the whole reaction network is necessary for understanding the chemical behavior of the Be/H system. This work aims to describe the fragmentation and reactivity of small beryllium deuterides in the temperature range 0-1000 K and is based on quantum-chemistry calculations and transition state theory. We report a stability analysis, thermodynamic data of neutral and charged BeD 1-3 molecules, the standard enthalpies and free energies of reaction for their possible dissociation channels, and their corresponding dissociation energies. Calculated transition states and activation energies can be used to estimate reaction rate constants. Furthermore, reaction schemes for the production of BeD from beryllium surfaces exposed to D irradiation as extracted from MD simulations are discussed.

Quantum-chemical calculations
The optimized structures and vibrational frequencies of the beryllium hydrides were obtained by the Gaussian-4 (G4) [36] method and by density functional theory (DFT). We compared different DFT functionals: the often used hybrid functional B3LYP-D [37,38], the B97D functional and the meta-GGA M06 functional [39]. The first two functionals contain Grimme's GD3 empirical dispersion parameters [40]. We also employed the double hybrid B2PLYPD [41,42] functional which includes dispersion and density corrections by second-order Møller-Plesset perturbation theory (MP2) [43]. Furthermore, we optimized the predicted structures by coupled cluster calculations with single and double substitutions and noniteratively included triple excitations (CCSD(T)) [44] with all electrons correlated. All calculations were performed with Dunning's correlation consistent core-valence quadruple zeta (aug-cc-pCVQZ) basis set [45], the sole exception being the pre-defined G4 method, and employed the GAUSSIAN 09 software package [46]. Different scaling factors are recommended for some of the methods used throughout this study [47]. However, we report unscaled frequencies here with the intention to introduce as little empiricism as possible.

Thermodynamics
In order to qualitatively examine the reactivity of small beryllium deuterides (hydrides), we first calculated the standard enthalpies and the standard free energies of the reactants and products of the dissociation channels (Eq. 1). The standard enthalpies (the free energies) of reaction were calculated as a difference of the sum of electronic energy ε 0 , zero-point energy ϵ ZPE and thermally corrected enthalpies (free energies) of products and reactants [48]: and similarly for Δ r G°(T). Only calculated values were used to obtain Δ r H°and Δ r G°. Subsequently, equilibrium constants, K EQ , were obtained from the standard free energies of reactions: The standard enthalpies of formation of molecules at 298.15 K, Δ f H°(M), were calculated using experimental enthalpies of formation for elements, Δ f H Å 0K (X) and their corresponding thermal corrections, H 0→298:15 corr ; (see Table 1) and the procedure suggested by McQuarrie [48]: where D 0 is the dissociation energy which is equal to the atomization energy for a number of x atoms of the type X in molecule M. We compared the selected methods by comparing thermodynamic data of the BeH and the BeH 2 molecules with experimental data from the NIST-JANAF database [49].
The results are shown in Table 2. All methods predict values of Δ f H°(BeH) rather close to the experimental value of 321 ± 30 kJ mol -1 (= 3.3 ±0.3 eV) [49], with CCSD(T) differing the most by ∼19 kJ mol -1 (0.2 eV). Δ f H°(BeH 2 ) is significantly overestimated compared to the experimental value of  Table 2 Comparison of calculated thermodynamic properties as obtained by various methods for BeH, BeH 2 , and H 2 (D 2 ) molecules to known experimental values. Δ f H°is in kJ mol -1 , R e in Å, and ω e in cm −1 . The enthalpy of formation for BeH and its estimation for BeH 2 are taken from ref. [50]. Experimental values for the hydrogen molecule are taken from the NIST-JANAF thermochemical Tables [49]. Spectroscopic values for BeH and BeH 2 are retrieved from ref. [32,33]. ω a e for BeH 2 corresponds to the asymmetric stretching vibration and ω b e to the bending vibration

Method
BeH ∼125.52 kJ mol -1 (1.3 eV) which is based on an empirical method [49,50]. For BeH 2 , the DFT functionals overall give values closer to the experimental value than the higher order methods, which differ by ∼33-43 kJ mol -1 (0.3-0.5 eV). The M06 functional predicts enthalpies of formation closest to the reference values. B3LYP-D and CCSD(T) harmonic frequencies are closest to the experimental data. The harmonic frequency for BeD obtained with B3LYP-D and CCSD(T) methods, 1530.4 cm −1 , agrees well with the experimental value of 1529.5 cm −1 [32]. Eventually, to obtain the total change of the free energy of reaction, a term dependent on initial concentrations or more generally on instantaneous activities a v j j of products and a v i i of reactants, referred to as reaction quotient, Q r , has to be added [48]: The ratio Q r / K EQ determines the direction of reaction: if Q r > K EQ , the reaction favors the reactants; if Q r < K EQ , the products are preferable. The reaction is in equilibrium for Q r = K EQ . The reaction quotients were not calculated nor otherwise included for studied dissociation channels in this work.

Localizing transition states and calculating rate constant
Approximate transition states geometries were at first guessed instead of using computational methods such as QST2 [51]. These structures were then optimized using the B3LYP-D functional. It was checked by vibrational frequency analysis that a transition state has only one imaginary frequency, with modes corresponding to the reaction path. In addition, IRC calculations were performed to ensure that the obtained transition states connect the local minima on the PES which refer to reactants and products for the considered reactions. Subsequently, obtained structures were optimized by the CCSD(T) method. The rate constants were determined only for dissociation channels with log(K EQ ) > −5 for any point in the temperature range from 0 to 1000 K because those are the reactions most affecting plasma-wall interactions. In case of reactions with log(K EQ ) < −5 nearly only reactants will appear in the equilibrium mixture [48]. We used transition state theory (TST) to estimate the reaction rate constants for temperature T by employing the Eyring-Polanyi Eq. [52]: where k B is the Boltzmann constant, R the gas constant, h the Planck constant, c the concentration, n the order of reaction, and ΔG ‡ is the free energy of activation. We set c to 1 for results in the present work. The linear form of this equation (Eq. 6.b), where ΔH ‡ and ΔS ‡ are the enthalpy and entropy of activation, was used to present the calculated rate constants. From the rate constants, we are able to determine the reaction schemes for  Table 3 Bond lengths and angles of the beryllium hydrides. The lengths are given in Å, angles in degrees. For the meaning of R 1 and R 2 see Fig. 3. The data for BeH and BeH 2 are given in Table 2 BeH unmixed reactants and products in their standard states at the pressure of 1 atm and thus can predict the feasibility of the studied reactions based on the data from computational electronic structure methods and the rules of chemical kinetics.

Molecular dynamics simulations
We studied the sputtering of BeD by low energy D irradiation from pure Be surfaces by means of molecular dynamics (MD) simulations using the same procedure as in ref. [15]. The D + bombardment was simulated with the DL_POLY 3.9 software [53] which was extended to include ABOP potentials [54]. The details and parameters of the Be-H potentials are given in ref. [17]. The hexagonal closed packed Be surface (0001) with 3718 atoms (30×30×40 Å) was equilibrated by slowly heating the samples to 300 K at a rate of 50 K/ps. Subsequently, 1000 cumulative D impacts with 7, 10, and 20 eV were performed from a distance of 5 Å perpendicular to the center of the surface. A single impact lasted 7 ps and was divided into two parts: the first 3 ps consist of the impact itself followed by 4 ps of relaxation of the cell to remove extra energy from the system. Each step lasted 0.5 fs. The surface was randomly shifted in x-and ydirections after each impact. We compared the sputtering yields with other work (see Fig. 1) and extracted data about single sputtering events to look closer at the mechanisms.

Stability analysis
We obtained optimized geometries of neutral and charged BeH, BeH 2 , and BeH 3 molecules from the various functionals, the G4, and CCSD(T) methods. Concerning BeH 2 the optimized structures all have negative electron affinities with an absolute value of 5-32 kJ mol -1 (0.1-0.3 eV), i.e., an energy is required to attach an electron. They are thus thermodynamically unstable and were removed from further analysis. Be, H, and H 2 are also included in this analysis. The structural properties of H 2 , BeH, and BeH 2 are given in Table 2, the ones for the remaining molecules are summarized in Table 3. Dissociation energies and enthalpies of formation for neutral and positive ions of beryllium deuterides are given in Tables 4 and 5, respectively. The bond lengths for BeH range from 1.338 to 1.367 Å, with CCSD(T) and M06 being closest to the experimental value of 1.342 Å [32]. The bond length of the respective cation is shorter by 0.03 ± 0.01 Å on average, whereas the bond length of the respective anion is longer by 0.08 ± 0.03 Å. The Be-H bond lengths of neutral BeH 2 range from 1.323 to 1.337 Å. Again, the CCSD(T) and M06 values agree excellently with the experimental bond length of 1.326 Å [33]. All methods predict two bent structures for BeH þ 2 , with H-Be-H angles of ∼90°(I) and ∼24°(II), and a multiple saddle point  Fig. 2. This is in contrast with the potential energy surface of BeH þ 2 produced at CMRCI/cm 3 -pVTZ level of theory in ref. [29], where only the linear asymmetric structure is reported beside the van der Waals minimum. However, only BeH þ 2 II is below the dissociation limit for Be + + H 2 . Furthermore, we optimized the equilibrium geometries of the other local minima for BeH 2 and BeH þ 2 of the same publication [29] with B3LYP-D and CCSD(T) to compare and validate our approach for transition state search. The the results are in good agreement (see Table 6), except for the fact mentioned above concerning the non-existence of an asymmetric linear structure of BeH þ 2 for B3LYP-D. All methods yield similar structures for neutral and ionic BeH 3 molecules (see Fig. 3). The angle formed between H-Be-H is notably different with ∼152°for neutral molecules, ∼166°for cations and 120°for anions. The experimental BeD bond length and bond length in BeD 2 are 1.342 and 1.326 Å, respectively [32,33].

Thermodynamics
Standard enthalpies Δ r H°and the free energies Δ r G°of reaction were calculated for all possible dissociation channels of the stable neutral and ionic beryllium deuterides in the temperature range from 0 to 1000 K. We selected reaction channels (Eqs. 7 and 8) with log(K EQ ) > −5 for a further analysis of reaction pathways. The temperature dependences of the free energies of reaction are shown in Fig. 4 for the dissociation of neutral and cationic beryllium deuterides as calculated with B3LYP-D and CCSD(T). We calculated the enthalpies and the free energies of reactions for BeD and BeD + based on very accurate data extracted from MR ACPF calculations by Koput [22] which serve as benchmark values for our results. Δ r H°at 298.15 K and Δ r G°at 298.15 K and 1000 K obtained by different methods are provided for neutral molecules and cations in Tables 7 and 8, respectively. The MR ACPF enthalpies and free energies of the reactions 7.a and 8.a are very close to G4 and CCSD(T) values yielding differences up to 5 kJ mol -1 for G4 and CCSD(T) indicating that BeD and BeD + do not yield strong multi-referential character of the wave functions  in equilibrium. However, this will be discussed in more detail in the next section which is dedicated to transition states.
BeD g ð Þ→Be g ð Þ þ D g ð Þ ð7:aÞ :bÞ Almost all studied dissociation channels have equilibrium constants very close to 0 (K EQ < < 1); the reactants dominate in the mixtures and an increase of the concentration of products leads to the production of more reactants. Reaction 7.c has a rather high Δ r H°value of 395 kJ mol -1 for G4 and 411 kJ mol -1 for B97D (4.1-4.3 eV) at 298.15 K. BeD 2 more likely dissociates into Be and D 2 with Δ r H°(298.15 K) = 188.0 kJ mol -1 for B3LYP-D. The change of standard free energy predicts log(K EQ ) > −5 for more than ∼800 K. For BeD 3 , both channels (7.d and 7.e) have similar characteristics with Δ r G°( T) close to 0 kJ mol -1 already at 0 K. Δ r G°(T) for dissociation of BeD lies in the range of 203-237 kJ mol -1 (2.1-2.5 eV) at 298.15 K. BeD + has the highest Δ r G°(T) of the reported cations. 8.b is the preferable channel of the two most likely ways of the dissociation of BeD þ 2 (Δ r G°< 0 already at 300 K). The second channel, reaction 8.c, has a lower Δ r H°than its neutral alternative: 76-106 kJ mol -1 at 298.15 K (0.7-1.1 eV) with log(K EQ ) < −5 at ∼500 K. There is only one channel (8.d) with Δ r H°= 91-100 kJ mol -1 (0.9-1.0 eV) for BeD þ 3 as the others all have a log(K EQ ) < −5. This reaction has log(K EQ ) > −5 from ∼500 K on.

Reactivity of beryllium deuterides
The structural properties and the free energy of activation ΔG ‡ of transition states obtained using the CCSD(T) method for the forward and reverse reactions along each channel are presented in Tables 9 and 10, respectively. Their structures are depicted in Fig. 5. Reaction rate constants can be calculated

BeD and BeD +
We did not find any transition state concerning the dissociation of BeH or BeH + in line with earlier studies [22][23][24].
BeD 2 and BeD þ 2 A more complex behavior is found for the reactivity of neutral and cationic BeD 2 . Beryllium dihydride can dissociate into Be + D 2 through the transition state TS 1, which corresponds to the one found in ref. [29]. No transition state was found for the other channel (7.c) in the ground state, but we localized one (TS 2 ) in the triplet state. CCSD(T) yields the same transition state structures as B3LYP-D. Both transition states are very close in energy and in the region where the ground state potential energy surface intersects the one of the triplet state. Therefore, their activation free energies are similar and rather high, ∼380 kJ mol -1 (3.9 eV), with regard to the BeD 2 ground state. TS 2 is about 55 kJ mol -1 (0.6 eV) higher than the triplet state local minima and about 30 kJ mol -1 (0.3 eV) below the plateau of the excited Be + D 2 complex. We assume that the channel resulting in Be + D 2 is preferable in general, however, in a fusion plasma environment the required excitation energy is easily reachable and dissociation into BeD + D is therefore values are calculated from data in ref. [22]. All values are in kJ mol -1 . The threshold temperature T k for crossing the log(K EQ ) > −5 limit is also given Method BeD(g)→ Be(g) + D(g) BeD 2 (g)→ Be(g) + D 2 (g) BeD 3 (g)→ BeD 2 (g) + D(g) Δ r H°Δ r G°Δ r H°Δ r G°Δ r H°Δ r G°Δ r H°Δ r G°Δ r H°Δ r G°T ---already at 0 ∼200 Table 8 Changes of the standard enthalpy of reaction at 298.15 K and changes of the standard free energy of reaction at 298.15 and 1000 K for dissociation of cationic beryllium deuterides obtained from various methods. MR ACPF values are calculated from data in ref. [23]. All values are in kJ mol -1 . The values for the bent structure of BeD þ 2 are used here. The threshold temperature T k f or crossing log(K EQ ) > −5 limit is also presented

Method
BeD also possible. The difference of ΔG ‡ for the two transition states, ∼20 kJ mol -1 (0.2 eV) in the range 0-1000 K, favors decomposition into BeD + D. The reverse reaction 7.b has a high barrier as well, ∼260 kJ mol -1 (2.7 eV) at 298.15 K and ∼320 kJ mol -1 (3.3 eV) at 1000 K. We could not identify any transition state for channels 8.b and 8.c with CCSD(T). Still, BeD þ 2 is predicted to dissociate into Be + + D 2 due to the negative free reaction energy already at 300 K.
The BeD 3 molecule is metastable. We found a two-step reaction mechanism for reaction 7.d. However, the intermediate (IM) and the transition state TS 4 (see Fig. 4) are very close in energy and similar in structure and are also lower in energy than the local minima at higher temperatures. Thus, they make the decomposition into BeD + D 2 more likely than the competing reaction 7.e. The free activation energies, ΔG ‡ , for these barriers are ∼14 kJ mol -1 (0.1 eV) at 298.15 K and ∼ −0.5 kJ mol -1 (0.01 eV) at 298.15 for TS 3 and TS 4 , respectively. B3LYP-D predicts the same characteristics of the transition states as CCSD(T). We did not identify any transition state for the reaction in Eq. 8.d. In fact, BeD 3 + seems to dissociate most likely into BeD + + D 2 as this channel is energetically preferable (Δ r G°∼0 at 1000 K).
The transition states found using B3LYP-D and CCSD(T) were also investigated for their multi-referential character using the D 1 and T1 diagnostics [57,58]. The results are in excellent agreement with those of higher-order correlation methods if D 1 < 0.03. If D 1 < 0.05, the method still performs well. However, a multi-reference character of the ground-state introduced by strong orbital relaxation effects is indicated by larger values of D 1 . Similarly, if T 1 > 0.02, the system should be investigated by a multi-reference electron correlation method. The conclusion of the D 1 and T 1 diagnostics for BeH molecules is as follows: all equilibrium structures have D 1 less than 0.03 and T 1 < 0.02, except asymmetric linear BeD þ 2 (T 1 = 0.02,D 1 = 0.08), thus CCSD(T) should describe their ground states reliably. Far from equilibrium and for some transition states these diagnostics yield higher values, indicating that single-reference methods could become inadequate for describing these states. This concerns only transition states for neutral and positive BeH 2 (T 1 ∼0.35, D 1 ∼0.08). Transition states related to BeH 3 molecules yield D 1 < 0.03, with minimum and maximum values of 0.015 and 0.029, respectively, and T 1 ∼0.01. Further investigations are required to scrutinize how the various PESs and energetics are affected by more accurate correlation and multi-referential character, which we plan to do in a following work.

Production of BeD
Analysis of our MD simulations yields that sputtering of BeD can often be described by the following reactions (9.a and 9.b) Table 9 Bond lengths and angles of the transition states obtained by CCSD(T) for studied dissociation channels in Eqs. 7 and 8. The lengths are given in Å, angles in degrees. These structures are presented in Fig  BeD 2 g ð Þ þ D g ð Þ→ BeD 3 g ð Þ→BeD g ð Þ þ D 2 g ð Þ ð9:aÞ Be g ð Þ þ D 2 g ð Þ→BeD 2 g ð Þ→BeD g ð Þ þ D g ð Þ ð9:bÞ The reaction 9.b without the intermediate BeD 2 was also suggested and analyzed by Nishijima et al. [55]. However, the ionic reaction Be + + D 2 was shown to dominate the production of BeD + inside a plasma column. Based on our calculated standard free energies of reaction and reaction schemes for BeD 2 and BeD 3 we can determine the reactivity of the proposed reactions in Eq. 9. BeD 2 and D can easily form the BeD 3 intermediate dissociating subsequently into BeD + D 2 . Impinging particles with small kinetic energy could possibly lead to BeD and D 2 via this reaction pathway. The other reaction, 9.b, has a high barrier for all studied temperatures, but it is a possible reaction nonetheless. Concerning the production of BeD from the surface, reaction 9.b is slower than 9.a and higher impact energies are needed.

Conclusions
We report reaction schemes of the dissociation of small beryllium deuterides based on calculated thermodynamic properties, standard enthalpies of formation, and standard enthalpies and free energies of reaction as obtained by quantum-chemical methods. Transition state theory was used to determine the rate constants for the considered reactions. BeD and BeD + are the most stable species, unlikely to further dissociate into their components. BeD 2 and BeD þ 2 are more likely to decompose into Be + D 2 than into BeD + D. BeD 3 and BeD þ 3 are metastable against their dissociation into BeD + D 2 . Concerning the source of beryllium hydride production, we performed MD simulations of low energy D irradiation on Be surfaces to obtain the details of the sputtering events and analyzed these events from thermodynamic and kinetic points of view. The analysis of the MD trajectories confirms that the formation of BeD occurs along the reaction pathways that have been suggested before.