21 Feb 2024

RIKEN Symposium “Fundamentals in Density Functional Theory”
20-22 Feb 2024, RIKEN Kobe Campus

Glass Transition Studied by
Molecular Dynamics
Simulations
K. SHIRAI1,2), K. Watanabe2), H. Momida2),and S. Hyun3)
1 VNU Vietnam Japan University, Hanoi, Vietnam
2 SANKEN, Osaka University, Ibaraki, Japan
3 Korea Institute of Ceramic Engineering and Technology, Jinju-si, Koria

Glasses are nonequilibrium states
History dependence(a)

VFT law(Vogel-Fulcher-Tammann)?
13U ≠ U(T, P)η = expLogηbut

True thermodynamictransition?

U = U(T, P, {history})
Tg

T0• Memory effect• Aging phenomena(b)glassCl• Rejuvenation

Ttransition
regionCD( T − T0 )liquidKinetic transition ?
Cg
Tg,1in DFT

U = Ē(T, P, {R̄j})

Tg,2

T(c)

(l)
(sl)

S(g)

K.S, J. Phys. Commun. 4, 085015 (2020)
J. Phys. Commun. 5, 015004 (2021)

(c)
TK

Tm

Nature of Glass Transition
glasstransition
regionsupercooled
liquidnormalliquid∆Tg

Cp

ΔCp = Cp(l) − Cp(g)Cex

Excess speci c heat

crystalTg

TmΠ : Prigogine-Defay (PD) ratio
Π=1∆Cpfi

no rst-order transition

Jump at Tg∆Cpfi

Absence of latent heat L

Kinetic origin?

(l)(c)
Cex = Cp − CpTΠ=ΔCpΔκTV(Δα)2>1

for the second-order transitionDif culties in calculating ∆Cp
1. Lack of standard theory of speci c heat for liquid(g)(l)ΔCp = Cp − Cp

2. Hysteresis of C-T curve in glass transition

3. Large separation in the timescale
between glass transition and MD simulations
Glass transition time10 psfifi> 1 minsimulation time of FP-MDTheory of speci c heatSpeci c heat of silica glass(experiment)

Cv: Classical limit

32gas?

liquid3solid4RichetSudo

Cp /R3in kB

21per atom

Tst

TanTsf00500

1000

1500

2000

T (K)

P. Richet, et al., Geochim. Cosmochim. Acta 46 2639 (1982)

fifi

H. Sudo, Practical Manual for Amorphous Siliceous Materials,
ed H. Kawazoe (Tokyo: REALIZE, 1999), p. 83.

First-principles calculation of speci c heatSpecific heat

∂UCx =( ∂T )Xhigh T

Microscopic energy

1Etot(t) = EGS({Rj(t)}) +Mjvj(t)22∑
jEKE phlow T
EEP,vibin thermodynamic equilibrium

Internal energy
1
U = Etot(t) = EGS({Rj(t)}) +Mjvj(t)22∑jHarmonic approx.

Uin adiabatic MDRj(t) = R̄j + uj(t)Est({R̄j}) + Eph(T ) + Ete(V )structural part(configurational)

Cp=

phonon

Cst

Est

X

thermal expansion+ Cph+ Cte
fiFor solids,

3kBT
2

Phonon contribution to speci c heatFor solids
Eph =∑(n̄k +kPhonon DOS of quartz1ℏωk)
20.20ℏωe ℏω/kBTCph(T ) = kBg(ω)dω∫ ( kBT ) (e ℏω/kBT − 1)22

0.15

0.10

0.05

Cp = Cst + Cph + Cte

0.00

structural partphonon0200

400600ω

thermalexpansion

800

1000

1200

1400(cm-1)

For liquids
Is the phonon-picture realistic for liquids?

U = Etot(t)No!
Effect of energy dissipationfi

Need of adiabatic MD simulation

Problem for liquids
No eigenstatesZ=∑eEnergy dissipation determines C−βϵi

