# DL_POLY_4 Two temperature model

Radiation damage has long been modelled using cascade simulations in classical molecular dynamics (MD) [1], applying a high velocity to a particular atom at the start and observing the resulting defects in structure. The effects of electrons in such systems have generally been ignored, but the use of higher energy impulses requires accounting for energy lost to inelastic scattering by electrons. Not only can cold electrons quench energy from highly excited atoms as a form of electronic stopping, the electrons and atoms can also transfer energy between each other by means of electron-phonon interactions.

A two-temperature model (TTM or 2T-MD for its implementation with MD simulations) has been developed [2], which splits the system into atomistic (ionic) and electronic sections. The atomistic system is modelled using a modified form of the Langevin thermostat with a heterogeneous friction term $\chi_i$, normally set to a default value $\chi_p$ but can be increased by $\chi_s$ to represent electronic stopping of high-speed atoms, and a stochastic force to feed energy from the electrons into the atomistic system [3] that depends on local electronic temperature $T_e$.
$$m_i \frac{\partial \mathbf{v}_i}{\partial t} = \mathbf{F}_i (t) - \chi_i m_i \mathbf{v}_i + \sqrt{\frac{2 k_B T_e m_i \chi_p}{\Delta t}} \mathbf{R}_i (t)$$
The electrons are represented as a continuum gas, which is divided into a grid of voxels. Each voxel carries a local electronic temperature, which propagates by a heat equation incorporating the heat capacity $C_e$ and thermal conductivity $\kappa_e$ of the material. The heat equation includes additional source terms for transfers of energy between atoms and electrons due to temperature differences (electron-phonon coupling) and for electronic stopping: the corresponding coupling coefficients $g_p$ and $g_s$ are related to $\chi_p$ and $\chi_s$ used in the Langevin equation.
$$C_e \frac{\partial T_e}{\partial t} = \nabla \left( \kappa_e \nabla T_e \right) - g_p \left(T_e - T_a\right) + g_s T_a^{\prime}$$

Applying the 2T-MD model requires dividing the MD simulation volume into a grid of smaller voxels to calculated localised atomistic temperature. For a domain-decomposed MD code such as DL_POLY_4, it makes sense to divide the atomistic temperature voxels spatially among the available processor cores in a similar fashion to the simulation volume, as well as assign electronic temperature voxels to correspond with them and those extending beyond the atomistic system (used to control energy transfer to/from the entire system with boundary conditions). A scheme was devised to determine the cores holding voxels that cross MD domain boundaries between them, illustrated in Figure 1, which dictates the required core-to-core communications to gather information required to calculate atomistic temperatures and source terms, as well as propagate electronic temperatures with the heat equation.

The above considerations were used for the implementation of TTM in DL_POLY_4 [4], which uses an explicit finite-difference method (FDM) solver to evolve the heat equation for electronic temperatures. By applying domain decomposition to the temperature grids as well as the MD atomistic system, the parallel scalability of this implementation is in line with DL_POLY_4's, ensuring good performance for cascade simulations of hundreds of millions of atoms on thousands of processor cores. To provide a baseline for coupled TTM/MD simulations, a version of the heterogeneous Langevin thermostat (using the system temperature in place of local electronic temperature) was also implemented to implement the effects of electronic stopping alone [5].

As well as cascades, 2T-MD is also suited to simulations of solid systems subjected to laser and swift heavy ion irradiation [6], both of which can be initiated by applying energy depositions to the electronic system. In cases of thin film expansion as a result of these depositions, the temperature cells within the MD volume may become inactive due to a lack of available atoms to reliably calculate atomistic temperatures: to maintain energy conservation, electronic energy from newly inactive voxels can be transferred to active neighbours. These additional functionalities, along with the ability to supply various properties (e.g. specific heat capacities, electron-phonon coupling coefficients) in tabulated form, have also been included in DL_POLY_4's TTM implementation.

 Defects from original structure Displacements from original structure

The implementation of 2T-MD in DL_POLY_4 has been applied to high-energy cascades in, among others, iron [7], zirconia [8] and tungsten [9]. (An example of a high-energy cascade in iron can be seen in Figure 2.) It has also been applied for laser irradiation of gold [10] and tungsten films [11], as well as swift heavy ion irradiation of germanium [12], silicon [13] and various metallic structures [14].

[1] Malerba L, J. Nucl. Mater. 351 28 (2006)
[2] Duffy DM and Rutherford AM, J. Phys.: Cond. Mat. 19 016207 (2007)
[3] Caro A and Victoria M, Phys. Rev. A 40 2287 (1989)
[4] Seaton MA, Todorov IT, Daraszewicz SL, Khara GS and Duffy DM, Mol. Simul. (2018) (doi: 10.1080/08927022.2018.1511902)
[5] Zarkadoula E, Daraszewicz SL, Duffy DM, Seaton MA, Todorov IT, Nordlund K, Dove MT and Trachenko K, J. Phys.: Cond. Mat. 25 125402 (2013)
[6] Darkins R and Duffy DM, Comput. Mater. Sci. 147 145 (2018)
[7] Zarkadoula E, Daraszewicz SL, Duffy DM, Seaton MA, Todorov IT, Nordlund K, Dove MT and Trachenko K, J. Phys.: Cond. Mat. 26 085401 (2014)
[8] Zarkadoula E, Devanathan R, Weber WJ, Seaton MA, Todorov IT, Nordlund K, Dove MT and Trachenko K, J. Appl. Phys. 115 083507 (2014)
[9] Zarkadoula E, Duffy DM, Nordlund K, Seaton MA, Todorov IT, Weber WJ and Trachenko K, J. Phys.: Cond. Mat. 27 135401 (2015)
[10] Daraszewicz SL, Giret Y, Naruse N, Murooka Y, Yang J, Duffy DM, Shluger AL and Tanimura K, Phys. Rev. B 88 184101 (2013)
[11] Murphy ST, Daraszewicz SL, Giret Y, Watkins M, Shluger AL, Tanimura K and Duffy DM, Phys. Rev. B 92 134110 (2015)
[12] Daraszewicz SL and Duffy DM, Nucl. Instrum. Methods Phys. Res. B 303 112 (2013)
[13] Khara GS, Murphy ST, Daraszewicz SL and Duffy DM, J. Phys.: Cond. Mat. 26 395201 (2016)
[14] Khara GS, Murphy ST and Duffy DM, J. Phys.: Cond. Mat. 29 285303 (2017)