An improved model of interatomic forces for large simulations of metals containing excited electrons

We derive a form for the non-conservative damping forces on metal ions due to their interactions with electrons, and present the result in the second-moment tight-binding approximation suitable for direct and efficient inclusion in a large-scale molecular dynamics simulation. We demonstrate that this form accurately captures the direction, velocity, temperature and local atomic environment dependence of the non-adiabatic force in quantum mechanical simulations in which electronic stopping is accurately calculated. No previous empirical damping force is able to reproduce this rich behaviour.


Introduction
An understanding of dynamic processes involving fast metal ions requires an accurate model for energy exchange between ions and electrons. In sputtering processes, the local non-uniform electronic temperature distribution can influence the sputtering yield [1,2]. In radiation damage cascades, ionic kinetic energy is absorbed by the electrons and transported out of the damage region [3,4]. Duffy et al investigated the effect of electron-ion energy exchange on the evolution of radiation damage phenomena. Their results suggest that this energy exchange can significantly impact the residual defect population following a collision cascade [3,5] and that the strength of the electron-ion interaction is a key determinant of the formation of tracks around channelled ions [6]. For ions channelling between crystal lattice planes with kinetic energies of tens or hundreds of keV, the principal energy loss mechanism is to the electrons [7,8]. All these processes take place on atomic timescales, making experimental observation difficult, and so classical molecular dynamics (MD) simulations are an important means of gaining insight.
In classical MD simulations of metals, the electrons are treated implicitly in the empirical potentials. Non-adiabatic effects, where included, are commonly modelled as an electronic friction-an additional force anti-parallel to the ion velocity [3,5,6], [9]- [13]. The use of such viscous drag forces is often justified with reference to models of electronic stopping power based on highly idealized treatments of a homogeneous electron gas [14] or of isolated binary collisions [15,16]. Such treatments clearly neglect any environmental dependence or directionality of the non-adiabatic forces by construction and so, while they can successfully predict the statistical behaviour of energy exchange processes [8], they cannot be expected to capture any microscopic detail. We should view the use of a drag force as a convenient means of effecting the removal of energy from the ionic subsystem on average over the course of a radiation damage event, rather than as an accurate model of the non-adiabatic forces exerted on the ions by the electrons.
We have previously demonstrated that the non-adiabatic force on a moving ion will have a significant dependence on the local atomic environment, on the direction of motion and on the electronic temperature [17]. In cascade simulations, we observed that a viscous drag force cannot capture the different rates of damping of different modes of ionic motion [18]. In this paper, we derive a form for the non-adiabatic force in a metal that reproduces the microscopic detail of the quantum mechanical force and can be incorporated at very low cost within a classical MD simulation.
In section 2, we derive a model for the non-adiabatic electronic damping force acting on an ion in a metal, starting with the quantum mechanical evolution of electrons coupled to classical ions. In the first instance we obtain a general, basis-independent, expression for the force. We then specialize to a tight-binding picture to arrive at a form suitable for inclusion in MD. In section 3, we demonstrate that our damping model used with an empirical potential reproduces extremely closely the results of time-dependent quantum mechanical simulations down to the detailed forces on individual atoms. Section 4 considers the implications of our new model for MD simulations of radiation damage in metals. 3 Finnis-Sinclair form, the total binding energy is written as the sum of repulsive and attractive components, where the former is a pairwise interionic repulsion and the latter is proportional to the square root of the second moment of the local density of states, computed from the sum of the squares of the hopping integrals (see for example [20]) where r is the 3N vector of all atomic positions and p the momentum, r ab ≡ r b − r a is the separation between atoms a and b and N a is the set of atoms in the neighbourhood of a, i.e. within the range of the potential. For the convenience of notation we assume that all ions have the same mass m, and write the pairwise contributions to the repulsive energyŨ ab ≡Ũ (|r ab |) and the hopping integrals H ab ≡ H (|r ab |). For this potential, the conservative force on an atom is The central result of this paper is our corresponding expression for an additional nonconservative damping force where x is a constant for fine tuning the damping to experimental data and √ 0 is the contribution to the attractive component of the binding energy for an atom in the perfect lattice. The role of the particle velocity in viscous damping models is taken in (3) by the rate of change of bond lengthsṙ ab .
To derive this form, we have considered the following: when an ion moves in a metal, we often assume that the electron response is infinitely rapid (the Born-Oppenheimer approximation), but in reality the ion collides with electrons at the Fermi level producing a shower of electron-hole pairs (excitons) 2 . These excitons give rise to a lag in the bond-orders dependent on the ion velocity. Excitons evolve via the Schrödinger equation, and so when many have been generated, we expect their phases to be incoherent-the net effect of the history of this exciton production is to raise the electronic temperature [21,23]-and only those created recently and in the vicinity of the moving ion will have a significant contribution to the nonadiabatic force. The time locality is captured by introducing an electron scattering time and the space locality is encapsulated in a local projection operator, derived in appendix A.
We derive our result starting from quantum-classical Ehrenfest dynamics, which has been shown to accurately capture the energy transfer from hot ions to cold electrons [24], but does not include the purely quantum mechanical effect in which the electron-hole recombination emits a quantum phonon. We will not attempt to derive an expression to represent spontaneous phonon emission. A productive approach to such a derivation might proceed from the correlated electron ion dynamics (CEID) formalism of Horsfield et al [25,26], but this is beyond the 4 scope of the present paper. Instead, we propose a stochastic form for the spontaneous phonon emission consistent with a Langevin formulation of energy exchange [27], which ensures that an equilibrium is reached between ions and electrons. The total force on an atom in an MD simulation is the Langevin force F a = (F 0 ) a + (F 1 ) a + η a , where η a is a stochastically drawn return force. A compatible form for η a is derived in appendix B.

