Con en s lis s a ailable a ScienceDi ec
Mechanics o Ma e ials
jou nal homepage: www.else ie .com/loca e/mecma
Resea ch pape
On- he- ly mean ield ansi ion-s a e heo y o di usi e molecula dynamics
M. Molinosa, M. O izb,c, M.P. A izaa,∗
aEscuela Técnica Supe io de Ingenie ía, Uni e sidad de Se illa, Camino de los descub imien os, s.n., 41092, Se illa, Spain
bDi ision o Enginee ing and Applied Science, Cali o nia Ins i u e o Technology, 1200 E. Cali o nia Bl d., Pasadena, CA 91125, USA
cCen e In e nacional de Mé odes Nume ics en Enginye ia (CIMNE), Uni e si a Poli ècnica de Ca alunya, Jo di Gi ona 1, Ba celona, 08034, Spain
A R T I C L E I N F O
Keywo ds:
Magnesium
Magnesium hyd ides
Angula -dependen in e a omic po en ials
Mean ield app oxima ion
T ansi ion s a e heo y
Mass anspo
A B S T R A C T
We apply ansi ion s a e heo y o de i e a omic-le el mas e equa ions o mass anspo om empi ical
in e a omic po en ials wi hin he Di usi e Molecula Dynamics (DMD) amewo k. We show ha mean ield
app oxima ion p o ides an exceedingly e icien and accu a e means o compu ing ee-ene gy ba ie s in
a bi a y local a omic con igu a ions, hus enabling long- e m DMD ‘on- he- ly’ and on he sole basis o an
unde lying in e a omic po en ial, wi hou addi ional modeling assump ions. We apply and alida e he esul ing
mean ield DMD pa adigm in simula ions o p ocesses o hyd ogena ion and dehyd ogena ion o Mg using
Angula -Dependen in e a omic Po en ials (ADP). We show ha mean ield DMD co ec ly p edic s hyd ogen
di usi i ies in hcp Mg and acancy di usi i ies in u ile MgH2. We demons a e he abili y o mean ield
DMD o p edic e olu ion h ough calcula ions conce ned wi h dilu e concen a ions o hyd ogen in hcp Mg,
and wi h dilu e concen a ions o hyd ogen acancies in u ile MgH2, including o -s oichiome y hyd ogen
concen a ions and empe a u e e ec s. Rema kably, he ime s eps equi ed by DMD a e up o six o de s o
magni ude la ge han hose equi ed by Molecula Dynamics (MD), which demons a es he o e whelming
supe io i y o he DMD pa adigm in simula ions o phenomena occu ing on he di usi e ime scale.
1. In oduc ion
Me al hyd ides a e a ac i e o ehicula hyd ogen (Pa ne ship,
2017), hyd ogen s o age (Moh adi and O imo, 2016; Yang e al., 2021),
and many o he eme ging applica ions, and emain he subjec o
ex ensi e ongoing esea ch (Yang e al., 2021; Tan and Ramak ishna,
2021). Hyd ogena ion and dehyd ogena ion p ocesses in me als a e
a e limi ed by hyd ogen up ake and mass anspo wi hin he hos
me al and, he e o e, ope a e on a di usi e ime scale (Luo e al., 2019;
Sh iniwasan and Ta ipa i, 2019; Pund , 2004; Shen and Aguey-Zinsou,
2016). Hyd ogen di usi i y in me als is a complex phenomenon due
o di e ences in c ys al s uc u e be ween he hos me al and i s
hyd ide phases, which ine i ably induces mic os uc u al e olu ion
go e ned by ans o ma ion kine ics. Magnesium hyd ide, he main
ocus o he p esen s udy, is a p ime example o a ans o ma ional
hyd ogen s o age sys em: hexagonal closely packed (hcp) s uc u e
o Mg and e agonal u ile s uc u e o 𝛼-MgH2. Owing o hese
complexi ies, hyd ogen di usi i y canno be easily in e ed om expe -
imen al da a (Bu ge e al., 1961; A ons e al., 1970, 1974; Co nell and
Seymou , 1975; Mazzolai and Zuchne , 1981; Nishimu a e al., 1999a;
Fe nandez and Sanchez, 2002; Ce mak and K al, 2008; Li e al., 2018).
A p edic i e unde s anding o he mechanisms unde lying hyd o-
gena ion and dehyd ogena ion p ocesses in me als he e o e equi es
∗Co esponding au ho .
E-mail add esses: [email p o ec ed] (M. Molinos), [email p o ec ed] (M. O iz), [email p o ec ed] (M.P. A iza).
me hods o analysis ha deli e a omis ic ideli y and, simul aneously,
he abili y o each ac oss o he di usi e ime scale. A he mos
basic le el, he a omis ic models mus accu a ely p edic equilib ium
p ope ies and ansi ion ene gies o hyd ogen anspo wi hin he
hos me al la ice. Densi y Func ional Theo y (DFT) (Tao e al., 2009;
Klyukin e al., 2015), ab ini io me hods (Klyukin e al., 2013; Junkaew
e al., 2014) and empi ical po en ials (Mishin and Lozo oi, 2006;
Smi no a e al., 2018; Zhou e al., 2019a; Molinos e al., 2024) ha e
been widely used o ha end.
Molecula dynamics (MD) has also been used o p edic hyd ogen
di usi i ies in me als and hei hyd ides (Zhou e al., 2016, 2017,
2018b,?; Spa a u e al., 2020). Howe e , MD is se e ely limi ed by he
need o esol e he he mal ib a ions o he a oms, wi h he esul ha
he calcula ions encompass imes ha a e exceedingly sho ela i e
o he di usi e ime scale o in e es . This limi a ion o MD has o en
been sides epped by ecou se o kine ic Mon e Ca lo (KMC) me hods
o pu poses o simula ing anspo phenomena media ed by hopping
ansi ions and o he a e e en s (Bo z e al., 1975; Vo e , 2005;
Ba aile, 2008; Ma inez e al., 2011; Reina e al., 2011). Howe e , KMC
me hods equi e he a p io i enume a ion o all possible ansi ion pa hs
o he sys em and he elucida ion o he co esponding ansi ion a es,
h ps://doi.o g/10.1016/j.mechma .2025.105380
Recei ed 31 Janua y 2025; Recei ed in e ised o m 29 Ap il 2025; Accep ed 29 Ap il 2025
Mechanics o Ma e ials 207 (2025) 105380
A ailable online 19 May 2025
0167-6636/© 2025 The Au ho s. Published by Else ie L d. This is an open access a icle unde he CC BY-NC license ( h p://c ea i ecommons.o g/licenses/by-
nc/4.0/ ).
M. Molinos e al.
which ende s he app oach in ac able when he ansi ion pa hs a e
nume ous and complex, as expec ed in ans o ma ional sys ems. The
same cu se o complexi y bese s e o s o lea n ansi ion pa hs and
ene gy ba ie s using machine lea ning me hods (Tang e al., 2024).
The e o e, he e emains a need o modeling pa adigms capable
o deli e ing ull a omis ic ideli y while simul aneously b idging he
ib a ional and di usi e ime scales. One such pa adigm, p oposed
by Kulka ni (2006), Kulka ni e al. (2008), Ven u ini (2011), Ven u ini
e al. (2014) and e med Di usi e Molecula Dynamics (DMD) by Li
e al. (2011), is p edica ed on a ep esen a ion o he s a e o he
sys em unde going mass anspo as a collec ion o si es ha can be
occupied o emp y. The ene ge ics o he sys em is desc ibed by means
o con en ional empi ical po en ials and a non-equilib ium s a is ical–
mechanical ea men o he ensemble allows he si es o be pa ially
occupied and accoun s o he mal e ec s, possibly including hea
anspo . The e olu ion o he s a e is go e ned by kine ic equa ions
in he spi i o Onsage kine ics, which can be iewed as a omic-le el
Fick and Fou ie laws. The me hodology has been success ully applied
o p oblems o hea anspo (Kulka ni, 2006; Kulka ni e al., 2008;
A iza e al., 2011; Ponga e al., 2015, 2016; Gup a e al., 2021) and
mass anspo (Ven u ini, 2011; Li e al., 2011; Sa ka e al., 2012;
Ven u ini e al., 2014; Ma in e al., 2015; Wang e al., 2015; Simpson
e al., 2016; Sun e al., 2017; Mendez e al., 2018; Mendez and Ponga,
2021).
The main compu a ional challenges inhe en o DMD a e: he com-
pu a ion o phase-space in eg als equi ed o e alua e local equilib ium
ela ions; and he o mula ion o e ec i e kine ic laws om empi ical
in e a omic po en ials. In p e ious wo k, including applica ions o
Mg-H sys ems (Molinos e al., 2024), we ha e shown ha he i s
challenge, he e alua ion o phase-space-in eg als, can be add essed
accu a ely and e icien ly by ecou se o mean ield app oxima ion and
nume ical quad a u e. In combina ion wi h Angula -Dependen in e -
a omic Po en ials (ADP) (Smi no a e al., 2018; Mishin and Lozo oi,
2006) he esul ing he malized and mixed ADP po en ials accu a ely
p edic equilib ium p ope ies o Mg and i s hyd ides, including ee
en opy, hea capaci y, he mal expansion, mola olumes, equa ion o
s a e and elas ic cons an s (Molinos e al., 2024).
In he p esen wo k, we appeal o ansi ion s a e heo y (Weine ,
2012) o de i e a omic-le el mas e equa ions o mass anspo om
empi ical in e a omic po en ials, allowing o a bi a y econ igu a-
ions o he hos me al la ice, wi hin he con ex o he DMD amewo k
o Ven u ini (2011). The undamen al ques ion o be elucida ed con-
ce ns he e icien cha ac e iza ion o a emp equencies and ene gy
ba ie s a endan o hyd ogen ansi ions occu ing in a bi a y –
possibly complex – local a omic con igu a ions, such as ee su aces,
phase bounda ies, g ain bounda ies, amo phized egions and o he s.
The nudged elas ic band me hod (Henkelman and Jonsson, 2000;
Henkelman e al., 2000; Nakano, 2008) is widely used o compu e
ansi ion ene gy ba ie s a 0K. Howe e , he calcula ions a e chal-
lenging and cos ly o complex high-dimensional ene gy landscapes.
In addi ion, o sys ems unde going displaci e phase ansi ions, o
gene ally in he icini y o ex ended la ice de ec s, he numbe and
complexi y o possible ansi ion pa hs is exceedingly la ge and di icul
– i no impossible – o pa ame e ize a p io i.
By con as , we show ha mean ield app oxima ion, combined wi h
he occupancy- a iable ep esen a ion, as ly simpli ies he implemen-
a ion o ansi ion-s a e heo y and supplies an exceedingly e icien
and accu a e means o compu ing ee-ene gy ba ie s in complex
en i onmen s, hus enabling DMD ‘on- he- ly’ and on he sole basis
o an unde lying in e a omic po en ial, wi hou addi ional modeling
assump ions.
By way o demons a ion, we apply and alida e he esul ing
mean ield DMD me hodology o he simula ion o p ocesses o hyd o-
gena ion and dehyd ogena ion o Mg using ADP po en ials (Smi no a
e al., 2018; Molinos e al., 2024). We alida e p edic ions o hyd ogen
di usi i y in hcp Mg agains expe imen ally measu ed a e aged bulk
di usi i ies (Nishimu a e al., 1999b; Renne and G abke, 1978) and
DFT calcula ions (Klyukin e al., 2015). The p edic ed di usi i ies a e
in good ag eemen wi h he expe imen al da a and co ec ly cap u e he
expe imen ally obse ed s ong aniso opy and hexagonal symme y.
We addi ionally alida e p edic ions o acancy di usi i y in u ile
MgH2. He e again, he expe imen ally obse ed s ong aniso opy o
he di usi i ies is cap u ed by he heo y. The p edic ed hyd ogen
di usi i ies in magnesium hyd ide a e much lowe han in magnesium,
also in ag eemen wi h expe imen al da a and he MD calcula ions
o Spa a u e al. (2020).
Finally, in o de o showcase and assess he abili y and e iciency
o mean ield DMD o p edic e olu ion, we p esen calcula ions con-
ce ned wi h dilu e concen a ions o hyd ogen in hcp Mg, and wi h
dilu e concen a ions o hyd ogen acancies in magnesium hyd ide,
including o -s oichiome y hyd ogen concen a ions and empe a u e
e ec s. Bo h p ocesses a e encoun e ed in p ac ice du ing he ope a ion
cycle o hyd ogen s o age ma e ials. In all cases, he compu a ional se -
up is eplica ed om Spa a u e al. (2020) o pu poses o compa ison
wi h MD. We ind ha mean ield DMD accu a ely cap u es bo h he
sho - e m and long- e m esponse o he sys ems, including a ansi ion
om classical o ballis ic di usion, la ice dis o ions esul ing om he
hyd ogen anspo , and la ice s abili y as a unc ion o empe a u e.
Rema kably, in he case o acancy di usion in u ile MgH2 he
ime s eps equi ed by DMD a e six o de s o magni ude la ge ha
hose equi ed by MD, which ep esen s a s agge ing accele a ion and
demons a es he o e whelming supe io i y o he DMD pa adigm o
simula ing phenomena occu ing on he di usi e ime scale.
2. Local equilib ium ela ions and mean ield app oxima ion
By way o backg ound and o se he amewo k, we summa-
ize salien aspec s o he non-equilib ium s a is ical mechanics heo y
o Ven u ini e al. (2014) and i s applica ion o Mg-H sys ems (Molinos
e al., 2024).
2.1. Non-equilib ium s a is ical mechanics
We conside a closed sys em consis ing o 𝑁 si es, e. g., a omic
posi ions o molecules, each o which can be o one o 𝑀 species,
including acancies. Fo each si e 𝑖= 1,…, 𝑁, and each species 𝑘=
1,…, 𝑀, we in oduce he occupancy unc ion
𝑛𝑖𝑘 ={1,i si e 𝑖 is occupied by species 𝑘,
0,o he wise,(1)
in o de o desc ibe he occupancy o each si e. We no e ha , om
de ini ion (1), we mus ha e
𝑀
∑
𝑘=1
𝑛𝑖𝑘 = 1 (2)
a e e y si e 𝑖. We addi ionally deno e by 𝑛𝑖= (𝑛𝑖𝑘)𝑀
𝑘=1 he local
occupancy a ay o si e 𝑖. The mic oscopic s a es o he sys em a e
de ined by he ins an aneous posi ion {𝑞} = (𝑞𝑖)𝑁
𝑖=1, momen a {𝑝} =
(𝑝𝑖)𝑁
𝑖=1, and occupancy a ays {𝑛} = (𝑛𝑖)𝑁
𝑖=1 o all 𝑁 si es in he sys em.
The occupancy unc ions 𝑛𝑖 ake alues in a se 𝑀 consis ing o he
elemen s o {0,1}𝑀 ha sa is y he cons ain (2). In addi ion, he
occupancy a ays {𝑛} ake alues in he se
𝑁𝑀 ={{𝑛} ∈ {0,1}𝑁𝑀 ∶
𝑛𝑖∈𝑀 o 𝑖= 1,…, 𝑁 }.
(3)
The expec ed o mac oscopic alue o a unc ion 𝐴({𝑞},
{𝑝},{𝑛}) is gi en by he phase a e age
⟨𝐴⟩=∑
{𝑛}∈𝑁𝑀
1
ℎ3𝑁×
∫𝛤
𝐴({𝑞},{𝑝},{𝑛}) 𝜌({𝑞},{𝑝},{𝑛}) 𝑑𝑞 𝑑𝑝,
(4)
Mechanics o Ma e ials 207 (2025) 105380
2
M. Molinos e al.
whe e 𝛤= (R3×R3)𝑁, ℎ is Planck’s cons an and ℎ−3𝑁 supplies he
na u al uni o phase olume o sys ems o dis inguishable pa icles (Hill,
1987; Gi i alco, 2000).
We assume ha he s a is ics o he sys em obeys Jaynes’ p inciple
o maximum en opy (Jaynes, 1957a,b; Zuba e , 1974; Callen, 1985),
which posi s ha he p obabili y densi y unc ion 𝜌({𝑞},{𝑝},{𝑛}), cha -
ac e izing he p obabili y o inding he sys em in a s a e ({𝑞},{𝑝},{𝑛}),
maximizes he in o ma ion- heo e ical en opy
[𝜌] = −𝑘𝐵⟨log 𝜌⟩,(5)
among all p obabili y measu es consis en wi h he cons ain s on he
sys em. In (5) and subsequen ly, 𝑘𝐵 deno es Bol zmann’s cons an .
We conside sys ems consis ing o dis inguishable pa icles whose
Hamil onians ha e he addi i e s uc u e
𝐻({𝑞},{𝑝},{𝑛}) =
𝑁
∑
𝑖=1
ℎ𝑖({𝑞},{𝑝},{𝑛}),(6)
whe e ℎ𝑖 is he local Hamil onian o pa icle 𝑖. Following Feynman
(1998), suppose ha he expec ed pa icle ene gies and he expec ed
pa icle a omic ac ions
⟨ℎ𝑖⟩=𝑒𝑖,⟨𝑛𝑖𝑘⟩=𝜒𝑖𝑘,(7)
a e known, espec i ely. We no e ha he local a omic ac ions sa is y
he iden i ies
𝑀
∑
𝑘=1
𝜒𝑖𝑘 = 1,(8)
which ollow om (2) and he second o (7). En o cing he cons ain s
(7) by means o Lag ange mul iplie s 𝑘𝐵{𝛽}≡(𝑘𝐵𝛽𝑖)𝑁
𝑖=1 and {𝛾}≡
((𝛾𝑖𝑘)𝑀
𝑘=1)𝑁
𝑖=1, leads o he Lag angian
[𝜌, {𝛽},{𝛾}] =
[𝜌] − 𝑘𝐵{𝛽}𝑇{⟨ℎ⟩} + 𝑘𝐵{𝛾}𝑇{⟨𝑛⟩},(9)
whe e we in e p e
𝑇𝑖=1
𝑘𝐵𝛽𝑖
,(10)
as he local empe a u e o pa icle 𝑖. In iew o iden i ies (8), we
addi ionally append he cons ain s
𝑀
∑
𝑘=1
𝛾𝑖𝑘 = 0,(11)
in o de o ende {𝛾} de e mina e. Maximizing [⋅,{𝛽},{𝛾}] among
p obabili y measu es gi es
𝜌({𝑞},{𝑝},{𝑛}) = 1
𝛯e−{𝛽}𝑇{ℎ}+{𝛾}𝑇{𝑛},(12)
wi h
𝛯=∑
{𝑛}∈𝑁𝑀
1
ℎ3𝑁∫𝛤
e−{𝛽}𝑇{ℎ}+{𝛾}𝑇{𝑛}𝑑𝑞 𝑑𝑝. (13)
The co esponding equilib ium alues o {𝛽} and {𝛾} ollow om
(7) as a unc ion o {𝑒} and {𝜒}. We in e p e (12) and (13) as
non-equilib ium gene aliza ions o he Gibbs g and-canonical p obabili y
densi y unc ion and he g and-canonical pa i ion unc ion, espec i ely.
2.2. Mean ield app oxima ion
The calcula ion o he he modynamic po en ials in closed o m
is gene ally in ac able and app oxima ion is he e o e equi ed. A
a ia ional mean ield heo y (S anley, 1971; Yeomans, 1992a) may
be o mula ed by es ic ing (9) o some class o p obabili y densi y
unc ions o he o m
𝜌0({𝑞},{𝑝},{𝑛}) = 1
𝛯0
𝑒−{𝛽}𝑇{
ℎ}+{𝛾}𝑇{𝑛},(14)
wi h
𝛯0=∑
{𝑛}∈𝑁𝑀
1
ℎ3𝑁∫𝛤
𝑒−{𝛽}𝑇{
ℎ}+ {𝛾}𝑇{𝑛}𝑑𝑞 𝑑𝑝. (15)
and {
ℎ} in some class 0 o local ial Hamil onians, possibly de ined
pa ame ically. The es ic ed Lag angian is
[𝜌0,{𝛽},{𝛾}] = [𝜌0]−
𝑘B{𝛽}𝑇{⟨ℎ⟩0} + 𝑘B{𝛾}𝑇{⟨𝑛⟩0},(16)
whe e
[𝜌0] = −𝑘𝐵⟨log 𝜌0⟩0(17)
and ⟨⋅⟩0 deno es a e aging wi h espec o 𝜌0. Inse ing (14) and (15)
in o (16) using (17) gi es
[𝜌0,{𝛽},{𝛾}] = 𝑘Blog 𝛯0−𝑘B{𝛽}𝑇{⟨ℎ−
ℎ⟩0}.(18)
The op imal ial Hamil onians a e de e mined by maximizing
[𝜌0,{𝛽},{𝛾}] wi h espec o 𝜌0, o some sui able pa ame iza ion
o {
ℎ} he eo . In addi ion, he co esponding mean ield equilib ium
alues o {𝛽} and {𝛾} ollow as a unc ion o {𝑒} and {𝜒} om he
Eule –Lag ange equa ions o [𝜌0,{𝛽},{𝛾}],
⟨
ℎ𝑖⟩0=𝑒𝑖,⟨𝑛𝑖𝑘⟩0=𝜒𝑖𝑘.(19)
2.3. In e s i ial hyd ogen in me als
As a special case o he gene al heo y jus ou lined, conside a
c ys al la ice whe e he base la ice si es a e always occupied by me al
a oms, while he in e s i ial si es a e ei he occupied by an H a om, o
unoccupied (Sun e al., 2017). We index he base la ice si es by an
index se 𝐼M, and he in e s i ial si es by 𝐼H. Unde hese assump ions,
we can cha ac e ize occupancy by a single occupancy numbe 𝑛𝑖 on he
in e s i ial si es, 𝑖∈𝐼H, aking he alue o 0 i he si e unoccupied and
1 i he si e occupied. We assume a Hamil onian o he o m
𝐻({𝑞},{𝑝},{𝑛}) =
∑
𝑖∈𝐼M
1
2𝑚M|𝑝𝑖|2+∑
𝑖∈𝐼H
1
2𝑚H|𝑝𝑖|2+𝑉({𝑞},{𝑛}),(20)
whe e 𝑚M and 𝑚H a e he a omic masses o he me al hos and hy-
d ogen, espec i ely, and 𝑉({𝑞},{𝑛}) is a mixed angula -dependen
in e a omic po en ial (ADP) o he o m
𝑉({𝑞},{𝑛}) = ∑
𝑖
𝑛𝑖(𝜌𝑖) + 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≠𝑖(𝑛𝑖𝑛𝑗1 𝜆𝑖𝑗1)(𝑛𝑖𝑛𝑗2 𝜆𝑖𝑗2),
(21)
which belongs o a class o many-body po en ials p oposed by Mishin
and Lozo oi (2006). In (21), he sub-index 𝑖 deno es each a omic si e
o he domain while sub-indexes 𝑗, 𝑗1 o 𝑗2 a e used o enume a e he
neighbo s o each si e 𝑖 and 𝑖(𝑟c) deno es he local neighbo hood o 𝑖
wi h cu o adius 𝑟c. The i s e m is he embedding ene gy, he second
e m is he pai wise-in e ac ion ene gy (𝑖−𝑗), he hi d o i h e ms a e
he con ibu ions o he ene gy due o angula in e ac ions (𝑖−𝑗1−𝑗2)
whe e
𝜇𝑖𝑗 =𝑢𝑖𝑗 𝐫𝑖𝑗 and 𝜆𝑖𝑗 =𝑤𝑖𝑗 𝐫𝑖𝑗 ⊗𝐫𝑖𝑗 (22)
a e espec i ely he dipole and quad upole. Whe e 𝐫𝑖𝑗 is he dis ance
ec o be ween a pa icle 𝑖 and i s neighbo 𝑗. The scala unc ions
Mechanics o Ma e ials 207 (2025) 105380
3
M. Molinos e al.
(𝜌𝑖), 𝜙(𝑟𝑖𝑗 ), 𝑢(𝑟𝑖𝑗 ) and 𝑤(𝑟𝑖𝑗 ) a e spline unc ions commonly ep e-
sen ed in abula o m e. g. Smi no a e al. (2018). These unc ions a e
exp essed in e ms o he dis ance no m be ween wo si es, |𝐫𝑖𝑗 |=𝑟𝑖𝑗 ,
and he elec on densi y 𝜌𝑖. See also Molinos e al. (2024) o u he
de ails o he nume ical ea men and o a alida ion assessmen o
he po en ial.
The po en ial ene gy (21) is a unc ion o he hyd ogen occupancies
and, in ha sense, may be ega ded as a h ee-dimensional Ising model.
Such models ha e been ex ensi ely s udied by a a ie y o means,
including mean ield heo y (c ., e. g., Yeomans (1992b)).
Building on ha backg ound, we assume ha he sys em is closed
and in he mal equilib ium, whence 𝛽𝑖=𝛽= 1∕𝑘B𝑇 o all pa icles,
and choose ial Hamil onians o he o m

ℎ𝑖(𝑞𝑖, 𝑝𝑖) = 𝑘B𝑇
2𝜎2
𝑖|𝑞𝑖−𝑞𝑖|2+1
2𝑚𝑖|𝑝𝑖−𝑝𝑖|2,(23)
whe e 𝑞𝑖, 𝑝𝑖 and 𝜎𝑖 a e pa ame e s ha cha ac e ize he ial space. The
co esponding mean ield p obabili y densi y unc ion (14) e alua es o
𝜌0({𝑞},{𝑝},{𝑛}) = 1
𝛯0
exp {
−∑
𝑖∈𝐼M∪𝐼H
1
2𝜎2
𝑖|𝑞𝑖−𝑞𝑖|2−∑
𝑖∈𝐼M
𝛽
2𝑚M|𝑝𝑖−𝑝𝑖|2
−∑
𝑖∈𝐼H(𝛽
2𝑚H|𝑝𝑖−𝑝𝑖|2−𝛾𝑖𝑛𝑖)},
(24)
and he mean ield g and-canonical pa i ion unc ion (15) o
𝛯0={∏
𝑖∈𝐼M(𝜎𝑖√𝑚M∕𝛽
ℏ)3}×
{∏
𝑖∈𝐼H(𝜎𝑖√𝑚H∕𝛽
ℏ)3(1+e𝛾𝑖)},
(25)
wi h ℏ he educed Planck cons an . Unde he assumed iso he mal
condi ions, he mean ield Lag angian (16) educes o
[𝜌0,{𝛽},{𝛾}] = 𝑘Blog 𝛯0−
𝑘𝐵𝛽(⟨𝑉⟩0+∑
𝑖∈𝐼M∪𝐼H
1
2𝑚𝑖|𝑝𝑖|2)+ 3 𝑁 𝑘B,(26)
whe e 𝑁= #𝐼M+ #𝐼H is he o al numbe o si es in he sys em. The
co esponding mean ield op imali y condi ions a e
−𝜕
𝜕 𝑞𝑖
=𝜕
𝜕 𝑞𝑖⟨𝑉⟩0=⟨𝜕𝑉
𝜕𝑞𝑖⟩0= 0,(27a)
−𝜕
𝜕 𝑝𝑖
=1
𝑚𝑖
𝑝𝑖= 0,(27b)
−𝜕
𝜕 𝜎𝑖
= − 3𝑘B
𝜎𝑖
+𝑘B𝛽𝜕
𝜕 𝜎𝑖⟨𝑉⟩0= 0,(27c)
and he mean ield Eule –Lag ange Eqs. (19) e alua e o
⟨𝑉𝑖⟩0+|𝑝𝑖|2
2𝑚𝑖
=𝑒𝑖, 𝑖 ∈𝐼M∪𝐼H,(28a)
e𝛾𝑖
1+e𝛾𝑖=𝜒𝑖, 𝑖 ∈𝐼H.(28b)
Al e na i ely, Eq. (28b) can be in e ed o yield
𝛾𝑖= log(𝜒𝑖
1 − 𝜒𝑖), 𝑖 ∈𝐼H.(29)
The equilib ium ela ion (29) is plo ed in Fig. 1 o a uni o m H a omic
ac ion in a pe ec hcp la ice and a ange o empe a u es.
We ecall om (4) ha
⟨𝑉⟩0=∑
{𝑛}∈𝑁𝑀
1
ℎ3𝑁×
∫𝛤
𝑉({𝑞},{𝑛})𝜌0({𝑞},{𝑝},{𝑛}) 𝑑𝑞 𝑑𝑝.
(30)
Fig. 1. Local equilib ium ela ion (29) exp essed in e ms o he local chemical
po en ial 𝜇𝑖∶= 𝛾𝑖∕𝛽𝑖 and he local hyd ogen mola ac ion 𝜒𝑖 o H a di e en
empe a u es.
In iew o (24) and wi h e e ence o Jensen’s inequali y, we app oxi-
ma e (30) as
⟨𝑉⟩0≈∫𝑉({𝑞},{𝜒})
{∏
𝑗∈𝐼M∪𝐼H
1
(√2𝜋 𝜎𝑗)3
exp(−1
2𝜎2
𝑗|𝑞𝑗−𝑞𝑗|2)}𝑑𝑞.
(31)
as p oposed by Ven u ini e al. (2014) o a oid occupancy sums o
combina o ial complexi y. In addi ion, we app oxima e he emaining
in eg al o e con igu a ion space by means o nume ical quad a u e,
see Molinos e al. (2024) o de ails o he nume ical implemen a ion
and e i ica ion he eo .
3. Con igu a ion-dependen hyd ogen anspo kine ics
In o de o o mula e a gene al amewo k o kine ics, including
hea and mass anspo , we begin by examining he balance o ene gy
a he pa icle le el (Ven u ini e al., 2014). The in e nal ene gy 𝑒𝑖 o
pa icle 𝑖 may be iden i ied wi h he expec ed alue o he co espond-
ing pa icle Hamil onian as in he i s o (7). Suppose, in addi ion,
ha he local Hamil onians ℎ𝑖 depend on mac oscopic a iables {} =
(1,…,𝜈), such as olume o de o ma ion (Weine , 2002). Then, he
balance o ene gy a pa icle 𝑖 can be exp essed as
𝑒𝑖=𝜇𝑇
𝑖𝜒𝑖+
𝑤𝑖+𝑟𝑖,
𝑤𝑖=
𝜈
∑
𝛼=1⟨𝜕ℎ𝑖
𝜕𝛼

𝛼⟩(32)
whe e
𝜇𝑖=𝛾𝑖∕𝛽𝑖(33)
is pa icle chemical po en ial, 
𝑤𝑖 is he ex e nal mechanical powe and
𝑟𝑖 is he hea low in o pa icle 𝑖. We assume ha he hea low 𝑟𝑖 in o
pa icle 𝑖 is o he o m
𝑟𝑖=∑
𝑗≠𝑖
𝑅𝑖𝑗 , 𝑅𝑖𝑗 = −𝑅𝑗𝑖,(34)
whe e he sum ex ends o all pa icles di e en om 𝑖, and 𝑅𝑖𝑗 is he
disc e e hea lux om pa icle 𝑗 o 𝑖. Likewise, we shall assume ha he
mass low in o pa icle 𝑖 may be exp essed as
𝑥𝑖=∑
𝑗≠𝑖
𝐽𝑖𝑗 , 𝐽𝑖𝑗 = −𝐽𝑗𝑖,(35)
Mechanics o Ma e ials 207 (2025) 105380
4
M. Molinos e al.
whe e 𝐽𝑖𝑗 is he disc e e mass lux a ay om pa icle 𝑗 o 𝑖. The in e nal
en opy p oduc ion a e o a pa icle pai is gi en by
𝛴𝑖𝑗 =𝐾𝑇
𝑖𝑗 𝐽𝑖𝑗 +𝑃𝑖𝑗 𝑅𝑖𝑗 ≥0,(36)
whe e we w i e
𝐾𝑖𝑗 =𝜇𝑖
𝑇𝑖
−𝜇𝑗
𝑇𝑗
=𝑘𝐵(𝛾𝑖−𝛾𝑗),
𝑃𝑖𝑗 =1
𝑇𝑖
−1
𝑇𝑗
=𝑘𝐵(𝛽𝑖−𝛽𝑗).
(37)
Following Onsage (1931a,b), he local dissipa ion inequali y (36) sug-
ges s kine ic laws o he gene al o m
𝐽𝑖𝑗 = − 𝜕𝛹
𝜕𝐾𝑖𝑗
({𝑃},{𝐾}),
𝑅𝑖𝑗 =𝜕𝛹
𝜕𝑃𝑖𝑗
({𝑃},{𝐾}),
(38)
whe e 𝛹({𝑃},{𝐾}) is a disc e e kine ic po en ial.
The abili y o models based on non-equilib ium s a is ical mechan-
ics o ep oduce he obse ed aniso opy, empe a u e and size de-
pendence o he he mal conduc i i y o Si nanowi es was es ablished
by Ma in e al. (2015) by way o alida ion o he heo y. The models
ha e also demons a ed p edic i e abili y in applica ions including
nano oid g ow h in me als a low and high s ain a es (A iza e al.,
2011; Ponga e al., 2017).
3.1. Mass anspo
Nex , we specialize he gene al heo y o mass anspo . To his
end, we conside a ansi ion om an ini ial equilib ium s a e 𝑠∶=
({𝑞},{𝑛}) o a inal equilib ium s a e 𝑠′∶= ({𝑞′},{𝑛′}). The wo s a es
di e in he occupancy numbe s, some o which may ha e lipped, and
he posi ions o he a oms. The ansi ion 𝑠′→𝑠 is assumed o ollow a
con inuous ansi ion pa h 𝑠(𝜉), pa ame e ized by 𝜉∈ [0,1], wi h 𝑠(0) =
𝑠′ and 𝑠(1) = 𝑠. Fo he ene gy o he sys em o e ol e con inuously
h oughou he ansi ion, pai s o si es 𝑖 and 𝑗, one ini ially occupied,
𝑛′
𝑖= 1, and one ini ially unoccupied, 𝑛′
𝑗= 0, mus exchange mass,
𝑛′
𝑖→𝑛𝑖= 0 and 𝑛′
𝑗→𝑛𝑗= 1, a he poin 𝜉𝑐 o hei ajec o y when
hei si es a e coinciden , 𝑞𝑖(𝜉𝑐) = 𝑞𝑗(𝜉𝑐).
Acco ding o ansi ion s a e heo y (Weine , 2012), he a e age
lipping a e is gi en by he A henius ela ion
𝑓𝑖→𝑗=𝜈𝑖e−𝛽𝐸𝑖→𝑗.(39)
whe e 𝜈𝑖 is he a emp equency o si e 𝑖 and 𝐸𝑖→𝑗 is he ene gy
ba ie o 𝑖→𝑗 hops. Fu he in e p e ing he a omic ac ions 𝑥𝑖
as occupancy p obabili ies, he p obabili y ha si e 𝑖 be occupied is
p ecisely 𝑥𝑖 and he p obabili y ha he si e 𝑗 be unoccupied is (1−𝑥𝑗),
gi ing a ansi ion p obabili y
𝜓𝑖→𝑗=𝑥𝑖(1 − 𝑥𝑗)𝑓𝑖→𝑗.(40)
The a omic ac ion a e is hen gi en by he mas e equa ion (No d-
sieck e al., 1940)
𝑥𝑖=∑
𝑗≠𝑖(𝜓𝑗→𝑖−𝜓𝑖→𝑗).(41)
We e i y om (41) ha
∑
𝑖∈𝐼𝐻
𝑥𝑖= 0,(42)
i. e., he o al mass o he sys em is conse ed in he absence o sou ces
and sinks, as equi ed o a chemically isola ed sys em.
3.2. Calcula ion o ene gy ba ie s
We see om (39) ha , wi hin he p esen amewo k, he o mu-
la ion o kine ic models o mass anspo equi es he speci ica ion o
sui able o ms o 𝜈𝑖 and 𝐸𝑖→𝑗 and, speci ically, o hei dependence on
he a omic con igu a ion, e. g., h ough he local en i onmen .
The calcula ion o a emp equencies and ene gy ba ie s is he
main ocus o ansi ion-s a e heo y (Weine , 2002) and has adi ional
been a main compu a ional bo leneck in p ac ice. Vineya d’s o -
mula (Vineya d, 1957), which is de i ed om a s a is ical–mechanical
analysis o ha monic app oxima ions o he po en ial ene gy abou
ene gy wells and saddle poin s, and nume ous ex ensions and enhance-
men s he eo , a e widely used in p ac ice. We show nex ha he use
o occupancy a iables and mean ield app oxima ion as ly simpli ies
he implemen a ion o ansi ion-s a e heo y and enables he on- he- ly
e alua ion o ene gy ba ie s and a emp equencies, hus ende ing
he app oach p ac ical.
3.2.1. Hessian algo i hm and Vineya d’s o mula
We begin by specializing Vineya d’s o mula (Vineya d, 1957)
o sys ems o pa icles desc ibed by means o occupancy a iables.
We con ine a en ion o ansi ion pa hs such ha he o al ene gy
𝐸({𝑞(𝜉)},{𝑛(𝜉)}), pa ame e ized by an o de pa ame e 𝜉∈ [0,1],
a ains local minima a he ini ial and inal s a es and he minima
a e sepa a ed by one single maximum, o ene gy ba ie . Among such
ansi ion pa hs, we seek o de e mine ha o which he ene gy ba ie
is smalles . By his condi ion, he ene gy maximum necessa ily occu s
a a saddle poin , ({𝑞(𝜉𝑐)},{𝑛(𝜉𝑐)}), o ansi ion s a e.
In o de o ind he posi ion o he saddle poin be ween wo local
minima 𝑖 and 𝑗, we look o he maximum be ween 𝑖 and 𝑗 using he
sea ch di ec ion gi en by he local g adien and Hessian o he o al
ene gy po en ial. Once a local maximum poin is ound, we es whe he
i is indeed a saddle poin o he po en ial.
The sough ene gy ba ie is hen
𝐸𝑏=𝐸({𝑞(𝜉𝑐)},{𝑛(𝜉𝑐)}) − 𝐸({𝑞(0)},{𝑛(0)}),(43)
and he a emp equency o si e 𝑖 is gi en by Vineya d’s o mula (Vine-
ya d, 1957)
𝜈𝑖=1
2𝜋
de [𝐷2𝐸({𝑞(0)},{𝑛(0)})]
de [𝐷2𝐸({𝑞(𝜉𝑐)},{𝑛(𝜉𝑐)})+].(44)
In his exp ession, 𝐷2𝐸({𝑞},{𝑛}) is he ma ix o second de i a i es,
o Hessian, o 𝐸 a ({𝑞},{𝑛}), and he subsc ip ()+ designa es i s posi-
i e componen , i. e., he componen ob ained by conside ing posi i e
eigen alues only.
3.2.2. Mean ield calcula ion o ene gy ba ie s
Despi e he appeal o Vineya d’s o mula, and ex ensions he eo , i s
main d awback is ha he de e mina ion o ansi ion pa hs be ween
local ene gy minima is exceedingly di icul in gene al and compu a-
ionally expensi e o la ge sys ems. We o e come his di icul y by
ins ead es ima ing he ansi ions by ecou se o a ia ional mean ield
heo y.
By consis ency wi h he mean ield p obabili y densi y (24), we
app oxima e he o al in e nal ene gy as
𝐸({𝑞},{𝑛}) = ∑
𝑖∈𝐼
𝑘B𝑇
2𝜎2
𝑖|𝑞𝑖−𝑞𝑖|2−∑
𝑖∈𝐼𝐻
𝑘B𝑇 𝛾𝑖𝑛𝑖,(45)
whe e {𝑞} and {𝜎} a e he mean ield a iables in oduced in Sec-
ion 2.2. We hen conside ansi ions be ween pai s o hyd ogen si es
𝑖, 𝑗∈𝐼𝐻, 𝑖≠𝑗, and assume ha he emaining occupancies 𝜒𝑘 emain
unchanged, and ha he emaining pa icle posi ions 𝑞𝑘 emain close
o hei a e age alues 𝑞𝑘, 𝑘∈𝐼𝐻, 𝑖≠𝑘≠𝑗, h ough he ansi ion.
Unde hese assump ions, om (45) we ha e
𝐸(𝜉)∼ 𝜅𝑖
2|𝑞𝑖(𝜉) − 𝑞𝑖|2−𝜑𝑖, o 𝜉∼ 0,(46a)
Mechanics o Ma e ials 207 (2025) 105380
5

M. Molinos e al.
𝐸(𝜉)∼ 𝜅𝑗
2|𝑞𝑗(𝜉) − 𝑞𝑗|2−𝜑𝑗, o 𝜉∼ 1,(46b)
whe e we w i e
𝜅𝑖=𝑘B𝑇
𝜎2
𝑖
, 𝜑𝑖=𝑘B𝑇 𝛾𝑖.(47)
Be ween hese wells, we in e pola e he in e nal ene gy by means o a
i h-o de polynomial
𝐸(𝜉) = 𝑐0+𝑐1𝜉+𝑐2𝜉2+𝑐3𝜉3+𝑐4𝜉4+𝑐5𝜉5,(48)
wi h coe icien s ob ained imposing he ze o, i s and second o de
consis ency be ween (48) and (46a) and (46b), i.e.
𝐸(0) = − 𝜑𝑖, 𝐸(1) = − 𝜑𝑗
𝐸′(0) = 0, 𝐸′(1) = 0
𝐸′′(0) = 𝜅𝑖𝑟2
𝑖𝑗 , 𝐸′′(1) = 𝜅𝑗𝑟2
𝑖𝑗
(49)
wi h he esul
𝑐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 ansi ion poin 𝜉𝑐 hen ollows as he maximum poin o 𝐸(𝜉) in
he in e al (0,1). The ene gy ba ie s ollow as
𝐸𝑖→𝑗=𝐸(𝜉𝑐) − 𝐸(0) = 𝐸(𝜉𝑐) + 𝜑𝑖,(51a)
𝐸𝑗→𝑖=𝐸(𝜉𝑐) − 𝐸(1) = 𝐸(𝜉𝑐) + 𝜑𝑗,(51b)
and he a emp equencies as
𝜈𝑖=1
2𝜋√𝜅𝑖
𝑚𝑖
, 𝜈𝑗=1
2𝜋√𝜅𝑗
𝑚𝑗
,(52)
which comple es he de ini ion o mas e Eq. (41).
4. Valida ion examples
We p esen se e al examples o alida ion ha es a ious aspec s
o he heo y, wi h pa icula ocus on kine ics o hyd ogen anspo
and long- e m simula ion o mic os uc u e e olu ion.
4.1. Ene gy ba ie s
We begin by es ing he accu acy o he ene gy ba ie es ima es se
o h in Sec ion 3.2. To his end, we speci ically conside he case o one
e ahed al (T) si e loca ed a he cen e o a (67 × 66 × 62 Å) Mg hcp
pe iodic cell, Fig. 2. Then, we compu e he ee ene gy o he sys em
along he ansi ion pa h be ween e ahed al o e ahed al (T- o-T)
and e ahed al o oc ahed al (T- o-O) con igu a ions o empe a u es
in he ange 100, ..., 600 K, Figs. 3and 4. We obse e ha he ene gy
ba ie s, exp esses in e ms o he A henius exponen ial ac o 𝛽𝐸𝑏,
dec ease mono onically wi h empe a u e, as expec ed. The compu ed
ene gy ba ie s a 10 K, 0.072 eV (T- o-T) and 0.17 (T- o-O) eV espec-
i ely, a e in excellen ag eemen wi h DFT calcula ions by Isme e al.
(2009), who epo 0.08 eV and 0.19 eV, espec i ely, and also wi h
MD calcula ions using NEB, see Vegge (2004).
In he case o u ile, we conside a 45 × 32 × 48 Å MgH2 pe iodic
cell wi h a hyd ogen acancy in he cen e o he domain Fig. 5.
The ee-ene gy a ia ion along ansi ion pa hs is shown in Fig. 6 in
e ms o he A henius exponen 𝛽𝐸𝑏. He e again, he compu ed ene gy
ba ie o 0.57 eV a 10 K is in good ag eemen wi h he alue o 0.65 eV
epo ed by Du e al. (2007) om DFT calcula ions.
We no e ha all pa hs exhibi one single ene gy ba ie along he
ansi ion. This condi ion in u n allows he sum Eq. (21) o be e-
s ic ed o nea es neighbo s. The locally quad a ic o m o he mean ield
po en ial (45) con enien ly de e mines such nea es -neighbo s as he
si es 𝑗 ha sha e wi h 𝑖 a ace in he Vo onoi essella ion o all si es.
Fig. 2. In e s i ial la ice si es and ansi ion pa hs o hyd ogen in hcp Mg. Magnesium
a oms shown in ed, e ahed al and oc ahed al in e s i ial hyd ogen si es shown in blue
and gold, espec i ely. Black a ows ep esen he ene gy di usion pa hs conside ed.
(Fo in e p e a ion o he e e ences o colo in his igu e legend, he eade is e e ed
o he web e sion o his a icle.)
Fig. 3. hcp Mg. Tempe a u e dependence o he A henius exponen ial ac o o
hyd ogen ansi ion be ween wo adjacen e ahed al si es.
Fig. 4. hcp Mg. Tempe a u e dependence o he A henius exponen ial ac o o
hyd ogen ansi ion be ween wo adjacen oc ahed al si es.
Mechanics o Ma e ials 207 (2025) 105380
6
M. Molinos e al.
Fig. 5. In e s i ial la ice si es and ansi ion pa hs o hyd ogen in u ile MgH2.
Magnesium a oms shown in ed, hyd ogen si es shown in gold. Black a ows ep esen
he ene gy di usion pa hs conside ed. (Fo in e p e a ion o he e e ences o colo in
his igu e legend, he eade is e e ed o he web e sion o his a icle.)
Fig. 6. Ru ile MgH2. Tempe a u e dependence o he A henius exponen ial ac o o
hyd ogen ansi ion be ween wo adjacen si es.
4.2. Di usi i y analysis
The mas e Eq. (41) is expec ed o se o h a di usi e p ocess o
mass anspo . Es ima es o e ec i e di usi i ies a e he e o e o en
used in o de o cha ac e ize mass anspo p ope ies o ma e ials
sys ems. In molecula dynamics calcula ions, an analogy o he sho -
e m asymp o ics o he di usion ke nel (mos ly a ia ions o he
celeb a ed o mula o Va adhan (1967)) is equen ly used o de ine
e ec i e di usi i ies (see Guinan e al. (1977), also Busch and Paschek
(2023) and e e ences he ein).
In he p esen se ing, local a om-wise di usi i ies can be con e-
nien ly de ined and compu ed di ec ly om he mas e Eq. (41) (Ven-
u ini e al., 2014). To his end, o a ixed si e 𝑖 we make he local
ansa z
𝑥𝑗∼𝑥𝑖+ ∇𝑥𝑖⋅𝑟𝑖𝑗 +1
2∇∇𝑥𝑖⋅(𝑟𝑖𝑗 ⊗ 𝑟𝑖𝑗 ) + h. o. .,(53)
whe e 𝑟𝑖𝑗 =𝑟𝑗−𝑟𝑖 is he ela i e posi ion ec o om si es 𝑖 o 𝑗 and
∇𝑥𝑖 and ∇∇𝑥𝑖 ep esen mac oscopic i s and second g adien s o he
a omic mola ac ion, espec i ely. Fo mally inse ing he ansa z in o
he mas e Eq. (41) and collec ing e ms, we ob ain he quasilinea
pa abolic equa ion
𝑥𝑖=𝑎𝑖⋅∇∇𝑥𝑖+𝑏𝑖⋅∇𝑥𝑖+𝑐𝑖,(54)
Fig. 7. A henius plo s o hcp Mg wi h 1 hyd ogen a om.
Fig. 8. A henius plo s o u ile MgH2 wi h 1 hyd ogen acancy.
wi h coe icien s
𝑎𝑖=1
2∑
𝑗≠𝑖((1 − 𝑥𝑖)𝑓𝑗→𝑖+𝑥𝑖𝑓𝑖→𝑗)𝑟𝑖𝑗 ⊗ 𝑟𝑖𝑗 ,(55a)
𝑏𝑖=∑
𝑗≠𝑖((1 − 𝑥𝑖)𝑓𝑗→𝑖+𝑥𝑖𝑓𝑖→𝑗)𝑟𝑖𝑗 ,(55b)
𝑐𝑖=∑
𝑗≠𝑖
𝑥𝑖(1 − 𝑥𝑖) (𝑓𝑖→𝑗−𝑓𝑗→𝑖),(55c)
which depend on he local s a e and a omic con igu a ion.
The e ms in (54) accoun o mass apping, bias and di usion. In
pa icula , 𝑐𝑖 is he local a e o apping, 𝑏𝑖 is a local d i eloci y and
𝑎𝑖 is he local di usi i y enso . We no e ha 𝑎𝑖 is symme ic and can be
aniso opic in gene al. Re e ence di usi i y p ope ies can be de ined
by conside ing he di usion o an ini ially ully occupied and isola ed
si e, co esponding o 𝑥𝑖= 1 and 𝑥𝑗= 0 o 𝑗≠𝑖. In his case, (55)
educes o
𝑎𝑖=1
2∑
𝑗≠𝑖
𝑓𝑖→𝑗𝑟𝑖𝑗 ⊗ 𝑟𝑖𝑗 ,(56a)
𝑏𝑖=∑
𝑗≠𝑖
𝑓𝑖→𝑗𝑟𝑖𝑗 , 𝑐𝑖= 0.(56b)
Mechanics o Ma e ials 207 (2025) 105380
7
M. Molinos e al.
Fig. 9. hcp Mg a 600 K. Time e olu ion o a omic mola ac ions o 1 hyd ogen loca ed a he cen e o he pe iodic cell. (a) 0.01 ns; (b) 0.03 ns; (c) 0.06 ns. Magnesium
si es in g ay, a omic mola ac ions o hyd ogen shown colo -coded. (Fo in e p e a ion o he e e ences o colo in his igu e legend, he eade is e e ed o he web e sion
o his a icle.)
Fig. 10. hcp Mg a 600 K. Time e olu ion o a omic mola ac ions o 5 hyd ogens loca ed a andom loca ions in he pe iodic cell. (a) 0.01 ns; (b) 0.03 ns; (c) 0.06 ns.
Magnesium si es in g ay, a omic mola ac ions o hyd ogen shown colo -coded. (Fo in e p e a ion o he e e ences o colo in his igu e legend, he eade is e e ed o he
web e sion o his a icle.)
I all he si es a e embedded in he same local a omic con igu a ion,
as is he case, e. g., o a pe ec la ice, hen a single di usi i y enso
cha ac e izes he sys em. We also no e ha di usion is unbiased a si e
𝑖 i
∑
𝑗≠𝑖
𝑓𝑖→𝑗𝑟𝑖𝑗 = 0,(57)
which is ensu ed, e. g., by cen osymme y o he c ys al la ice.
We obse e ha o simple la ices hese e e ence p ope ies a e
independen o 𝑖, as expec ed om ansla ion in a iance, and depend
only on he local la ice s uc u e and he anspo coe icien s (39).
Fo complex la ices, he e ec i e di usi i y 𝑎𝑖 may a y om si e o
si e in acco dance wi h he local a omic con igu a ion o he si es.
In Fig. 7, we alida e p edic ions o he case o an isola ed hyd o-
gen a om in a Mg hcp cell agains expe imen ally measu ed a e aged
bulk di usi i ies (Nishimu a e al., 1999b), Renne and G abke (1978)
and DFT calcula ions (Klyukin e al., 2015). We obse e om he
igu e ha he di usion coe icien s a e in good ag eemen wi h he
expe imen al da a and co ec ly exhibi he expec ed hexagonal sym-
me y o hcp Mg. Fig. 8 shows co esponding di usion cons an s o
an isola ed acancy in u ile MgH2. He e again, he aniso opy o he
di usion enso is e iden in he igu e. In pa icula , bo h he hcp Mg-
H and e agonal u ile phase o MgH2, a e p edic ed o exhibi a ying
di usi i ies pa allel and no mal o he basal plane. Rema kably, he
p edic ed hyd ogen di usi i ies in magnesium hyd ide a e much lowe
han in magnesium, in ag eemen wi h expe imen al da a and he MD
calcula ions o Spa a u e al. (2020).
4.3. Hyd ogena ion and dehyd ogena ion kine ics
In o de o showcase and assess he abili y and e iciency o he
heo y o p edic mic os uc u al e olu ion, especially o long imes,
we conside p ocesses in ol ing hyd ogena ion and dehyd ogena ion
kine ics a a ious s ages o e olu ion o bo h magnesium and mag-
nesium hyd ide as a unc ion o s oichiome y and empe a u e. Fo
pu poses o compa ison, we speci ically choose con igu a ions iden ical
o hose analyzed by Spa a u e al. (2020) using MD simula ions.
The chosen con igu a ions a e conce ned wi h dilu e concen a ions
o hyd ogen in hcp Mg, and wi h dilu e concen a ions o hyd ogen
acancies in magnesium hyd ide. Bo h con igu a ions a e encoun e ed
in p ac ice du ing he ope a ion cycle o hyd ogen s o age ma e ials.
The wo k o Spa a u e al. (2020) can be consul ed o an in-dep h
discussion o he expe imen ally obse ed beha io and i s physical
o igins. In he p esen wo k, we ocus speci ically on assessing he
e iciency o DMD ela i e o MD, especially as ega ds long- e m
beha io .
Fo each a ge empe a u e and p essu e, we ini ially equilib a e
he sys em unde NPT condi ions and cons an ini ial hyd ogen occu-
pancies, and subsequen ly swi ch o 𝜇VT condi ions. We hen upda e
he sys em in ime by means o a s agge ed app oach in which: he
mean ield he momechanical p oblem (27) is i s sol ed a cons an
hyd ogen occupancy using he c i ical poin line sea ch (B une e al.,
2015) implemen ed in he PETSc/TAO lib a y (Balay e al., 2023);
and he mass anspo Eq. (41) is hen in eg a ed a cons an a omic
posi ions using a backwa d-Eule implici in eg a ion scheme. The ime
s ep is selec ed so as o esol e he mos es ic i e elaxa ion ime
deduced om he luxes 𝑓𝑖→𝑗 in (41).
4.3.1. Sho - e m analysis
We begin by conside ing he sho - e m hyd ogena ion o hcp Mg.
The compu a ional pe iodic cell con ains 42 (11
20) planes in he 𝑥
di ec ion, 24 (1
100) planes in he 𝑦 di ec ion, which coincides wi h he
𝑐 axis, and 12 (0001) planes in he 𝑧 di ec ion. The dimension o he
esul ing compu a ional cell is app oxima ely 67 × 66 × 62 Å, con-
aining 12096 Mg a oms. We conside h ee dilu e ini ial dis ibu ions
o hyd ogen, consis ing o 1, 5 and 121 occupied in e s i ials in he
Mechanics o Ma e ials 207 (2025) 105380
8
M. Molinos e al.
Fig. 11. hcp Mg a 600 K. Time e olu ion o a omic mola ac ions o 121 hyd ogens loca ed a andom loca ions in he pe iodic cell. (a) 0.01 ns; (b) 0.03 ns; (c) 0.06 ns.
Magnesium si es in g ay, a omic mola ac ions o hyd ogen shown colo -coded. (Fo in e p e a ion o he e e ences o colo in his igu e legend, he eade is e e ed o he
web e sion o his a icle.)
pe iodic compu a ional cell and co esponding oughly o s oichiome-
ies X𝐻= 8.27 × 10−5, 4.13 × 10−4 and 0.01, espec i ely. In o de
o asce ain he empe a u e dependence, we epea he calcula ions
o se en empe a u es in he ange o 400 o 700 K; and he u ile
calcula ions o nine empe a u es in he ange o 600 o 800 K.
Sequences o snapsho s o he ea ly s ages o he e olu ion o he
Mg-H sys em a 600K a e shown in Fig. 9, 10 and 11 a imes 0.01,
0.03 and 0.06 ns. The ime s ep used in he mass- anspo equa ion
calcula ions is 0.01 ns, which is 2 × 104 la ge han he ime s ep o
0.5 s used in MD calcula ions by Spa a u e al. (2020), o a ou o de -
o -magni ude gain. Thus, he DMD calcula ions ope a e on he di usi e
ime scale, in sha p con as o he ime s eps necessi a ed by molecula
dynamics which a e con olled by he pe iod o he mal ib a ion o he
la ice.
I bea s emphasis ha his s agge ing gain is achie ed a no loss
o ideli y. Indeed, he DMD calcula ions exhibi he expec ed phe-
nomenology o mass anspo on a la ice, which akes he o m o
aniso opic disc e e di usion. Thus, in he ea ly s ages o di usion, he
ini ial hyd ogen in e s i ial occupancies de ocus in o de ocus a omic
mola ac ion dis ibu ions o s eadily inc easing bu ini e ex en , in
con as wi h classical di usion in which a dis u bance o change in
concen a ion sp eads ins an ly h oughou he en i e domain, albei
wi h Gaussian decay.
Fo su icien ly sho imes, he a omic mola ac ion dis ibu ions
induced by each o he ini ial hyd ogen in e s i ials do no o e lap, bu
ne e heless in e ac weakly h ough hei elas ic ields, Figs. 9and 10.
The p edic ed aniso opy o hyd ogen di usi i y is also e iden om
he elonga ed shape o he e ol ing a omic mola ac ion dis ibu ions,
which e lec s [0 0 0 1] as he p e e en ial di ec ion o di usion, see
Fig. 7.
Fo a su icien ly dense ini ial dis ibu ion o hyd ogen in e s i ials,
o o su icien ly long imes, he a omic mola ac ion dis ibu ions
induced by each o he ini ial hyd ogen in e s i ials o e lap and in e ac
s ongly, bo h elas ically and en opically. Fig. 11. The esul ing e olu-
ion can be cha ac e ized by means o a clus e analysis o he a omic
mola ac ion dis ibu ions. Fig. 12. As expec ed om he disc e eness
o he p ocess, he e olu ion o he size o he a omic mola ac ion
clus e size de ia es om classical g ow h 𝑟∼√𝐷𝑡, and ins ead ini ially
exhibi s g ow h 𝑟∼𝑡3∕4 ollowed by ballis ic-di usion g ow h 𝑟∼𝑡 a
in e media e imes.
4.3.2. Long- e m dehyd ogena ion
Finally, we assess he abili y o DMD o cha ac e ize long- e m
kine ics. To his end, we conside he p ocess o dehyd ogena ion o
u ile MgH𝑋. Assuming a e agonal la ice alignmen , he esul ing
compu a ional cell o he u ile s uc u e o MgH2 con ains 16 (010)
planes in he 𝑥 di ec ion, 10 (011) planes in he 𝑦 di ec ion, and 14
(0
11) planes in he 𝑧 di ec ion. The dimension o he esul ing pe iodic
cell is 50 × 30 × 45 Å, con aining 2240 Mg a oms and 4480 H a oms.
Fig. 12. Time dependence o he s anda d-de ia ion adius o hyd ogen densi y a ound
ini ially occupied si es and in e ed ime exponen s showing non-classical di usion.
We conside wo hyd ogen s oichiome ies: X𝐻= 1.95 and 1.7. These
s oichiome ies a e ini ially ealized by andomly se ing hyd ogen
occupancies a in e s i ial si es o 0 (unoccupied) o 1 (occupied).
Speci ically, we ill 112 and 672 in e s i ial si es, espec i ely.
As al eady no ed in Sec ion 4.2, acancy di usi i ies in u ile MgH𝑋
a e much smalle han hyd ogen di usi i ies in hcp Mg, esul ing in
a much slowe e olu ion o he acancy dis ibu ion and, co espond-
ingly, in he need o ex end he analysis o much longe imes. I
is, he e o e, c ucial ha he ime s ep equi ed in calcula ions o
in eg a e he mass- anspo Eq. (41) is 1 ns. Again, i bea s emphasis
ha his ime s ep is chosen o esol e he di usi e ime scale and is a
s agge ing six o de s o magni ude la ge han he ime s eps equi ed
by molecula dynamics (c ., e. g., Spa a u e al. (2020)), which a e in
he s scale.
Fig. 13, 14 and 15 show he e olu ion o he acancy dis ibu ion
up o 10 ns a 𝑇= 600 K o he sys ems o 1, 112 and 672 hyd ogen
acancies, espec i ely. As in he hcp Mg hyd ogena ion case, he
e olu ion o he 1- acancy sys em e lec s clea ly he aniso opy o
he di usi i y enso , see Fig. 8. The de ocusing o he dis ibu ions
o he indi idual acancies is localized a ea ly imes, as expec ed
om disc e e di usion, and e ol e independen ly. A longe imes,
he dis ibu ions o he indi idual acancies become delocalized and
in e ac s ongly, wi h a gene al end owa ds a uni o m dis ibu ion,
as expec ed.
The s uc u al s abili y o he sys ems is showcased in Figs. 16 and
17, which show [
110]- iews o con igu a ions a imes 1 ns and 10 ns
Mechanics o Ma e ials 207 (2025) 105380
9