On-the-fly meanfield transition-state theory for diffusive molecular dynamics

Contents lists available at ScienceDirect Mechanics of Materials journal homepage: www.elsevier.com/locate/mecmat Research paper On-the-fly meanfield transition-state theory for diffusive molecular dynamics M. Molinosa, M. Ortizb,c, M.P. Arizaa,∗ aEscuela Técnica Superior de Ingeniería, Universidad de Sevilla, Camino de los descubrimientos, s.n., 41092, Sevilla, Spain bDivision of Engineering and Applied Science, California Institute of Technology, 1200 E. California Blvd., Pasadena, CA 91125, USA cCentre Internacional de Métodes Numerics en Enginyeria (CIMNE), Universitat Politècnica de Catalunya, Jordi Girona 1, Barcelona, 08034, Spain A R T I C L E I N F O Keywords: Magnesium Magnesium hydrides Angular-dependent interatomic potentials Mean field approximation Transition state theory Mass transport A B S T R A C T We apply transition state theory to derive atomic-level master equations for mass transport from empirical interatomic potentials within the Diffusive Molecular Dynamics (DMD) framework. We show that meanfield approximation provides an exceedingly efficient and accurate means of computing free-energy barriers in arbitrary local atomic configurations, thus enabling long-term DMD ‘on-the-fly’ and on the sole basis of an underlying interatomic potential, without additional modeling assumptions. We apply and validate the resulting meanfield DMD paradigm in simulations of processes of hydrogenation and dehydrogenation of Mg using Angular-Dependent interatomic Potentials (ADP). We show that meanfield DMD correctly predicts hydrogen diffusivities in hcp Mg and vacancy diffusivities in rutile MgH2. We demonstrate the ability of meanfield DMD to predict evolution through calculations concerned with dilute concentrations of hydrogen in hcp Mg, and with dilute concentrations of hydrogen vacancies in rutile MgH2, including off-stoichiometry hydrogen concentrations and temperature effects. Remarkably, the time steps required by DMD are up to six orders of magnitude larger than those required by Molecular Dynamics (MD), which demonstrates the overwhelming superiority of the DMD paradigm in simulations of phenomena occurring on the diffusive time scale. 1. Introduction Metal hydrides are attractive for vehicular hydrogen (Partnership, 2017), hydrogen storage (Mohtadi and Orimo, 2016; Yang et al., 2021), and many other emerging applications, and remain the subject of extensive ongoing research (Yang et al., 2021; Tan and Ramakrishna, 2021). Hydrogenation and dehydrogenation processes in metals are rate limited by hydrogen uptake and mass transport within the host metal and, therefore, operate on a diffusive time scale (Luo et al., 2019; Shriniwasan and Tatiparti, 2019; Pundt, 2004; Shen and Aguey-Zinsou, 2016). Hydrogen diffusivity in metals is a complex phenomenon due to differences in crystal structure between the host metal and its hydride phases, which inevitably induces microstructural evolution governed by transformation kinetics. Magnesium hydride, the main focus of the present study, is a prime example of a transformational hydrogen storage system: hexagonal closely packed (hcp) structure for Mg and tetragonal rutile structure for 𝛼-MgH2. Owing to these complexities, hydrogen diffusivity cannot be easily inferred from experimental data (Burger et al., 1961; Arons et al., 1970, 1974; Cornell and Seymour, 1975; Mazzolai and Zuchner, 1981; Nishimura et al., 1999a; Fernandez and Sanchez, 2002; Cermak and Kral, 2008; Li et al., 2018). A predictive understanding of the mechanisms underlying hydrogenation and dehydrogenation processes in metals therefore requires ∗Corresponding author. E-mail addresses: [email protected] (M. Molinos), [email protected] (M. Ortiz), [email protected] (M.P. Ariza). methods of analysis that deliver atomistic fidelity and, simultaneously, the ability to reach across to the diffusive time scale. At the most basic level, the atomistic models must accurately predict equilibrium properties and transition energies for hydrogen transport within the host metal lattice. Density Functional Theory (DFT) (Tao et al., 2009; Klyukin et al., 2015), ab initio methods (Klyukin et al., 2013; Junkaew et al., 2014) and empirical potentials (Mishin and Lozovoi, 2006; Smirnova et al., 2018; Zhou et al., 2019a; Molinos et al., 2024) have been widely used to that end. Molecular dynamics (MD) has also been used to predict hydrogen diffusivities in metals and their hydrides (Zhou et al., 2016, 2017, 2018b,?; Spataru et al., 2020). However, MD is severely limited by the need to resolve the thermal vibrations of the atoms, with the result that the calculations encompass times that are exceedingly short relative to the diffusive time scale of interest. This limitation of MD has often been sidestepped by recourse to kinetic Monte Carlo (KMC) methods for purposes of simulating transport phenomena mediated by hopping transitions and other rare events (Bortz et al., 1975; Voter, 2005; Battaile, 2008; Martinez et al., 2011; Reina et al., 2011). However, KMC methods require the a priori enumeration of all possible transition paths of the system and the elucidation of the corresponding transition rates, https://doi.org/10.1016/j.mechmat.2025.105380 Received 31 January 2025; Received in revised form 29 April 2025; Accepted 29 April 2025 Mechanics of Materials 207 (2025) 105380 Available online 19 May 2025 0167-6636/© 2025 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY-NC license ( http://creativecommons.org/licenses/bync/4.0/ ). M. Molinos et al. which renders the approach intractable when the transition paths are numerous and complex, as expected in transformational systems. The same curse of complexity besets efforts to learn transition paths and energy barriers using machine learning methods (Tang et al., 2024). Therefore, there remains a need for modeling paradigms capable of delivering full atomistic fidelity while simultaneously bridging the vibrational and diffusive time scales. One such paradigm, proposed by Kulkarni (2006), Kulkarni et al. (2008), Venturini (2011), Venturini et al. (2014) and termed Diffusive Molecular Dynamics (DMD) by Li et al. (2011), is predicated on a representation of the state of the system undergoing mass transport as a collection of sites that can be occupied or empty. The energetics of the system is described by means of conventional empirical potentials and a non-equilibrium statistical– mechanical treatment of the ensemble allows the sites to be partially occupied and accounts for thermal effects, possibly including heat transport. The evolution of the state is governed by kinetic equations in the spirit of Onsager kinetics, which can be viewed as atomic-level Fick and Fourier laws. The methodology has been successfully applied to problems of heat transport (Kulkarni, 2006; Kulkarni et al., 2008; Ariza et al., 2011; Ponga et al., 2015, 2016; Gupta et al., 2021) and mass transport (Venturini, 2011; Li et al., 2011; Sarkar et al., 2012; Venturini et al., 2014; Martin et al., 2015; Wang et al., 2015; Simpson et al., 2016; Sun et al., 2017; Mendez et al., 2018; Mendez and Ponga, 2021). The main computational challenges inherent to DMD are: the computation of phase-space integrals required to evaluate local equilibrium relations; and the formulation of effective kinetic laws from empirical interatomic potentials. In previous work, including applications to Mg-H systems (Molinos et al., 2024), we have shown that the first challenge, the evaluation of phase-space-integrals, can be addressed accurately and efficiently by recourse to meanfield approximation and numerical quadrature. In combination with Angular-Dependent interatomic Potentials (ADP) (Smirnova et al., 2018; Mishin and Lozovoi, 2006) the resulting thermalized and mixed ADP potentials accurately predict equilibrium properties of Mg and its hydrides, including free entropy, heat capacity, thermal expansion, molar volumes, equation of state and elastic constants (Molinos et al., 2024). In the present work, we appeal to transition state theory (Weiner, 2012) to derive atomic-level master equations for mass transport from empirical interatomic potentials, allowing for arbitrary reconfigurations of the host metal lattice, within the context of the DMD framework of Venturini (2011). The fundamental question to be elucidated concerns the efficient characterization of attempt frequencies and energy barriers attendant to hydrogen transitions occurring in arbitrary – possibly complex – local atomic configurations, such as free surfaces, phase boundaries, grain boundaries, amorphized regions and others. The nudged elastic band method (Henkelman and Jonsson, 2000; Henkelman et al., 2000; Nakano, 2008) is widely used to compute transition energy barriers at 0K. However, the calculations are challenging and costly for complex high-dimensional energy landscapes. In addition, for systems undergoing displacive phase transitions, or generally in the vicinity of extended lattice defects, the number and complexity of possible transition paths is exceedingly large and difficult – if not impossible – to parameterize a priori. By contrast, we show that meanfield approximation, combined with the occupancy-variable representation, vastly simplifies the implementation of transition-state theory and supplies an exceedingly efficient and accurate means of computing free-energy barriers in complex environments, thus enabling DMD ‘on-the-fly’ and on the sole basis of an underlying interatomic potential, without additional modeling assumptions. By way of demonstration, we apply and validate the resulting meanfield DMD methodology to the simulation of processes of hydrogenation and dehydrogenation of Mg using ADP potentials (Smirnova et al., 2018; Molinos et al., 2024). We validate predictions of hydrogen diffusivity in hcp Mg against experimentally measured averaged bulk diffusivities (Nishimura et al., 1999b; Renner and Grabke, 1978) and DFT calculations (Klyukin et al., 2015). The predicted diffusivities are in good agreement with the experimental data and correctly capture the experimentally observed strong anisotropy and hexagonal symmetry. We additionally validate predictions of vacancy diffusivity in rutile MgH2. Here again, the experimentally observed strong anisotropy of the diffusivities is captured by the theory. The predicted hydrogen diffusivities in magnesium hydride are much lower than in magnesium, also in agreement with experimental data and the MD calculations of Spataru et al. (2020). Finally, in order to showcase and assess the ability and efficiency of meanfield DMD to predict evolution, we present calculations concerned with dilute concentrations of hydrogen in hcp Mg, and with dilute concentrations of hydrogen vacancies in magnesium hydride, including off-stoichiometry hydrogen concentrations and temperature effects. Both processes are encountered in practice during the operation cycle of hydrogen storage materials. In all cases, the computational setup is replicated from Spataru et al. (2020) for purposes of comparison with MD. We find that meanfield DMD accurately captures both the short-term and long-term response of the systems, including a transition from classical to ballistic diffusion, lattice distortions resulting from the hydrogen transport, and lattice stability as a function of temperature. Remarkably, in the case of vacancy diffusion in rutile MgH2 the time steps required by DMD are six orders of magnitude larger that those required by MD, which represents a staggering acceleration and demonstrates the overwhelming superiority of the DMD paradigm for simulating phenomena occurring on the diffusive time scale. 2. Local equilibrium relations and meanfield approximation By way of background and to set the framework, we summarize salient aspects of the non-equilibrium statistical mechanics theory of Venturini et al. (2014) and its application to Mg-H systems (Molinos et al., 2024). 2.1. Non-equilibrium statistical mechanics We consider a closed system consisting of 𝑁 sites, e. g., atomic positions or molecules, each of which can be of one of 𝑀 species, including vacancies. For each site 𝑖= 1,…, 𝑁, and each species 𝑘= 1,…, 𝑀, we introduce the occupancy function 𝑛𝑖𝑘 ={1,if site 𝑖 is occupied by species 𝑘, 0,otherwise,(1) in order to describe the occupancy of each site. We note that, from definition (1), we must have 𝑀 ∑ 𝑘=1 𝑛𝑖𝑘 = 1 (2) at every site 𝑖. We additionally denote by 𝑛𝑖= (𝑛𝑖𝑘)𝑀 𝑘=1 the local occupancy array of site 𝑖. The microscopic states of the system are defined by the instantaneous position {𝑞} = (𝑞𝑖)𝑁 𝑖=1, momenta {𝑝} = (𝑝𝑖)𝑁 𝑖=1, and occupancy arrays {𝑛} = (𝑛𝑖)𝑁 𝑖=1 of all 𝑁 sites in the system. The occupancy functions 𝑛𝑖 take values in a set 𝑀 consisting of the elements of {0,1}𝑀 that satisfy the constraint (2). In addition, the occupancy arrays {𝑛} take values in the set 𝑁𝑀 ={{𝑛} ∈ {0,1}𝑁𝑀 ∶ 𝑛𝑖∈𝑀for 𝑖= 1,…, 𝑁 }. (3) The expected or macroscopic value of a function 𝐴({𝑞}, {𝑝},{𝑛}) is given by the phase average ⟨𝐴⟩=∑ {𝑛}∈𝑁𝑀 1 ℎ3𝑁× ∫𝛤 𝐴({𝑞},{𝑝},{𝑛}) 𝜌({𝑞},{𝑝},{𝑛}) 𝑑𝑞 𝑑𝑝, (4) Mechanics of Materials 207 (2025) 105380 2 M. Molinos et al. where 𝛤= (R3×R3)𝑁, ℎ is Planck’s constant and ℎ−3𝑁 supplies the natural unit of phase volume for systems of distinguishable particles (Hill, 1987; Girifalco, 2000). We assume that the statistics of the system obeys Jaynes’ principle of maximum entropy (Jaynes, 1957a,b; Zubarev, 1974; Callen, 1985), which posits that the probability density function 𝜌({𝑞},{𝑝},{𝑛}), characterizing the probability of finding the system in a state ({𝑞},{𝑝},{𝑛}), maximizes the information-theoretical entropy [𝜌] = −𝑘𝐵⟨log 𝜌⟩,(5) among all probability measures consistent with the constraints on the system. In (5) and subsequently, 𝑘𝐵 denotes Boltzmann’s constant. We consider systems consisting of distinguishable particles whose Hamiltonians have the additive structure 𝐻({𝑞},{𝑝},{𝑛}) = 𝑁 ∑ 𝑖=1 ℎ𝑖({𝑞},{𝑝},{𝑛}),(6) where ℎ𝑖 is the local Hamiltonian of particle 𝑖. Following Feynman (1998), suppose that the expected particle energies and the expected particle atomic fractions ⟨ℎ𝑖⟩=𝑒𝑖,⟨𝑛𝑖𝑘⟩=𝜒𝑖𝑘,(7) are known, respectively. We note that the local atomic fractions satisfy the identities 𝑀 ∑ 𝑘=1 𝜒𝑖𝑘 = 1,(8) which follow from (2) and the second of (7). Enforcing the constraints (7) by means of Lagrange multipliers 𝑘𝐵{𝛽}≡(𝑘𝐵𝛽𝑖)𝑁 𝑖=1 and {𝛾}≡ ((𝛾𝑖𝑘)𝑀 𝑘=1)𝑁 𝑖=1, leads to the Lagrangian [𝜌, {𝛽},{𝛾}] = [𝜌] − 𝑘𝐵{𝛽}𝑇{⟨ℎ⟩} + 𝑘𝐵{𝛾}𝑇{⟨𝑛⟩},(9) where we interpret 𝑇𝑖=1 𝑘𝐵𝛽𝑖 ,(10) as the local temperature of particle 𝑖. In view of identities (8), we additionally append the constraints 𝑀 ∑ 𝑘=1 𝛾𝑖𝑘 = 0,(11) in order to render {𝛾} determinate. Maximizing [⋅,{𝛽},{𝛾}] among probability measures gives 𝜌({𝑞},{𝑝},{𝑛}) = 1 𝛯e−{𝛽}𝑇{ℎ}+{𝛾}𝑇{𝑛},(12) with 𝛯=∑ {𝑛}∈𝑁𝑀 1 ℎ3𝑁∫𝛤 e−{𝛽}𝑇{ℎ}+{𝛾}𝑇{𝑛}𝑑𝑞 𝑑𝑝. (13) The corresponding equilibrium values of {𝛽} and {𝛾} follow from (7) as a function of {𝑒} and {𝜒}. We interpret (12) and (13) as non-equilibrium generalizations of the Gibbs grand-canonical probability density function and the grand-canonical partition function, respectively. 2.2. Meanfield approximation The calculation of the thermodynamic potentials in closed form is generally intractable and approximation is therefore required. A variational meanfield theory (Stanley, 1971; Yeomans, 1992a) may be formulated by restricting (9) to some class of probability density functions of the form 𝜌0({𝑞},{𝑝},{𝑛}) = 1 𝛯0 𝑒−{𝛽}𝑇{ ℎ}+{𝛾}𝑇{𝑛},(14) with 𝛯0=∑ {𝑛}∈𝑁𝑀 1 ℎ3𝑁∫𝛤 𝑒−{𝛽}𝑇{ ℎ}+ {𝛾}𝑇{𝑛}𝑑𝑞 𝑑𝑝. (15) and { ℎ} in some class 0 of local trial Hamiltonians, possibly defined parametrically. The restricted Lagrangian is [𝜌0,{𝛽},{𝛾}] = [𝜌0]− 𝑘B{𝛽}𝑇{⟨ℎ⟩0} + 𝑘B{𝛾}𝑇{⟨𝑛⟩0},(16) where [𝜌0] = −𝑘𝐵⟨log 𝜌0⟩0(17) and ⟨⋅⟩0 denotes averaging with respect to 𝜌0. Inserting (14) and (15) into (16) using (17) gives [𝜌0,{𝛽},{𝛾}] = 𝑘Blog 𝛯0−𝑘B{𝛽}𝑇{⟨ℎ− ℎ⟩0}.(18) The optimal trial Hamiltonians are determined by maximizing [𝜌0,{𝛽},{𝛾}] with respect to 𝜌0, or some suitable parametrization of { ℎ} thereof. In addition, the corresponding meanfield equilibrium values of {𝛽} and {𝛾} follow as a function of {𝑒} and {𝜒} from the Euler–Lagrange equations of [𝜌0,{𝛽},{𝛾}], ⟨ ℎ𝑖⟩0=𝑒𝑖,⟨𝑛𝑖𝑘⟩0=𝜒𝑖𝑘.(19) 2.3. Interstitial hydrogen in metals As a special case of the general theory just outlined, consider a crystal lattice where the base lattice sites are always occupied by metal atoms, while the interstitial sites are either occupied by an H atom, or unoccupied (Sun et al., 2017). We index the base lattice sites by an index set 𝐼M, and the interstitial sites by 𝐼H. Under these assumptions, we can characterize occupancy by a single occupancy number 𝑛𝑖 on the interstitial sites, 𝑖∈𝐼H, taking the value of 0 if the site unoccupied and 1 if the site occupied. We assume a Hamiltonian of the form 𝐻({𝑞},{𝑝},{𝑛}) = ∑ 𝑖∈𝐼M 1 2𝑚M|𝑝𝑖|2+∑ 𝑖∈𝐼H 1 2𝑚H|𝑝𝑖|2+𝑉({𝑞},{𝑛}),(20) where 𝑚M and 𝑚H are the atomic masses of the metal host and hydrogen, respectively, and 𝑉({𝑞},{𝑛}) is a mixed angular-dependent interatomic potential (ADP) of the form 𝑉({𝑞},{𝑛}) = ∑ 𝑖 𝑛𝑖(𝜌𝑖) + 1 2∑ 𝑖∑ 𝑗∈𝑖(𝑟c) 𝑗≠𝑖 𝑛𝑖𝑛𝑗𝜙𝑖𝑗 +1 2∑ 𝑖∑ 𝑗1,𝑗2∈𝑖(𝑟c) 𝑗1,𝑗2≠𝑖(𝑛𝑖𝑛𝑗1𝜇𝑖𝑗1)⋅(𝑛𝑖𝑛𝑗2𝜇𝑖𝑗2), +1 2∑ 𝑖∑ 𝑗1,𝑗2∈𝑖(𝑟c) 𝑗1,𝑗2≠𝑖(𝑛𝑖𝑛𝑗1𝜆𝑖𝑗1)∶(𝑛𝑖𝑛𝑗2𝜆𝑖𝑗2) −1 6∑ 𝑖∑ 𝑗1,𝑗2∈𝑖(𝑟c) 𝑗1,𝑗2≠𝑖(𝑛𝑖𝑛𝑗1tr𝜆𝑖𝑗1)(𝑛𝑖𝑛𝑗2tr𝜆𝑖𝑗2), (21) which belongs to a class of many-body potentials proposed by Mishin and Lozovoi (2006). In (21), the sub-index 𝑖 denotes each atomic site of the domain while sub-indexes 𝑗, 𝑗1 or 𝑗2 are used to enumerate the neighbors of each site 𝑖 and 𝑖(𝑟c) denotes the local neighborhood of 𝑖 with cutoff radius 𝑟c. The first term is the embedding energy, the second term is the pairwise-interaction energy (𝑖−𝑗), the third to fifth terms are the contributions to the energy due to angular interactions (𝑖−𝑗1−𝑗2) where 𝜇𝑖𝑗 =𝑢𝑖𝑗 𝐫𝑖𝑗 and 𝜆𝑖𝑗 =𝑤𝑖𝑗 𝐫𝑖𝑗 ⊗𝐫𝑖𝑗 (22) are respectively the dipole and quadrupole. Where 𝐫𝑖𝑗 is the distance vector between a particle 𝑖 and its neighbor 𝑗. The scalar functions Mechanics of Materials 207 (2025) 105380 3 M. Molinos et al. (𝜌𝑖), 𝜙(𝑟𝑖𝑗 ), 𝑢(𝑟𝑖𝑗 ) and 𝑤(𝑟𝑖𝑗 ) are spline functions commonly represented in tabular form e. g. Smirnova et al. (2018). These functions are expressed in terms of the distance norm between two sites, |𝐫𝑖𝑗 |=𝑟𝑖𝑗 , and the electron density 𝜌𝑖. See also Molinos et al. (2024) for further details of the numerical treatment and for a validation assessment of the potential. The potential energy (21) is a function of the hydrogen occupancies and, in that sense, may be regarded as a three-dimensional Ising model. Such models have been extensively studied by a variety of means, including meanfield theory (cf., e. g., Yeomans (1992b)). Building on that background, we assume that the system is closed and in thermal equilibrium, whence 𝛽𝑖=𝛽= 1∕𝑘B𝑇 for all particles, and choose trial Hamiltonians of the form  ℎ𝑖(𝑞𝑖, 𝑝𝑖) = 𝑘B𝑇 2𝜎2 𝑖|𝑞𝑖−𝑞𝑖|2+1 2𝑚𝑖|𝑝𝑖−𝑝𝑖|2,(23) where 𝑞𝑖, 𝑝𝑖 and 𝜎𝑖 are parameters that characterize the trial space. The corresponding meanfield probability density function (14) evaluates to 𝜌0({𝑞},{𝑝},{𝑛}) = 1 𝛯0 exp { −∑ 𝑖∈𝐼M∪𝐼H 1 2𝜎2 𝑖|𝑞𝑖−𝑞𝑖|2−∑ 𝑖∈𝐼M 𝛽 2𝑚M|𝑝𝑖−𝑝𝑖|2 −∑ 𝑖∈𝐼H(𝛽 2𝑚H|𝑝𝑖−𝑝𝑖|2−𝛾𝑖𝑛𝑖)}, (24) and the meanfield grand-canonical partition function (15) to 𝛯0={∏ 𝑖∈𝐼M(𝜎𝑖√𝑚M∕𝛽 ℏ)3}× {∏ 𝑖∈𝐼H(𝜎𝑖√𝑚H∕𝛽 ℏ)3(1+e𝛾𝑖)}, (25) with ℏ the reduced Planck constant. Under the assumed isothermal conditions, the meanfield Lagrangian (16) reduces to [𝜌0,{𝛽},{𝛾}] = 𝑘Blog 𝛯0− 𝑘𝐵𝛽(⟨𝑉⟩0+∑ 𝑖∈𝐼M∪𝐼H 1 2𝑚𝑖|𝑝𝑖|2)+ 3 𝑁 𝑘B,(26) where 𝑁= #𝐼M+ #𝐼H is the total number of sites in the system. The corresponding meanfield optimality conditions are −𝜕 𝜕 𝑞𝑖 =𝜕 𝜕 𝑞𝑖⟨𝑉⟩0=⟨𝜕𝑉 𝜕𝑞𝑖⟩0= 0,(27a) −𝜕 𝜕 𝑝𝑖 =1 𝑚𝑖 𝑝𝑖= 0,(27b) −𝜕 𝜕 𝜎𝑖 = − 3𝑘B 𝜎𝑖 +𝑘B𝛽𝜕 𝜕 𝜎𝑖⟨𝑉⟩0= 0,(27c) and the meanfield Euler–Lagrange Eqs. (19) evaluate to ⟨𝑉𝑖⟩0+|𝑝𝑖|2 2𝑚𝑖 =𝑒𝑖, 𝑖 ∈𝐼M∪𝐼H,(28a) e𝛾𝑖 1+e𝛾𝑖=𝜒𝑖, 𝑖 ∈𝐼H.(28b) Alternatively, Eq. (28b) can be inverted to yield 𝛾𝑖= log(𝜒𝑖 1 − 𝜒𝑖), 𝑖 ∈𝐼H.(29) The equilibrium relation (29) is plotted in Fig. 1 for a uniform H atomic fraction in a perfect hcp lattice and a range of temperatures. We recall from (4) that ⟨𝑉⟩0=∑ {𝑛}∈𝑁𝑀 1 ℎ3𝑁× ∫𝛤 𝑉({𝑞},{𝑛})𝜌0({𝑞},{𝑝},{𝑛}) 𝑑𝑞 𝑑𝑝. (30) Fig. 1. Local equilibrium relation (29) expressed in terms of the local chemical potential 𝜇𝑖∶= 𝛾𝑖∕𝛽𝑖 and the local hydrogen molar fraction 𝜒𝑖 of H at different temperatures. In view of (24) and with reference to Jensen’s inequality, we approximate (30) as ⟨𝑉⟩0≈∫𝑉({𝑞},{𝜒}) {∏ 𝑗∈𝐼M∪𝐼H 1 (√2𝜋 𝜎𝑗)3 exp(−1 2𝜎2 𝑗|𝑞𝑗−𝑞𝑗|2)}𝑑𝑞. (31) as proposed by Venturini et al. (2014) to avoid occupancy sums of combinatorial complexity. In addition, we approximate the remaining integral over configuration space by means of numerical quadrature, see Molinos et al. (2024) for details of the numerical implementation and verification thereof. 3. Configuration-dependent hydrogen transport kinetics In order to formulate a general framework for kinetics, including heat and mass transport, we begin by examining the balance of energy at the particle level (Venturini et al., 2014). The internal energy 𝑒𝑖 of particle 𝑖 may be identified with the expected value of the corresponding particle Hamiltonian as in the first of (7). Suppose, in addition, that the local Hamiltonians ℎ𝑖 depend on macroscopic variables {} = (1,…,𝜈), such as volume or deformation (Weiner, 2002). Then, the balance of energy at particle 𝑖 can be expressed as 𝑒𝑖=𝜇𝑇 𝑖𝜒𝑖+ 𝑤𝑖+𝑟𝑖, 𝑤𝑖= 𝜈 ∑ 𝛼=1⟨𝜕ℎ𝑖 𝜕𝛼  𝛼⟩(32) where 𝜇𝑖=𝛾𝑖∕𝛽𝑖(33) is particle chemical potential,  𝑤𝑖 is the external mechanical power and 𝑟𝑖 is the heat flow into particle 𝑖. We assume that the heat flow 𝑟𝑖 into particle 𝑖 is of the form 𝑟𝑖=∑ 𝑗≠𝑖 𝑅𝑖𝑗 , 𝑅𝑖𝑗 = −𝑅𝑗𝑖,(34) where the sum extends to all particles different from 𝑖, and 𝑅𝑖𝑗 is the discrete heat flux from particle 𝑗 to 𝑖. Likewise, we shall assume that the mass flow into particle 𝑖 may be expressed as 𝑥𝑖=∑ 𝑗≠𝑖 𝐽𝑖𝑗 , 𝐽𝑖𝑗 = −𝐽𝑗𝑖,(35) Mechanics of Materials 207 (2025) 105380 4 M. Molinos et al. where 𝐽𝑖𝑗 is the discrete mass flux array from particle 𝑗 to 𝑖. The internal entropy production rate for a particle pair is given by 𝛴𝑖𝑗 =𝐾𝑇 𝑖𝑗 𝐽𝑖𝑗 +𝑃𝑖𝑗 𝑅𝑖𝑗 ≥0,(36) where we write 𝐾𝑖𝑗 =𝜇𝑖 𝑇𝑖 −𝜇𝑗 𝑇𝑗 =𝑘𝐵(𝛾𝑖−𝛾𝑗), 𝑃𝑖𝑗 =1 𝑇𝑖 −1 𝑇𝑗 =𝑘𝐵(𝛽𝑖−𝛽𝑗). (37) Following Onsager (1931a,b), the local dissipation inequality (36) suggests kinetic laws of the general form 𝐽𝑖𝑗 = − 𝜕𝛹 𝜕𝐾𝑖𝑗 ({𝑃},{𝐾}), 𝑅𝑖𝑗 =𝜕𝛹 𝜕𝑃𝑖𝑗 ({𝑃},{𝐾}), (38) where 𝛹({𝑃},{𝐾}) is a discrete kinetic potential. The ability of models based on non-equilibrium statistical mechanics to reproduce the observed anisotropy, temperature and size dependence of the thermal conductivity of Si nanowires was established by Martin et al. (2015) by way of validation of the theory. The models have also demonstrated predictive ability in applications including nanovoid growth in metals at low and high strain rates (Ariza et al., 2011; Ponga et al., 2017). 3.1. Mass transport Next, we specialize the general theory to mass transport. To this end, we consider a transition from an initial equilibrium state 𝑠∶= ({𝑞},{𝑛}) to a final equilibrium state 𝑠′∶= ({𝑞′},{𝑛′}). The two states differ in the occupancy numbers, some of which may have flipped, and the positions of the atoms. The transition 𝑠′→𝑠 is assumed to follow a continuous transition path 𝑠(𝜉), parameterized by 𝜉∈ [0,1], with 𝑠(0) = 𝑠′ and 𝑠(1) = 𝑠. For the energy of the system to evolve continuously throughout the transition, pairs of sites 𝑖 and 𝑗, one initially occupied, 𝑛′ 𝑖= 1, and one initially unoccupied, 𝑛′ 𝑗= 0, must exchange mass, 𝑛′ 𝑖→𝑛𝑖= 0 and 𝑛′ 𝑗→𝑛𝑗= 1, at the point 𝜉𝑐 of their trajectory when their sites are coincident, 𝑞𝑖(𝜉𝑐) = 𝑞𝑗(𝜉𝑐). According to transition state theory (Weiner, 2012), the average flipping rate is given by the Arrhenius relation 𝑓𝑖→𝑗=𝜈𝑖e−𝛽𝐸𝑖→𝑗.(39) where 𝜈𝑖 is the attempt frequency of site 𝑖 and 𝐸𝑖→𝑗 is the energy barrier for 𝑖→𝑗 hops. Further interpreting the atomic fractions 𝑥𝑖 as occupancy probabilities, the probability that site 𝑖 be occupied is precisely 𝑥𝑖 and the probability that the site 𝑗 be unoccupied is (1−𝑥𝑗), giving a transition probability 𝜓𝑖→𝑗=𝑥𝑖(1 − 𝑥𝑗)𝑓𝑖→𝑗.(40) The atomic fraction rate is then given by the master equation (Nordsieck et al., 1940) 𝑥𝑖=∑ 𝑗≠𝑖(𝜓𝑗→𝑖−𝜓𝑖→𝑗).(41) We verify from (41) that ∑ 𝑖∈𝐼𝐻 𝑥𝑖= 0,(42) i. e., the total mass of the system is conserved in the absence of sources and sinks, as required for a chemically isolated system. 3.2. Calculation of energy barriers We see from (39) that, within the present framework, the formulation of kinetic models of mass transport requires the specification of suitable forms of 𝜈𝑖 and 𝐸𝑖→𝑗 and, specifically, of their dependence on the atomic configuration, e. g., through the local environment. The calculation of attempt frequencies and energy barriers is the main focus of transition-state theory (Weiner, 2002) and has traditional been a main computational bottleneck in practice. Vineyard’s formula (Vineyard, 1957), which is derived from a statistical–mechanical analysis of harmonic approximations of the potential energy about energy wells and saddle points, and numerous extensions and enhancements thereof, are widely used in practice. We show next that the use of occupancy variables and meanfield approximation vastly simplifies the implementation of transition-state theory and enables the on-the-fly evaluation of energy barriers and attempt frequencies, thus rendering the approach practical. 3.2.1. Hessian algorithm and Vineyard’s formula We begin by specializing Vineyard’s formula (Vineyard, 1957) to systems of particles described by means of occupancy variables. We confine attention to transition paths such that the total energy 𝐸({𝑞(𝜉)},{𝑛(𝜉)}), parameterized by an order parameter 𝜉∈ [0,1], attains local minima at the initial and final states and the minima are separated by one single maximum, or energy barrier. Among such transition paths, we seek to determine that for which the energy barrier is smallest. By this condition, the energy maximum necessarily occurs at a saddle point, ({𝑞(𝜉𝑐)},{𝑛(𝜉𝑐)}), or transition state. In order to find the position of the saddle point between two local minima 𝑖 and 𝑗, we look for the maximum between 𝑖 and 𝑗 using the search direction given by the local gradient and Hessian of the total energy potential. Once a local maximum point is found, we test whether it is indeed a saddle point of the potential. The sought energy barrier is then 𝐸𝑏=𝐸({𝑞(𝜉𝑐)},{𝑛(𝜉𝑐)}) − 𝐸({𝑞(0)},{𝑛(0)}),(43) and the attempt frequency of site 𝑖 is given by Vineyard’s formula (Vineyard, 1957) 𝜈𝑖=1 2𝜋 det[𝐷2𝐸({𝑞(0)},{𝑛(0)})] det[𝐷2𝐸({𝑞(𝜉𝑐)},{𝑛(𝜉𝑐)})+].(44) In this expression, 𝐷2𝐸({𝑞},{𝑛}) is the matrix of second derivatives, or Hessian, of 𝐸 at ({𝑞},{𝑛}), and the subscript ()+ designates its positive component, i. e., the component obtained by considering positive eigenvalues only. 3.2.2. Meanfield calculation of energy barriers Despite the appeal of Vineyard’s formula, and extensions thereof, its main drawback is that the determination of transition paths between local energy minima is exceedingly difficult in general and computationally expensive for large systems. We overcome this difficulty by instead estimating the transitions by recourse to variational meanfield theory. By consistency with the meanfield probability density (24), we approximate the total internal energy as 𝐸({𝑞},{𝑛}) = ∑ 𝑖∈𝐼 𝑘B𝑇 2𝜎2 𝑖|𝑞𝑖−𝑞𝑖|2−∑ 𝑖∈𝐼𝐻 𝑘B𝑇 𝛾𝑖𝑛𝑖,(45) where {𝑞} and {𝜎} are the meanfield variables introduced in Section 2.2. We then consider transitions between pairs of hydrogen sites 𝑖, 𝑗∈𝐼𝐻, 𝑖≠𝑗, and assume that the remaining occupancies 𝜒𝑘 remain unchanged, and that the remaining particle positions 𝑞𝑘 remain close to their average values 𝑞𝑘, 𝑘∈𝐼𝐻, 𝑖≠𝑘≠𝑗, through the transition. Under these assumptions, from (45) we have 𝐸(𝜉)∼ 𝜅𝑖 2|𝑞𝑖(𝜉) − 𝑞𝑖|2−𝜑𝑖,for 𝜉∼ 0,(46a) Mechanics of Materials 207 (2025) 105380 5 M. Molinos et al. 𝐸(𝜉)∼ 𝜅𝑗 2|𝑞𝑗(𝜉) − 𝑞𝑗|2−𝜑𝑗,for 𝜉∼ 1,(46b) where we write 𝜅𝑖=𝑘B𝑇 𝜎2 𝑖 , 𝜑𝑖=𝑘B𝑇 𝛾𝑖.(47) Between these wells, we interpolate the internal energy by means of a fifth-order polynomial 𝐸(𝜉) = 𝑐0+𝑐1𝜉+𝑐2𝜉2+𝑐3𝜉3+𝑐4𝜉4+𝑐5𝜉5,(48) with coefficients obtained imposing the zero, first and second order consistency between (48) and (46a) and (46b), i.e. 𝐸(0) = − 𝜑𝑖, 𝐸(1) = − 𝜑𝑗 𝐸′(0) = 0, 𝐸′(1) = 0 𝐸′′(0) = 𝜅𝑖𝑟2 𝑖𝑗 , 𝐸′′(1) = 𝜅𝑗𝑟2 𝑖𝑗 (49) with the result 𝑐5=1 2𝑟2 𝑖𝑗 (𝜅𝑗−𝜅𝑖) + 6( 𝜑𝑖−𝜑𝑗),(50a) 𝑐4=𝑟2 𝑖𝑗 (3 2𝜅𝑖−𝜅𝑗) − 15( 𝜑𝑖−𝜑𝑗),(50b) 𝑐3=1 2𝑟2 𝑖𝑗 (𝜅𝑗− 3 𝜅𝑖) + 10( 𝜑𝑖−𝜑𝑗),(50c) 𝑐2=1 2𝜅𝑖𝑟2 𝑖𝑗 , 𝑐1= 0, 𝑐0= − 𝜑𝑖(50d) The transition point 𝜉𝑐 then follows as the maximum point of 𝐸(𝜉) in the interval (0,1). The energy barriers follow as 𝐸𝑖→𝑗=𝐸(𝜉𝑐) − 𝐸(0) = 𝐸(𝜉𝑐) + 𝜑𝑖,(51a) 𝐸𝑗→𝑖=𝐸(𝜉𝑐) − 𝐸(1) = 𝐸(𝜉𝑐) + 𝜑𝑗,(51b) and the attempt frequencies as 𝜈𝑖=1 2𝜋√𝜅𝑖 𝑚𝑖 , 𝜈𝑗=1 2𝜋√𝜅𝑗 𝑚𝑗 ,(52) which completes the definition of master Eq. (41). 4. Validation examples We present several examples of validation that test various aspects of the theory, with particular focus on kinetics of hydrogen transport and long-term simulation of microstructure evolution. 4.1. Energy barriers We begin by testing the accuracy of the energy barrier estimates set forth in Section 3.2. To this end, we specifically consider the case of one tetrahedral (T) site located at the center of a (67 × 66 × 62 Å) Mg hcp periodic cell, Fig. 2. Then, we compute the free energy of the system along the transition path between tetrahedral to tetrahedral (T-to-T) and tetrahedral to octahedral (T-to-O) configurations for temperatures in the range 100, ..., 600 K, Figs. 3and 4. We observe that the energy barriers, expresses in terms of the Arrhenius exponential factor 𝛽𝐸𝑏, decrease monotonically with temperature, as expected. The computed energy barriers at 10 K, 0.072 eV (T-to-T) and 0.17 (T-to-O) eV respectively, are in excellent agreement with DFT calculations by Ismer et al. (2009), who report 0.08 eV and 0.19 eV, respectively, and also with MD calculations using NEB, see Vegge (2004). In the case of rutile, we consider a 45 × 32 × 48 Å MgH2 periodic cell with a hydrogen vacancy in the center of the domain Fig. 5. The free-energy variation along transition paths is shown in Fig. 6 in terms of the Arrhenius exponent 𝛽𝐸𝑏. Here again, the computed energy barrier of 0.57 eV at 10 K is in good agreement with the value of 0.65 eV reported by Du et al. (2007) from DFT calculations. We note that all paths exhibit one single energy barrier along the transition. This condition in turn allows the sum Eq. (21) to be restricted to nearest neighbors. The locally quadratic form of the meanfield potential (45) conveniently determines such nearest-neighbors as the sites 𝑗 that share with 𝑖 a face in the Voronoi tessellation of all sites. Fig. 2. Interstitial lattice sites and transition paths for hydrogen in hcp Mg. Magnesium atoms shown in red, tetrahedral and octahedral interstitial hydrogen sites shown in blue and gold, respectively. Black arrows represent the energy diffusion paths considered. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 3. hcp Mg. Temperature dependence of the Arrhenius exponential factor for hydrogen transition between two adjacent tetrahedral sites. Fig. 4. hcp Mg. Temperature dependence of the Arrhenius exponential factor for hydrogen transition between two adjacent octahedral sites. Mechanics of Materials 207 (2025) 105380 6 M. Molinos et al. Fig. 5. Interstitial lattice sites and transition paths for hydrogen in rutile MgH2. Magnesium atoms shown in red, hydrogen sites shown in gold. Black arrows represent the energy diffusion paths considered. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 6. Rutile MgH2. Temperature dependence of the Arrhenius exponential factor for hydrogen transition between two adjacent sites. 4.2. Diffusivity analysis The master Eq. (41) is expected to set forth a diffusive process of mass transport. Estimates of effective diffusivities are therefore often used in order to characterize mass transport properties of materials systems. In molecular dynamics calculations, an analogy to the shortterm asymptotics of the diffusion kernel (mostly variations of the celebrated formula of Varadhan (1967)) is frequently used to define effective diffusivities (see Guinan et al. (1977), also Busch and Paschek (2023) and references therein). In the present setting, local atom-wise diffusivities can be conveniently defined and computed directly from the master Eq. (41) (Venturini et al., 2014). To this end, for a fixed site 𝑖 we make the local ansatz 𝑥𝑗∼𝑥𝑖+ ∇𝑥𝑖⋅𝑟𝑖𝑗 +1 2∇∇𝑥𝑖⋅(𝑟𝑖𝑗 ⊗ 𝑟𝑖𝑗 ) + h. o. t.,(53) where 𝑟𝑖𝑗 =𝑟𝑗−𝑟𝑖 is the relative position vector from sites 𝑖 to 𝑗 and ∇𝑥𝑖 and ∇∇𝑥𝑖 represent macroscopic first and second gradients of the atomic molar fraction, respectively. Formally inserting the ansatz into the master Eq. (41) and collecting terms, we obtain the quasilinear parabolic equation 𝑥𝑖=𝑎𝑖⋅∇∇𝑥𝑖+𝑏𝑖⋅∇𝑥𝑖+𝑐𝑖,(54) Fig. 7. Arrhenius plots for hcp Mg with 1 hydrogen atom. Fig. 8. Arrhenius plots for rutile MgH2 with 1 hydrogen vacancy. with coefficients 𝑎𝑖=1 2∑ 𝑗≠𝑖((1 − 𝑥𝑖)𝑓𝑗→𝑖+𝑥𝑖𝑓𝑖→𝑗)𝑟𝑖𝑗 ⊗ 𝑟𝑖𝑗 ,(55a) 𝑏𝑖=∑ 𝑗≠𝑖((1 − 𝑥𝑖)𝑓𝑗→𝑖+𝑥𝑖𝑓𝑖→𝑗)𝑟𝑖𝑗 ,(55b) 𝑐𝑖=∑ 𝑗≠𝑖 𝑥𝑖(1 − 𝑥𝑖) (𝑓𝑖→𝑗−𝑓𝑗→𝑖),(55c) which depend on the local state and atomic configuration. The terms in (54) account for mass trapping, bias and diffusion. In particular, 𝑐𝑖 is the local rate of trapping, 𝑏𝑖 is a local drift velocity and 𝑎𝑖 is the local diffusivity tensor. We note that 𝑎𝑖 is symmetric and can be anisotropic in general. Reference diffusivity properties can be defined by considering the diffusion of an initially fully occupied and isolated site, corresponding to 𝑥𝑖= 1 and 𝑥𝑗= 0 for 𝑗≠𝑖. In this case, (55) reduces to 𝑎𝑖=1 2∑ 𝑗≠𝑖 𝑓𝑖→𝑗𝑟𝑖𝑗 ⊗ 𝑟𝑖𝑗 ,(56a) 𝑏𝑖=∑ 𝑗≠𝑖 𝑓𝑖→𝑗𝑟𝑖𝑗 , 𝑐𝑖= 0.(56b) Mechanics of Materials 207 (2025) 105380 7 M. Molinos et al. Fig. 9. hcp Mg at 600 K. Time evolution of atomic molar fractions for 1 hydrogen located at the center of the periodic cell. (a) 0.01 ns; (b) 0.03 ns; (c) 0.06 ns. Magnesium sites in gray, atomic molar fractions of hydrogen shown color-coded. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Fig. 10. hcp Mg at 600 K. Time evolution of atomic molar fractions for 5 hydrogens located at random locations in the periodic cell. (a) 0.01 ns; (b) 0.03 ns; (c) 0.06 ns. Magnesium sites in gray, atomic molar fractions of hydrogen shown color-coded. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) If all the sites are embedded in the same local atomic configuration, as is the case, e. g., of a perfect lattice, then a single diffusivity tensor characterizes the system. We also note that diffusion is unbiased at site 𝑖 if ∑ 𝑗≠𝑖 𝑓𝑖→𝑗𝑟𝑖𝑗 = 0,(57) which is ensured, e. g., by centrosymmetry of the crystal lattice. We observe that for simple lattices these reference properties are independent of 𝑖, as expected from translation invariance, and depend only on the local lattice structure and the transport coefficients (39). For complex lattices, the effective diffusivity 𝑎𝑖 may vary from site to site in accordance with the local atomic configuration of the sites. In Fig. 7, we validate predictions for the case of an isolated hydrogen atom in a Mg hcp cell against experimentally measured averaged bulk diffusivities (Nishimura et al., 1999b), Renner and Grabke (1978) and DFT calculations (Klyukin et al., 2015). We observe from the figure that the diffusion coefficients are in good agreement with the experimental data and correctly exhibit the expected hexagonal symmetry of hcp Mg. Fig. 8 shows corresponding diffusion constants for an isolated vacancy in rutile MgH2. Here again, the anisotropy of the diffusion tensor is evident in the figure. In particular, both the hcp MgH and tetragonal rutile phase of MgH2, are predicted to exhibit varying diffusivities parallel and normal to the basal plane. Remarkably, the predicted hydrogen diffusivities in magnesium hydride are much lower than in magnesium, in agreement with experimental data and the MD calculations of Spataru et al. (2020). 4.3. Hydrogenation and dehydrogenation kinetics In order to showcase and assess the ability and efficiency of the theory to predict microstructural evolution, especially for long times, we consider processes involving hydrogenation and dehydrogenation kinetics at various stages of evolution for both magnesium and magnesium hydride as a function of stoichiometry and temperature. For purposes of comparison, we specifically choose configurations identical to those analyzed by Spataru et al. (2020) using MD simulations. The chosen configurations are concerned with dilute concentrations of hydrogen in hcp Mg, and with dilute concentrations of hydrogen vacancies in magnesium hydride. Both configurations are encountered in practice during the operation cycle of hydrogen storage materials. The work of Spataru et al. (2020) can be consulted for an in-depth discussion of the experimentally observed behavior and its physical origins. In the present work, we focus specifically on assessing the efficiency of DMD relative to MD, especially as regards long-term behavior. For each target temperature and pressure, we initially equilibrate the system under NPT conditions and constant initial hydrogen occupancies, and subsequently switch to 𝜇VT conditions. We then update the system in time by means of a staggered approach in which: the meanfield thermomechanical problem (27) is first solved at constant hydrogen occupancy using the critical point line search (Brune et al., 2015) implemented in the PETSc/TAO library (Balay et al., 2023); and the mass transport Eq. (41) is then integrated at constant atomic positions using a backward-Euler implicit integration scheme. The time step is selected so as to resolve the most restrictive relaxation time deduced from the fluxes 𝑓𝑖→𝑗 in (41). 4.3.1. Short-term analysis We begin by considering the short-term hydrogenation of hcp Mg. The computational periodic cell contains 42 (11 20) planes in the 𝑥 direction, 24 (1 100) planes in the 𝑦 direction, which coincides with the 𝑐 axis, and 12 (0001) planes in the 𝑧 direction. The dimension of the resulting computational cell is approximately 67 × 66 × 62 Å, containing 12096 Mg atoms. We consider three dilute initial distributions of hydrogen, consisting of 1, 5 and 121 occupied interstitials in the Mechanics of Materials 207 (2025) 105380 8 M. Molinos et al. Fig. 11. hcp Mg at 600 K. Time evolution of atomic molar fractions for 121 hydrogens located at random locations in the periodic cell. (a) 0.01 ns; (b) 0.03 ns; (c) 0.06 ns. Magnesium sites in gray, atomic molar fractions of hydrogen shown color-coded. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) periodic computational cell and corresponding roughly to stoichiometries X𝐻= 8.27 × 10−5, 4.13 × 10−4 and 0.01, respectively. In order to ascertain the temperature dependence, we repeat the calculations for seven temperatures in the range of 400 to 700 K; and the rutile calculations for nine temperatures in the range of 600 to 800 K. Sequences of snapshots of the early stages of the evolution of the Mg-H system at 600K are shown in Fig. 9, 10 and 11 at times 0.01, 0.03 and 0.06 ns. The time step used in the mass-transport equation calculations is 0.01 ns, which is 2 × 104 larger than the time step of 0.5 fs used in MD calculations by Spataru et al. (2020), or a four orderof-magnitude gain. Thus, the DMD calculations operate on the diffusive time scale, in sharp contrast to the time steps necessitated by molecular dynamics which are controlled by the period of thermal vibration of the lattice. It bears emphasis that this staggering gain is achieved at no loss of fidelity. Indeed, the DMD calculations exhibit the expected phenomenology of mass transport on a lattice, which takes the form of anisotropic discrete diffusion. Thus, in the early stages of diffusion, the initial hydrogen interstitial occupancies defocus into defocus atomic molar fraction distributions of steadily increasing but finite extent, in contrast with classical diffusion in which a disturbance or change in concentration spreads instantly throughout the entire domain, albeit with Gaussian decay. For sufficiently short times, the atomic molar fraction distributions induced by each of the initial hydrogen interstitials do not overlap, but nevertheless interact weakly through their elastic fields, Figs. 9and 10. The predicted anisotropy of hydrogen diffusivity is also evident from the elongated shape of the evolving atomic molar fraction distributions, which reflects [0 0 0 1] as the preferential direction of diffusion, see Fig. 7. For a sufficiently dense initial distribution of hydrogen interstitials, or for sufficiently long times, the atomic molar fraction distributions induced by each of the initial hydrogen interstitials overlap and interact strongly, both elastically and entropically. Fig. 11. The resulting evolution can be characterized by means of a cluster analysis of the atomic molar fraction distributions. Fig. 12. As expected from the discreteness of the process, the evolution of the size of the atomic molar fraction cluster size deviates from classical growth 𝑟∼√𝐷𝑡, and instead initially exhibits growth 𝑟∼𝑡3∕4 followed by ballistic-diffusion growth 𝑟∼𝑡 at intermediate times. 4.3.2. Long-term dehydrogenation Finally, we assess the ability of DMD to characterize long-term kinetics. To this end, we consider the process of dehydrogenation of rutile MgH𝑋. Assuming a tetragonal lattice alignment, the resulting computational cell for the rutile structure of MgH2 contains 16 (010) planes in the 𝑥 direction, 10 (011) planes in the 𝑦 direction, and 14 (0 11) planes in the 𝑧 direction. The dimension of the resulting periodic cell is 50 × 30 × 45 Å, containing 2240 Mg atoms and 4480 H atoms. Fig. 12. Time dependence of the standard-deviation radius of hydrogen density around initially occupied sites and inferred time exponents showing non-classical diffusion. We consider two hydrogen stoichiometries: X𝐻= 1.95 and 1.7. These stoichiometries are initially realized by randomly setting hydrogen occupancies at interstitial sites to 0 (unoccupied) or 1 (occupied). Specifically, we fill 112 and 672 interstitial sites, respectively. As already noted in Section 4.2, vacancy diffusivities in rutile MgH𝑋 are much smaller than hydrogen diffusivities in hcp Mg, resulting in a much slower evolution of the vacancy distribution and, correspondingly, in the need to extend the analysis to much longer times. It is, therefore, crucial that the time step required in calculations to integrate the mass-transport Eq. (41) is 1 ns. Again, it bears emphasis that this time step is chosen to resolve the diffusive time scale and is a staggering six orders of magnitude larger than the time steps required by molecular dynamics (cf., e. g., Spataru et al. (2020)), which are in the fs scale. Fig. 13, 14 and 15 show the evolution of the vacancy distribution up to 10 ns at 𝑇= 600 K for the systems of 1, 112 and 672 hydrogen vacancies, respectively. As in the hcp Mg hydrogenation case, the evolution of the 1-vacancy system reflects clearly the anisotropy of the diffusivity tensor, see Fig. 8. The defocusing of the distributions of the individual vacancies is localized at early times, as expected from discrete diffusion, and evolve independently. At longer times, the distributions of the individual vacancies become delocalized and interact strongly, with a general trend towards a uniform distribution, as expected. The structural stability of the systems is showcased in Figs. 16 and 17, which show [ 110]-views of configurations at times 1 ns and 10 ns Mechanics of Materials 207 (2025) 105380 9