How good is damped molecular dynamics as a method to simulate radiation damage in metals?

Classical molecular dynamics (MD) is a frequently used technique in the study of radiation damage cascades because it provides information on very small time and length scales inaccessible to experiment. In a radiation damage process, energy transfer from ions to electrons may be important, yet there is continued uncertainty over how to accurately incorporate such effects in MD. We introduce a new technique based on the quantum mechanical Ehrenfest approximation to evaluate different methods of accounting for electronic losses. Our results suggest that a damping force proportional to velocity is sufficient to model energy transfer from ions to electrons in most low energy cascades. We also find, however, that a larger rate of energy transfer is seen when the ionic kinetic energy is confined to a focused sequence of collisions. A viscous damping coefficient dependent on the local atomic environment is shown to be an excellent model for electronic energy losses in low energy cascades in metals.


Introduction
In a radiation damage process, incoming high energy particles collide with the ions of the solid. The primary knock-on atoms (PKAs), hit by the incoming particles, recoil and go on to cause damage themselves [1]. At present, numerical study of the evolution of the damage due to these PKAs typically uses classical molecular dynamics (MD) simulations. The crucial advantage of MD is that the evolution of the ions can be studied on experimentally inaccessible picosecond timescales and nanometre length scales. MD simulations have established that the damage process first involves a cascade of colliding ions, causing ionic displacements up to some point of maximum damage. The behaviour then becomes qualitatively different as kinetic and potential energy is repartitioned and the majority of the ions return to lattice sites. The focus of interest is on the density and structure of the remaining defects as these influence the material's long term behaviour. For reviews of a wide range of MD results see [2,3].
The propagation of a swift PKA through a solid also causes a response by the electronic system [4]. Valence electrons are excited into higher energy states and the electronic density reacts to dynamically screen the PKA; the electrons act as an energy sink and also provide an interatomic potential through which the ions move.
Numerical evolution of quantum mechanical electrons and ions on the time and length scales of a radiation damage cascade is impossible with current computing facilities. Attempts are therefore made to incorporate the required electronic effects within MD phenomenologically, often by treating the electrons as hidden classical degrees of freedom. In this paper, we introduce a method, based on quantum mechanical Ehrenfest dynamics [5], to investigate the use of such an approach. We examine the energy exchange between the ionic and electronic subsystems and find that, for the majority of the low energy cascades simulated, a damping force proportional to velocity is an effective way to include the energy transfer. However, in those cascades where the ionic kinetic energy is confined to a focused sequence of collisions, we see a higher rate of energy transfer and thus would require a larger damping force.
We begin in section 2 by highlighting the limitations of current techniques for including electronic effects in MD simulations of radiation damage. The outline of our Ehrenfest approach is given in section 3 and the details of the simulations performed are given in section 4.

3
Our results and the subsequent analysis are covered in section 5. Finally, in section 6, we present our conclusions and discuss possible extensions to our approach.