Relaxation processesidU(ω) = C(ω)dT(ω)ω1Destruction
ω2Creation

ω3

ω4N. O. Birge and S. R. Nagel, Phys. Rev. Lett.
54, 2674 (1985).

C0C(ω) =1 − iωτs

Adiabatic MD simulations automatically takethe energy dissipation process into account.

Tini

Energy dissipation
Relaxation
Relaxation
U=constT
Equilibrium

Glass Transition of Glycerol
0Exp.

2Log D (Å /ps)

Diffusion
constant

Tg = 185 K-1∆Cp = 0.70 R-2400
Heating
Cooling300

Reheating

E st (K)

Structuralenergy200

0.70

ΔCst = 0.58 R

0.1210000200≈ ΔCp(exp)TmTg

400600800

1000

T (K)

K. Shirai, K. Watanabe, and H. Momida, J. Phys.: Condens. Matters 34, 375902 (2022)

Glass Transition of Silica

2Diffusion
constantLog D (Å /ps)

1

Coolingprocess

01Exp.
fastTg = 1480 K(A)(B) slow2∆Cp = 0.32 R

3Structuralenergy

E st (K)

6000

4000

ΔCst = 0.5 R

2000

001000

2000

3000

4000

5000

T (K)

K. Shirai, et al., J. Phys.: CM 35 505401 (2023)

Nature of Glass Transition
ΔCp ≈ ΔCst

For Silica glass
more than 99% of ∆Cp comes from ∆Cst

Almost all the part ∆Cp is determined bythe structural change.
ΔCp(Tm − Tg) ≈ Lessentially the structural transition

fibut the randomness breaks the latent
heat into the excess speci c heat.

LCex
∆CpTg

TmT

Interpretation of Prigogine-Defay ratio
Π=ΔCpΔκ

Π>1TV(Δα)2There are more than one order parameters

TV 2Cte =ακΔκ Δα=
κα

Π=ΔCp
ΔCte=

Π>1

Nat{R̄j}

ΔCp = ΔCst + ΔCph + ΔCte
(Change in the total energy)(Part of isotropic volume change)

Structural transition

Summary

Glass state

Present results

Traditional view

Equilibrium state

Nonequilibrium state

U = U(T, V, {R̄j})∆Cp

Thermodynamicallywell-de ned quantity

U ≠ U(T, V )Kinetic parameter?

mostly determinedby ∆Est

Π>1fiGlasstransition

Structural transition

due to order parameters

Kinetic transition

The present view for glass
The traditional view

Experimental fact

Glass has hysteresis.

All solids have hysteresis.

NEQ

U = U(T, {R̄j})

Present theory

All the static states of a solid are equilibrium statesfrozenglass

General theory of state variables

arXiv:1812.08977

Thermodynamic method for glass

J. Phys. Commun. 4, 085015 (2020)
5, 015004 (2021)Activation energy of GT

J. Phys. Commun. 5 095013 (2021)

FP-MD simulation of speci c heat jump
FP-MD study on the Prigogine-Defay ratio

J. Phys.: CM 34 375902 (2022)

The third law

arXiv:2207.11421

fiApplications

within the time scale in which the solid maintains its current structure.

J. Phys.: CM 35 505401 (2023)

History-dependent properties
Nonequilibrium!

• Memory effect
The Kovacs effectU ≠ U(T, P)
but

U = U(T, P, {history})

Memory andrejuvenation inpolymethyl
methacrylate

K. Fukao and A. Sakamoto,Phys. Rev. E 71, 041803 (2005)• Aging phenomena

in DFTU = Ē(T, P, {R̄j})K.S, J. Phys. Commun. 4, 085015 (2020)• Rejuvenation

J. Phys. Commun. 5, 015004 (2021)

Energy landscape of spin glass
Conventional spin glass
Spin JamhierarchicalnonhierarchicalA. M. Samarakoon, et al., Proc. Natl.
Acad. Sci. U.S.A. 113 11806 (2016)