The non-adiabatic force
To derive the form of non-adiabatic forces, we consider energy transfer processes second order in time via a generalization of the familiar velocity Verlet algorithm. This shows us the direction and magnitude of a force that can be identified as non-adiabatic.
Starting with the quantum-classical non-self-consistent spin-degenerate Hamiltonian whereρ is a single-particle density matrix andĤ is a non-self-consistent one-electron spindegenerate Hamiltonian. U is a repulsive interionic potential deriving from ion-core electrostatic and exchange-correlation energies [28]. 3 Equation (4) differs from the empirical equation (1) in two respects. Firstly, the band energy is a trace over the electronic Hamiltonian; a moment expansion is not taken here. Secondly, the electronic state,ρ, is not necessarily the ground state, but instead evolves with time, as described below. In a time-dependent evolution, electrons and ions are distinct particles and moving one will merely exert a Coulomb force on the other:ρ is not an explicit function of position.
We can write the state at time t as s(t) = (r, p,ρ) T (the superscript denoting transpose) and its time derivative F(s) = (v, F, i/h[ρ,Ĥ ]) T . v = p/m is the ion velocity and F = −∇U − 2 Tr(ρ∇Ĥ ) is the Hellman-Feynman force. This scheme, known as Ehrenfest dynamics, is discussed in detail in [29].
The state at time t + δt is given by s(t + δt) = s(t) + (F(s) + F(s))δt/2, wheres = s + F(s)δt. This time-stepping algorithm, known as the Heun method, reduces to velocity Verlet if the state is a function of position and momentum only.
In order to separate out the non-adiabatic part of the evolution, we split the density matrix into two semi-independently evolving components,ρ ≡ρ 0 +δ, whereρ 0 represents 5 the electronic ground state assumed in classical MD. For transferability, the repulsive part of an empirical potential, which is partly electronic in origin, must be fitted at a specified electronic temperature (corresponding to the temperature associated with the data used to fit the coefficients of the potential). For consistency, the attractive electronic bonding should make the same assumption. Whether an empirical potential is fitted to zero or finite temperature, it is nevertheless the case that MD tacitly assumes the existence of an electronic reservoir, or other electron-electron interactions, which act to instantaneously thermalize the electronic subsystem [30]. The electronic ground state density matrixρ 0 that we use as the reference for the adiabatic evolution is therefore a canonical temperature ground state.
is the Fermi-Dirac occupation at electronic temperature T e , chemical potential µ, of the ith eigenstate of the instantaneous electron Hamiltonian (i.e.
ρ 0 therefore commutes with the instantaneous Hamiltonian, and is an explicit function of the atomic positions.δ is the non-adiabatic remainder, evolving as If we defineF 0 by writing F =F 0 − 2 Tr(δ∇Ĥ ), and the state σ (t) = (r, p,δ) T , then This equation has a particularly simple form when we can assume the non-adiabatic part of the density matrixδ starts negligibly small: The top line is recognized as the velocity Verlet evolution of position and momentum for motion under a conservative potential whose force isF 0 , the Hellman-Feynman force for the 6 ground state. Motion under this force alone, often referred to as Born-Oppenheimer dynamics, conserves free energy is the electronic entropy. Work is done on a reservoir continually resetting the electrons to a ground state. Internal energy can be conserved if the conservative force −2 Tr(∇ρ 0Ĥ ) is added toF 0 , but the difference between evolution in an open or closed system vanishes at zero temperature and is generally neglected in MD. The second line of equation (8) is the non-adiabatic part of the evolution. We see thatδ accumulates in the direction −v · ∇ρ 0 , producing a non-adiabatic force in the direction Tr((v · ∇ρ 0 )∇Ĥ ). In this way,δ is generated from the sum of all the individual electron-hole pairs excited during the evolution history and Ehrenfest dynamics is deterministic and reversible. The term in v · ∇ρ 0 can be interpreted as a contribution toδ due to the evolved density matrixρ(t) lagging behind the changing ground stateρ 0 . However, as each exciton evolves according to the Schrödinger equation, the phases of old excitations are incoherent and only the non-adiabatic density built up recently over a (slowly spatially varying) time τ approximating the electron scattering time will contribute to the non-adiabatic force. This suggests that we can write our non-adiabatic force as The gradient of the density matrix is extremely difficult to compute exactly, and so we need to find an approximation to use equation (9) in an MD simulation. Start by writing the identity where we have used Contributions to the first term will be largest for pairs of states close to the Fermi level so that we can approximate We also have where we have assumed that there is no systematic variation of the energy of eigenstates with position. Thus, we arrive at an approximate form for a metallic system at finite temperature: where indicates a sum that includes only states close to the Fermi level, and we have defined the projection operatorP ≡ i | φ i φ i |.

The non-adiabatic force in a tight-binding picture
To write the non-adiabatic force for a tight-binding potential, we will expand the trace in equation (9) in an orthonormal atomic orbital basis {| a l }, where | a l is the lth orbital on atom a, and define | a = l | a l . In this basis, F 1 ≈ 2τ ab a |v · ∇ρ 0 | b b |∇Ĥ | a , and so we can write the force on a particular atom c as In an atomic basis, we see the non-adiabatic force is made up of contributions acting along the bonds between ion pairs, proportional to the adiabatic forces b |∇ cĤ | a and to the degree of non-adiabaticity in the bond orders τ c a |v · ∇ρ 0 | b , where τ c is a local measure of the electronic scattering time at ion c. This latter term is non-local since the non-adiabaticity in bond a-b is a function of the rate of change of all the hopping integrals in the system. We can quantify the degree of non-adiabatic coupling to bond a-b with a single vector with dimensions of an inverse length: which allows us to write the non-adiabatic force on atom c as The electron scattering time appropriate for atom a must be dependent on the characteristic local electronic time scaleh D a , where D a is the local density of states at the Fermi level. It should also depend inversely on the electronic temperature: when the temperature is high there are more electronic excitations within the system to contribute to scattering. We therefore write τ a ∼h D a /Dk B T e , with D being the total density of states per atom at the Fermi level.
In appendix A, we derive a simple form forP, and argue that when the density of states is low P ab ≈ δ ab 2D a k B T e . We can further localize the force calculation by noting thatĤ ab is a function of |r ab |, so Note that the non-adiabatic force at atom a depends only on the bonds between a and its neighbours. It is highly sensitive to the electron density, scaling roughly with the second power of the density of states. Note also that the force is correctly zero for a rigid translation or rotation, in contrast to a viscous drag model. Equation (16) becomes the empirical form equation (3) under the assumptions of the second-moment tight-binding model H ab = (Ĥ ) ba and a = 1/D 2 a . An order unity scalar fitting parameter x is introduced to tune to available damping data. We can find this parameter by direct calculation, using time-dependent density functional calculations, or by fitting to the stopping and range of ions in matter (SRIM) [8] code output or experimental stopping data-for example, 8 by equating the expected energy loss rate for a single moving atom in a perfect lattice position to the loss rate from the tabulated effective viscous damping coefficient β eff : The simple example of a fast head-on collision will help illuminate the form of the model. In such collisions, the true bonding response between the collision partners lags behind the adiabatic response: the atoms are under-bonded in the first half of the collision and over-bonded in the second. If we assume an atom a collides with an initially stationary atom b, then up to the point of closest approach our model predicts a force that will act to decelerate a and accelerate b, consistent with a reduced attractive bonding interaction when compared with the instantaneous ground state. After the point of closest approach, the situation is reversed and our model predicts an accelerating force on a and a retarding force on b. These model forces agree with the picture of a retarded response of the bonding interaction.

The performance of the model
We can assess the performance of our proposed model for the non-adiabatic force by comparing its predictions with the forces derived from quantum-classical Ehrenfest simulations of idealized cases and radiation damage cascades. In all the results shown below, the model non-adiabatic force derived from the tight-binding model of copper used in the Ehrenfest simulations is exactly that which would also correspond to a Sutton-Chen copper potential [32]. We start with an idealized case, the damping of a single Einstein oscillator. Following the procedure set out in [17], we take a perfect static lattice of 2048 fcc atoms with periodic boundary conditions and an electronic temperature of 500 K. The formalism exploits a particularly simple single s-band tight-binding model metal [31] as an efficient way of combining an explicit quantum mechanical electronic system and a set of classical ions. The simplicity of this tight-binding model allows the direct simulation of collision cascades and makes clear the effects of non-adiabaticity (which might otherwise be obscured by the chemistry of a more complex model). One atom is displaced, and made to oscillate at a fixed frequency (here 0.5 rad fs −1 ), in a given direction. The damping can then be directly computed from the rise in electronic energy after each cycle, for a moving atom in different electronic environments moving in different directions. The result is shown in figure 1. The vertical scale is arbitrary energy transfer or damping units, reflecting the need for a scaling parameter. It is clear that our new model for the non-adiabatic force (16) does not just capture the coarse features of the damping, but reproduces accurately fine detail of the direction and environment dependence. Superimposed on the plot is the damping model of Caro and Victoria [9], which reproduces the environment dependence well, but simply sets the drag force direction anti-parallel to the velocity.
We have previously given details [17,18,21] of a time-dependent tight-binding formalism suitable for radiation damage simulations. We have conducted 24 cascade simulations in which the primary knock-on atom (PKA) was given 1 keV of kinetic energy in one of 24 directions evenly distributed over the 1/48th irreducible solid angle of the face-centred cubic unit cell.  The simulation cells had periodic boundaries and contained 2016 atoms (9 × 7 × 8 unit cells), all initially static at their perfect lattice sites. These initial conditions ensured that we were able to focus on the effects of non-adiabatic forces on atoms directly involved in a collision cascade rather than on thermal atoms in the surrounding material. We have simulated a relatively low PKA energy for computational expeditiousness, but we expect that our force model will remain valid whenever the single band tight-binding model is a good approximation (i.e. a similar range of validity to that of empirical potentials based on a single valence band). The model will begin to break down when excitations from core states become active, at kinetic energies of several tens of keV in the case of copper. Because a direct diagonalization of the electronic Hamiltonian is necessary in order to obtain details of the non-adiabatic force we were restricted to studying a snapshot of the forces every 0.05 fs for the first 25 fs of each cascade. However, it is in these early stages that the effects of energy exchange between electrons and ions should depend most strongly on the details of the cascade evolution and so our results provide a stricter test of the validity of our force model than if forces were averaged over a longer time frame. Figures 2-4 show the work done by the non-adiabatic forces acting on a sample of atoms from our cascade simulations and the Cartesian components of those forces. In each case, we compare the force derived from our simulations, the force predicted by the model presented in equation (16), and two drag forces anti-parallel to the ion velocity-a simple viscous drag Cartesian components of the non-adiabatic force. In each case, data are shown for the force derived from the Ehrenfest simulation, the force calculated using our new model, a simple viscous drag and the density-dependent drag of Caro and Victoria [9] (labelled C-V model). and the environment-dependent drag of Caro and Victoria [9]. A best-fit damping coefficient, discussed below, is used for each model. In figure 2, we show data for the PKA in a simulation in which it undergoes a glancing collision with another atom early in the cascade. Our force model not only captures the work done by the non-adiabatic force in the simulations, but also reproduces much of the fine detail in the direction and magnitude of the force. In figure 3, we present the same data for the second atom in a replacement collision sequence (RCS) initiated in a simulation in which the PKA is directed along the [110] close-packed direction. Again, our force model captures details of the non-adiabatic force that cannot be replicated by a simple drag, in particular, reproducing the 'double well' in the interaction of the ion with its neighbours. Finally, in figure 4, we show data for an atom set in motion later in a cascade. By this stage, many atoms are in motion and the approximation of locality made in the derivation of our force model will be less good, but we can see that the details of the non-adiabatic force in the simulations are still well reproduced. Where the model fails to capture finer details in the force (probably due to non-local effects on the density matrix), it tends to evolve as an average over the fine fluctuations.
For a statistical measure of the performance of our model, we can consider the work done on each individual atom in our simulations over the course of the 25 fs cascades. Figure 5(a) shows comparisons of our model adiabatic force, a simple drag force and the model of Caro  and Victoria [9] with the work done by the non-adiabatic force derived from the simulations. In each case, the model data have been scaled by a best-fit damping coefficient fitted using linear least-squares regression over data for all atoms for which the maximum kinetic energy in the simulations exceeded 1.0 eV. Our model significantly outperforms the drag models, showing greatly reduced scatter and achieving an R 2 goodness-of-fit measure of 0.983 compared with 0.903 for the simple drag and 0.906 for the Caro and Victoria model. In figure 5(b), we have aggregated the data of figure 5(a) so that each data-point represents the total work done by the non-adiabatic force on all the ions in each simulation. This cascade level data shows that different PKA directions can give significantly different values for this work (which is equal to the energy transferred irreversibly into the electronic system). In particular, the RCS is initiated when the PKA is directed along [110] shows a greatly enhanced damping. All these variations are well captured by our model, but are completely unaccounted for by a simple viscous damping force. The density dependence in the model of Caro and Victoria captures some of the enhancement of the damping of the RCS (as we have previously found [18]) 4 . We have previously published data from cascade simulations of longer duration (200 fs), which suggested that a simple damping force is able to capture ionic energy loss due to the non-adiabatic forces as a sum over all atoms over the full duration of a cascade [18]. Our new results show that the simple damping force is unable to capture fine details of the energy loss. A viscous damping force in MD should therefore be regarded as a mechanism for implementing the average energy loss rate and not as an accurate model for the non-adiabatic force.

Discussion and conclusions
Starting from the Ehrenfest equations of motion, we have derived a form for the non-adiabatic force in terms of the response of the density matrix to changes in the electronic Hamiltonian. We argued that the greatest contribution to the many-body terms in our expression for the evolution of the bond orders comes from simple two-body terms. Ultimately, the non-adiabatic force is determined by the lag in the response of the electronic orbitals to changes in the ionic positions. In Ehrenfest dynamics, in which the electronic system is closed in the quantum mechanical sense, the effects of different periods of the system history accumulate in the non-adiabatic part of the density matrix, but they tend to contribute incoherently and so we have been able to write a successful non-adiabatic force model in a time-local form. In an open quantum mechanical system, the historical information in the density matrix would decay with some characteristic decoherence time, but we would expect this time to be larger than the correlation time assumed in our model. The contribution to the non-adiabatic force of the 'memory' in the density matrix is limited by the former effect.
The approximations that we applied to arrive at our force model are necessary in order to obtain a local expression dependent only on the hopping integrals and other quantities accessible in any tight-binding empirical potential. Such a form is required if our model is to be implemented with minimal computational cost. Our non-adiabatic force model can be directly parameterized to correspond to any empirical potential model that includes the physics of electronic hopping integrals in its formulation-most obviously potentials of a Finnis-Sinclair form or bond-order potentials. In simulations that make use of potentials that do not take such a form, it would not be inconsistent to apply a parameterization of our model derived from a valid second-moment potential for the same metal-such a parameterization would correctly capture the form of the non-adiabatic force at the level of the approximations in our model. We have also proposed a stochastic force to model the effects of spontaneous phonon emission in returning energy from the electrons to the ions (see appendix A).
In the Lindhard picture of electronic stopping at low velocity, a moving ion induces a lagged electronic response, the linear part of which can be described using the dielectric function [14]. It is common practice to use a homogeneous, isotropic dielectric function for stopping calculations, in which case the only directionality present is given by the velocity vector of the moving ion. In our new tight-binding picture of electronic stopping, the Schrödinger equation is used to move the electrons, with non-equilibrium bond orders coming from the lagged redistribution of charge from high-energy regions to low regions energy regions. In this picture, bonds may be strengthened or weakened in front of, to the side of, or behind the moving ion, producing a directional non-adiabatic force. We take the limit where delta is small, that is to say that the non-adiabatic part of the density matrix does not accumulate coherently, but instead only appears as a gradual rise in electronic temperature [21]. In this sense, we are making the same linear approximation as is used in the dielectric theory of stopping, but with an inhomogeneous, anisotropic response, and accounting for the motion of all the ions, one pair at a time. These many-atom effects are absent from the standard dielectric stopping theory by construction (although density functional theory (DFT) linear response calculations have been performed that incorporate the effects of inhomogeneity and anisotropy in a static target lattice [22]).
We have tested our form for the non-adiabatic force against data from time-dependent tight-binding simulations of collision cascades. In these simulations, we employed a highly simplified tight-binding model in order to more cleanly separate the effects of non-adiabaticity from details of the tight-binding force. We found convincing agreement, not just of the work done by the force, but also of the detailed variation in the Cartesian components of the force acting on individual atoms throughout the early stages of collision cascades. These simulation results provide justification for our simplifying assumptions-in particular, the existence of a short correlation time τ and the spatial localization of the force model. A simple drag force is not able to capture the work done at the level of individual atoms, but it also fails to distinguish the different rates of damping of collision cascades with different PKA directions.
We have previously shown [21] that the effect of accumulating excitations on the conservative electronic forces in collision cascades could be described with an electronic temperature-dependent potential in which the evolving electronic temperature will correspond via the electronic heat capacity to the work done by the non-adiabatic forces. By implementing our proposed non-adiabatic force model in a classical MD code, it is now possible to accurately reproduce the direction and magnitude of the quantum mechanical non-adiabatic force at the atomic level. where we takeĤ 0 to mean the identity1. We can find the coefficients α n using the Cyrot-Lackmann moments theorem, Tr(Ĥ n ) = µ (n) , where µ (n) is the nth central moment of the density of states: To produce a useable expression, we need to truncate the polynomial expansion in equation (A.1), even though in doing so we inevitably lose the idempotency ofP. To first order P = α 01 + α 1Ĥ , with At low temperature α 0 → 2k B T e D(ε F ), α 1 → 2k B T e ε F D(ε F )/µ (2) , and at high temperature we have α 0 → N e , α 1 → 0, where N e is the number of orbitals per atom. To interpolate between these forms, we note that the square root of the second moment is a measure of the bandwidth of the density of states and so must be inversely proportional to the local density of states. We can therefore write an approximate interpolating form P aa = N e tanh 2D a k B T e N e , where D ab = √ D a D b is the symmetrized local density of states on atoms a and b at the Fermi level.
Note that equations (A.4) define an analytic expression for the projection operator requiring only the electronic temperature, hopping integrals and the local density of states (available through the second moment), at the sacrifice of the idempotency ofP. Until k B T e reaches the order of the bandwidth the matrix elements ofP scale linearly with temperature.
In the case of a noble metal, we can go even further in reducing the range ofP. Noting that the density of states at the Fermi level is small, we can approximate α 1 ≈ 0, and sõ P ab ≈ δ ab 2D a k B T e .

Appendix B. A stochastic return force
Spontaneous phonon emission is a purely quantum mechanical phenomenon, and so is difficult to convincingly capture in a quantum-classical picture. We might expect the true solution to couple specific electronic transitions to phonon modes-something very hard to compute correctly. However, we can make two important simplifying approximations. Firstly, we note that the time scale for a electronic de-excitation is too fast to couple to long wavelength phonon modes, and so a model that returns energy stochastically to individual atoms is sensible. Secondly, note that the total rate of energy transfer is easy to compute-it must generate an equilibrium between ions and electrons when they have matched temperatures.
Our proposed method is to add a Langevin-style stochastic return. The spatial variation of energy returned must match the variation of the strength of the electron-phonon coupling. This electron-phonon coupling strength is responsible for the environmentally dependent nonadiabatic force derived above. So we insert by hand into equation (5) a stochastic force, F → F + η, such that η = 0 and whose magnitude is chosen to produce a local thermal equilibrium.
In the Heun method algorithm, we draw a random force twice per timestep. Therefore, we can see from equation (5) that we will have an expected rise per timestep in the kinetic energy of atom a due to the stochastic force of |η a | 2 δt 2 /4m. Assuming that τ is small, the expected change in kinetic energy per timestep due to the non-adiabatic force is