Including electronic effects in MD
To the extent that empirical interatomic potentials are fitted to a range of material properties, classical MD simulations already include electronic contributions to the ionic forces at the level of the Born-Oppenheimer approximation. The majority of MD simulations of cascades reported in the literature do not include any electronic effects beyond this. The adiabatic forces are conservative and cannot represent a frictional loss of ionic energy to the electron system. However, given the high speed of the ions involved in a radiation damage cascade, nonadiabatic effects are bound to occur and may well influence the evolution of the cascade. A range of different methods has been suggested to include these effects.
At high ionic velocities, the PKA can be treated as a perturbation of the electronic system. In such cases, theories focus on the energy transfer to the electrons per unit distance along the PKA path, a quantity referred to as the electronic stopping power. At lower velocities, the ionic motion is more phonon-like and energy transfer is treated in terms of the electron-phonon interaction. We emphasize that both regimes of energy transfer arise from the interaction between electrons and ions in the time-dependent Schrödinger equation [6], but the distinction is convenient and is made throughout the literature.
A simple practical technique is to treat the electrons as a thermalizing heat bath. Caro and Victoria [7] use this method to include both electronic stopping power effects and the electron-phonon interaction. They model the ions using Langevin equations of motion, introducing a random force and an inhomogeneous viscous damping force. By making the damping force on an ion-dependent on the local electronic density, the energy loss from ions to electrons can be treated over a wide range of ionic energies. The work of Duffy and Rutherford [8,9] extends that of Caro and Victoria, explicitly modelling the transport of heat through the electron gas via a diffusion equation. In the work of Nordlund et al [10,11], a velocity proportional damping force is intended to represent the electronic stopping power and is applied only to ions with kinetic energy higher than 10 eV. Exchanges of energy between ions and electrons in the low-energy electron-phonon regime are neglected.
These three approaches use the theoretical methods of Lindhard et al (see [12] for a history of the technique and relevant references) as a basis for the inclusion of electronic stopping effects in MD. The original work aimed to calculate the range and straggling of swift ions in matter and extensions give increasingly accurate predictions of these quantities [13]- [16]. All such methods couch the problem in terms of a single ion interacting with an electron system in an otherwise static potential. At high energies, where the ion loses most of its energy to electrons, this model is expected to be accurate. For low energies, the Lindhard approach yields a stopping power proportional to velocity, suggesting that the inclusion of such effects in MD by a damping force proportional to velocity might be valid. Furthermore, extensive simulation of cascades within the binary collision approximation (BCA) has shown that an energy loss proportional to ionic kinetic energy is consistent with experimental data for ion implantation ranges and stopping by thin films [17,18]. However, the failure of Lindhard-like calculations to include either the disruption of the lattice or phonons, both of which will play a part for the low velocity ions of a cascade, raises questions about its direct application to radiation damage. 4 Finnis et al [19] study the energy transfer to electrons diffusing randomly through the hot cascade and scattering from the phonon system. This leads to a damping force that is dependent on the local ionic temperature. The method has been applied by a number of other authors [20,21]. There is no attempt to include electron effects beyond the electron-phonon description, although these are known to be important at higher ionic velocities.
The techniques described above all have some theoretical justification, but it is difficult to quantify the shortcomings of the different approaches. Does it matter that Lindhard-like calculations of stopping power do not include disruption of the lattice? How large are the errors introduced by neglecting effects beyond phonon scattering in the work of Finnis et al? The next section introduces an Ehrenfest-based approach that can be used to study these issues.

A semi-classical technique for validating methods of incorporating electronic effects in MD
Attempting to investigate methods of incorporating nonadiabatic electron effects into MD by comparison with the full quantum mechanical evolution of the electrons and ions would be computationally intractable. However, we have previously shown that the quantumclassical Ehrenfest approach accurately approximates the true electron-ion interaction when the temperature of the ionic subsystem is much greater than the temperature of the electronic system [22]. The Ehrenfest approach therefore provides a tractable way of including relevant quantum mechanical effects in a cascade. Our approach will be to compare energy transfer in an Ehrenfest simulation with the transfer in a classical MD simulation in which the ions are forced to move identically.
In the Ehrenfest approach, the ions are modelled as a collection of classical particles moving in a 'mean-field' generated by the electrons. In turn the electrons are modelled using full quantum mechanics evolving in a time-dependent potential due to the ions. We therefore represent the state of the electronic system by some density matrix,ρ(t), and the system of ions by a set of trajectories in position space, R(t) := {R 1 (t), R 2 (t), . . . , R N ion (t)}. For some initial electronic density matrix,ρ 0 , and an initial set of ionic positions and momenta, the evolution of the system is given by the following coupled equations: where R(t) denotes a dependence on all ionic positions. Equation (1) is known as the quantum Liouville equation, the density matrix formulation of the time-dependent Schrödinger equation. For simplicity, we assume the mass of each ion, denoted M, is the same. The electronic Hamiltonian,Ĥ (R(t)), is composed of electronic kinetic and potential energy terms as well as the electron-ion potential. It is this term that is dependent on the ionic positions. The operator −∇ R nĤ (R(t)) corresponds to the force on the nth ion due to its interaction with the electrons. The inter-ion force is given by is the classical potential energy arising from ion-ion interactions. We also denote byÛ (t) the unitary evolution propagator corresponding to the HamiltonianĤ (R(t)), thusρ(t) =Û (t)ρ 0Û + (t).

5
In the absence of electron-electron interactions the total energy of the system is This energy is conserved in Ehrenfest dynamics.
For comparison with MD simulations we must separate the adiabatic portion of the electronic energy from the nonadiabatic part. To do so we must first define, and then calculate, the adiabatic evolution of the electrons. Derivations of adiabatic/Born-Oppenheimer behaviour rely on the small electron-ion mass ratio [23,24]. Thus, in order to reach the adiabatic limit for some given ionic evolution, we might move ions along the same trajectories in position-space with a smaller electronic mass (or equivalently rescale the electronic time [25,26]). Adiabatic evolution is then defined as that obtained in the zero electron mass limit. Assuming that such a limit is valid we denote the corresponding propagator byÂ(t). We then define the nonadiabatic energy transfer in an Ehrenfest simulation by It is this quantity that we compare with the energy loss from the ions of MD simulations including models of electronic damping. We view the incorporation of electron-ion energy transfer in an MD simulation of a cascade as successful if, for a comparable set of ionic trajectories, it gives a similar energy loss to the Ehrenfest evolution. The internal energy of a classical MD simulation, which is conserved only in the microcanonical ensemble, is determined as where the two terms correspond to total kinetic and potential energy. We then ask if the total change in energy, E(t) − E(0), is approximately equal to the nonadiabatic energy transfer from the ions obtained using the Ehrenfest approach, i.e. − E ehr (t).
A direct simulation comparison between equation (4) and its classical analogue is possible only if classical MD and Ehrenfest methods give equal trajectories in ionic position-space. This is precluded by the Lyapunov instability of atomic systems which will exponentially magnify any infinitesimal difference in the two methods. Instead we perform an Ehrenfest simulation and record the resulting ionic trajectories. We assume a potential energy surface exists for which classical MD incorporating nonadiabatic effects through some specified technique produces the same position-space trajectories. We may then compute the energy loss that would occur applying the classical technique, should the Ehrenfest trajectory be followed on this potential energy surface. If the classical MD technique and Ehrenfest simulations give a similar energy change for a range of different simulation times we say that the technique used to damp the ions in the MD simulation accounts for nonadiabatic energy loss.
In the limit where the Ehrenfest approach is accurate, i.e. when the electronic temperature tends to zero, the nonadiabatic effects in the approaches of Nordlund et al [10,11], Caro and Victoria [7] and Finnis et al [19] all reduce to a damping proportional to velocity. When a damping proportional to velocity is used to incorporate nonadiabatic effects in an MD 6 simulation the loss of energy from the classical ions can be written as an integral over time. To see this, we note that in such an MD simulation the evolution of the nth ion obeys where b n (t) is the damping coefficient for the nth ion. By differentiating the total energy in equation (5) with respect to time we obtaiṅ Multiplying equation (6) byṘ n (t) and rearranging, we havė The first classical model we consider sets b n (t) = b, a single time-independent scalar. This is equivalent to the Finnis et al [19] model at zero temperature. The second model sets b n (t) again to a time-independent scalar b, but applied only where 1 2 MṘ n (t) 2 > T cut-off . Nordlund et al [10,11] use a similar model with a kinetic energy threshold T cut-off = 10 eV. The third model uses the environmentally dependent damping coefficient of Caro and Victoria, with b n (t) = A log(αρ 1/3 n + β), where the local electronic density ρ n is estimated using an embedded atom potential [27]. We define a comparable constant for model 3, b = Alog(αρ where ρ 0 is the electronic density computed for an atom in the perfect lattice. b for model 3 is the damping experienced by an atom at a perfect lattice site.
If, for the evolution of an electron-ion system by the Ehrenfest approach, the relation holds, this suggests that the model incorporating a damping proportional to velocity can be used to approximate the nonadiabatic energy exchange. It is worthwhile stating the logic behind our approach. If there exists some potential with a velocity proportional damping that can reproduce the Ehrenfest trajectories and energy loss then equation (9) must hold. Therefore, if equation (9) is not true, such a classical potential does not exist. Yet, if equation (9) is true, then it is only indicative that such a potential may exist, not a signal that it must, and certainly not that it is possible to find. Thus, strictly, our study is focused on ruling out the impossibility of using a velocity proportional damping to include nonadiabatic effects. Nonetheless, validation of the equation (9) is a strong indication that such an approach is effective. We note that this work treats the electron system as an energy sink and is not concerned with the changes in the interatomic forces that may arise as the electrons evolve. The damping force specified in equation (6) is always opposed to the motion of the ion and will miss any directional dependence of the force due to electronic excitation. This is an intrinsic flaw in using a scalar damping coefficient to account for nonadiabatic effects. One may go on to study the nonadiabatic changes in the forces, by changing the observable of interest in equation (4), but this is not done in this paper.

Computational method
We now apply the method discussed above to a set of Ehrenfest simulations of cascades in a metal. In this section, we summarize the computational methods used and the initial conditions of our simulations.

7
To simplify the modelling of the electronic system we use the tight-binding formalism, which approximates the electronic wavefunctions as superpositions of atomic orbitals and substantially simplifies, or even neglects, electron-electron interactions [28]. For a detailed description of the combination of the Ehrenfest approach and the tight-binding formalism see [29].
To ease the computational burden, we base our simulations on the copper parameterization of an empirical tight-binding model proposed by Sutton et al [30]. The model reproduces the structural stability and elastic constants of copper. However, the electronic density of states at the Fermi level, a prominent quantity in calculations of the electron-ion energy transfer [31], is not accurately represented. Therefore we view our results as a quantitative exploration of the nonadiabatic energy exchange rather than an accurate prediction in copper. The use of local charge neutrality to capture the effects of metallic screening is suggested by Sutton et al but we choose to neglect this feature to ease computation of the nonadiabatic energy transfer (see below). Our tight-binding model therefore neglects Hartree-like interactions that arise when there is charge transfer. Tests show that such Coulomb and screening effects are negligible at the velocities we are considering.
We consider a simulation cell of 9 × 8 × 7 face-centred-cubic unit cells containing N ion = 2016 ions in total. Periodic boundary conditions are imposed along the edges of the simulation cell. Previous work suggests that, over the timescales considered, such a simulation cell is sufficiently large to give an energy transfer nearly equivalent to simulations with an electronic system existing in an infinitely large crystal [31]. The ionic system is initially at zero temperature with all ions stationary at their lattice sites. Although finite ion temperatures are likely to have an effect on cascade evolution (see section 6), we believe it sensible to consider the zero temperature case first. The electronic and ionic systems are evolved simultaneously using the fourth order Runge-Kutta algorithm with a time-step of 0.01 fs.
The electronic system begins with a finite temperature of T el = 500 K in a potential set by the initial (perfect) lattice of ions. A small finite temperature improves numerical stability by avoiding discontinuous changes in the occupation of states at the Fermi energy. We have found that our results are insensitive to the choice of electronic temperature, provided k B T el is greater than the average spacing between electronic energy levels and much less than the Fermi temperature. The initial electronic density matrix is therefore given bŷ where |k(0) and ε k (0) denote the eigenvectors and eigenvalues of the initial Hamiltonian H ({R m (0)}) and the multiple of 2 is due to spin degeneracy. The Fermi-Dirac function gives the occupancy of each eigenstate of the electron system. The chemical potential, µ, is chosen so that Tr{ρ 0 } = 2νN ion , where the band-filling, ν, is a parameter of the model. Our strict periodic boundary conditions imply that we include only eigenstates corresponding to the -point of the supercell. Earlier work indicated that multiple k-point sampling is of little benefit when calculating energy transfer [31].
One ion is chosen as the PKA and given some initial recoil energy in a specified direction. A total of 10 recoil energies are modelled, from 100 eV to 1 keV. For each energy 24 different 8 initial directions are considered, chosen such that each direction corresponds to an equal solid angle within the irreducible 1 48 th part of the face-centred-cubic unit cell. The scheme used to make this selection of initial directions was adapted by the authors from the work of Tegmark [32,33]. The different initial directions sample a range of different cascades at each recoil energy. The cascade is evolved for 200 fs 2 .
Because of the computational demands of time-dependent tight binding, we are not able to explore the full evolution of a radiation damage cascade and are restricted in the size of the simulation cell we use. At high energies the resulting ion cascade will propagate across the periodic boundary within the first 50-80 fs and collide with itself, causing unrealistic ion evolution. We do not believe this invalidates our simulations as a basis for exploring the use of a velocity proportional damping in MD. We are interested in comparing the energy lost using a damping force with that given by the Ehrenfest approach and it is of little concern that the ionic evolution is not that of a true cascade, as long as it is representative.
In order to determine the nonadiabatic energy transfer it is necessary to calculate the energy of the adiabatically evolved electron system. In section 3, we defined the adiabatic evolution of the electrons for a given set of ionic trajectories as that which occurs in the zero-electron mass limit.
For the system sizes we consider here, an alternative computational method presents itself when we note that in the limit of zero electron mass an eigenstate of the initial Hamiltonian will evolve into a corresponding eigenstate of the instantaneous Hamiltonian [25,26], i.e.

A(t)|k(0) = |k(R(t)) ,
where |k(R(t)) is an eigenstate ofĤ (R(t)). This then giveŝ Therefore, for a configuration of ions at some instant of a cascade, R(t), we determine the adiabatically evolved electronic density matrix by diagonalizing the corresponding Hamiltonian to get the instantaneous eigenvectors and match them to the corresponding initial occupancies. Problems may arise in determining the instantaneous eigenstate into which each initial eigenstate has evolved. For a time-evolving system, when all symmetries are rapidly broken by the ionic motion, eigenvalue crossings are unlikely to occur and thus the correspondence between instantaneous and initial eigenstates can be determined by their positions in an eigenvalue ordered list. Preserving the initial occupations of the eigenvectors yields the adiabatically evolved density matrix.

Simulation results and analysis
In section 3, we argued that a demonstration of the equality of E ehr (t) and E model (t), where both quantities are calculated from the same Ehrenfest simulation, would offer strong evidence that the damping model can be used to incorporate the energy transfer from ions to electrons in MD simulations. We plot figure 1, where E model (t) has been calculated for each of the three classical models considered.
In each classical model the damping coefficient b scales the total energy transfer. We can thus use a linear least squares fitting procedure to see if E model and E ehr are proportional to one another and to calculate an estimator for the constant of proportionality. This constant is then used to determine the damping coefficient b. To analyse the quality of our assumptions and goodness of fit we will use the familiar techniques of linear regression. Analysis of the residuals shows that the difference between E ehr (t) and E model (t), caused by the complicated interplay of electronic forces, has an autocorrelation time small compared to the simulation time. Furthermore, when averaged over initial PKA directions, the expected magnitude of the residuals does not vary significantly with time (except close to t = 0 when the residuals must, of course, be nearly zero). We begin by performing the fitting procedure for each simulation individually, i.e. for a given initial direction and recoil energy. The fitting procedure therefore involves a total of 80 points (since E ehr (t) is calculated every 2.5 fs) and a total of 10 × 24 × 3 fittings are performed. Figure 2 shows the average value and spread of these parameters as a function of PKA energy. The damping coefficient for model 2 (KE cut-off) is more than twice that for model 1 (constant damping). As slow moving ions are ignored in model 2, energy must be transferred at a higher rate from the fast moving ions to give the same total. The damping coefficient for model 3 (inhomogeneous damping) is lower than that for model 1, because the coefficients correspond to the damping of an ion at an equilibrium lattice site and most of the energy transfer in a cascade is from ions experiencing a higher than normal electronic density.
In a typical cascade the PKA causes a broad region of ions to be displaced from their equilibrium lattice sites. However, an anomalously high damping coefficient is found in the cascades where the PKA is fired directly at one of its nearest neighbours, a single focused sequence of collisions is caused along a close packed direction of the lattice. Such ion behaviour is known as a replacement collision sequence (RCS). In an RCS all the initial energy of the PKA is concentrated in a few ions and they therefore move more rapidly and approach each other more closely than in a generic cascade. The classical phenomenological models we have considered are therefore seen to fail to capture the true dependence of the electronic damping on the atomic environment and direction of motion. We note that our simulation method has highlighted these collision sequences by using a low temperature ionic lattice, and that at higher lattice temperatures such collision sequences will be shorter. If we exclude the 110 RCS, then for PKA energies above 500 eV we see that the computed average damping coefficient for each classical model converges to an energy independent limit. This is the best estimate for a single value ofb which may be used in such a classical simulation. The values ofb together with the spread of the values computed are given in table 1. For each model the damping coefficient for a 110 RCS is significantly higher than for the other cascades. In model 1 (constant damping) the ratio is 1.92. The ratio is lower for model 3, because the increased damping of an RCS is partly taken into account by the electronic density dependence in the model. The ratio is lowest for model 2 as a consequence of the kinetic energy cut-off; the damping coefficient is determined by high energy collisions such as the RCS head-on collisions.
We use the computed value ofb, excluding the 110 RCS to find the energy transfer as a function of time for each cascade using the classical models. The distribution of the difference between these energies and those computed using the Ehrenfest method is then used to construct an R 2 measure of the goodness of fit, shown in figure 3. Excluding the 110 RCS simulations, the R 2 values for models 1 and 3 are seen to converge well at high PKA energies. Model 1 gives a limiting value of R 2 = 0.95, model 3 does slightly better with R 2 = 0.97. In figure 4, we plot the R 2 values computed for model 2 as a function of the kinetic energy cut-off T cut-off at each PKA energy. If the simulation data supported a nonzero T cut-off we should see this as a peak in R 2 at this value. However, we must remember that the introduction of the cut-off also introduces an extra degree of freedom which must allow for some improvement of the fit in some regions. We find the optimal T cut-off to be rather low (0.6 ± 0.1 eV), and that the goodness of fit is not very much better than that using no cut-off, being only raised to R 2 = 0.97. A cut-off value of T cut-off of the order of several electron volts, as suggested in the literature, does not give a good model for the energy transfer rate in our cascades. We note that our model has a rather flat electronic density of states. A more complex electronic Hamiltonian may better justify use of a cut-off if band structure effects are pronounced.

Conclusions and discussion
The quantum mechanical Ehrenfest method is an accurate and efficient technique for computing the nonadiabatic energy transfer between ions and electrons in radiation damage cascades. The methodology we present here is designed to assess the quality of a classical approximation, by finding the goodness of fit of a particular implementation of stopping power to evolution using the time-dependent Schrödinger equation. It is not our present intention to emphasize the damping coefficient computed within each model; a full calculation and comparison to experiment can be found in [22]. In the present work, we have applied the method to compare and assess various classical models for electronic damping used in MD simulations of radiation damage. The three models we consider: an isotropic homogeneous viscous damping coefficient; a viscous damping applied only to fast-moving ions; and an electronic density-dependent inhomogeneous damping, are the low electronic temperature limits of important methods for including nonadiabatic effects in classical MD today. We have shown that the simple viscous damping coefficient, analogous to the approach of Finnis et al is quite adequate to compute the average energy transfer in low-energy cascades in metals. Applying a very low ionic kinetic energy cut-off to the damping rate slightly improves the fit, but this effect is not dramatic and we have no physical justification for the optimized cut-off at 0.6 eV. If a finite kinetic energy cut-off T cut-off is used, then the dynamics will tend towards that of the corresponding microcanonical ensemble when the average KE per atom is small compared with T cut-off . This could be claimed as a practical justification for using such a cut-off. It would be of course be incorrect to equate this to dynamics within the canonical ensemble without applying a thermostatting mechanism such as a Nosé-Hoover or Langevin thermostat to equilibrate electronic and ionic temperatures [34]. The correct handling of electron-ion . The average of the high energy (>500 eV) cascades is shown as the dashed line. A modest improvement in the accuracy of the linear fit, indicating a better match between classical and quantum mechanical models is possible by optimizing the cut-off kinetic energy above which damping is applied to the moving atoms. This optimum we find is T cut-off = 0.6 ± 0.1 eV for the high energy cascades, raising the R 2 goodness of fit from 0.95 to 0.97. equilibration is extremely challenging within a semi-classical formalism. The Ehrenfest method does not give a finite temperature steady state as spontaneous phonon emission is missing from the equations [35]. We have looked only at short timescales (of the order of a few atomic vibration periods), where the missing return of heat from the electrons is negligible. We have also been careful to test classical models with a zero-temperature steady state, so that the comparison with the Ehrenfest steady state is fair. Our method should therefore be seen as a test for validating the damping terms within classical models; we make no claims for validating mechanisms for incorporating energy return from hot electrons to cold ions. We have shown that, while the energy transfer rate averaged over cascade energies and directions is reproduced well by the simple classical models, the instantaneous damping of a particular atom is not. The equivalent damping coefficient for atoms colliding head-on, generating a 110 RCS, should be double the average for glancing collisions. We find that it is possible to improve the classical model by taking some account of the variation in damping with the local atomic environment, as suggested by Caro and Victoria, as this goes some way to accounting for the increased energy transfer during head-on collisions.
Finite computer resources have limited our study to low energy cascades and a simple tightbinding model. With a more sophisticated Hamiltonian the technique described here could be 14 used to find accurate quantitative estimates of damping coefficients for classical MD simulation. Our methodology has not only been able to illustrate deficiencies within existing classical models of electronic damping, but will also act as a benchmark for new classical models as they appear.