Cooling-rate dependence of kinetic and mechanical stability of simulated glasses

0.35 MB·English

by  Hannah Staley

Cooling-rate dependence of kinetic and mechanical stability of simulated glasses
Hannah Staley,1 Elijah Flenner,2 and Grzegorz Szamel2
1Department of Physics, Colorado State University, Fort Collins, Colorado 80523, USA
2Department of Chemistry, Colorado State University, Fort Collins, Colorado 80523, USA
(Dated: January 15, 2015)
Recently,ultrastableglasseshavebeencreatedthroughvapordeposition. Subsequently,computer
simulation algorithms have been proposed that mimic the vapor deposition process and result in
simulated glasses with increased stability. In addition, random pinning has been used to generate
very stable glassy configurations without the need for lengthy annealing or special algorithms in-
spired by vapor deposition. Kinetic and mechanical stability of experimental ultrastable glasses is
compared to those of experimental glasses formed by cooling. We provide the basis for a similar
5 comparison for simulated stable glasses: we analyze the kinetic and mechanical stability of simu-
1 latedglassesformedbycoolingataconstantratebyexaminingthetransformationtimetoaliquid
0 upon rapid re-heating, the inherent structure energies, and the shear modulus. The kinetic and
2 structuralstabilityincreasesslowlywithdecreasingcoolingrate. Themethodsoutlinedherecanbe
n used to assess kinetic and mechanical stability of simulated glasses generated by using specialized
a algorithms.
J
4 PACSnumbers: 64.70.Q-,81.05.Kf
1

] Vapordepositiononasubstrateheldataround85%of and de Pablo [6] and Lyubimov, Ediger, and de Pablo
t the glass transition temperature is used to create glasses [7] modeled the vapor deposition of a binary mixture of
f
o that have a higher kinetic stability [1] and larger elastic Lennard-Jonesparticlestosimulatethecreationofanul-
s
moduli [2] than glasses created by cooling from a liq- trastableglass. Theyfoundanenhancedkineticstability
.
t
a uid. Due to these and other desirable properties of these compared to the glass obtained by cooling from a liquid
m so-called ultrastable glasses, it is of great interest to un- at a constant rate, but the enhanced stability was less
- derstand the differences between glasses created through pronouncedthaninexperiments. Furthermore,thesimu-
d vapor deposition and by cooling from a liquid. Since va- latedvapordepositionresultedinglasseswithstructural
n
por deposited glasses are created one layer at a time, it anisotropy, a higher density, and non-uniform composi-
o
issuspectedthatenhancedmobilityatthesurfaceallows tion compared to the simulated glasses created by cool-
c
[ the molecules to find more stable configurations. ing from a liquid. Very recently, Hocky et al. [8] created

two dimensional binary Lennard-Jones mixture glasses
1 Discoveryofultrastableglassesinthelabinspiredsev-
by randomly pinning particles. Glasses formed by ran-
v eralcomputersimulationstudies. Thegoalofthesestud-
9 dom pinning are, by construction, isotropic and uniform
ies has been two-fold. First, modeling the vapor depo-
4 (on average). Hocky et al. found that these glasses have
sition process can provide insight into the origin of the
4 increasedkineticstabilitybutthepresenceofpinnedpar-
3 stability of vapor deposited glasses. Second, the experi-
ticles prevented them from investigating whether these
0 mentalprocedurecanbeusedasaninspirationforthede-
glasses also have a larger shear modulus.
1. velopment of computer simulation algorithms that could
0 generate very stable simulated glasses without the need The above mentioned simulational studies focused on
5 for lengthy annealing. L´eonard and Harrowell [3] mod- the properties of stable glasses created by using different
1
eledultrastableglassesusingathree-spinfacilitatedIsing specialized algorithms. However, it is difficult to assess
:
v model. Theirstudyprovidedevidencethatenhancedsur- the stability of these glasses without knowledge of the
i facemobilityduringvapordepositionresultsinimproved stability of simulated glasses generated in a conventional
X
kinetic stability, but the simplicity of their model pre- way, i.e. by cooling at a constant rate. We emphasize
r
a ventedthemfromstudyingotherpropertiesofultrastable that simulated glasses generated by cooling are isotropic
glasses, such as the shear modulus. On the other end of anduniform. Thus,studyingtheseglassesmayalsoshed
the complexity spectrum, Singh et al. [4] modeled vapor somelightontheimportanceoftheanisotropyandcom-
depositionoftheorganicmoleculetrehalosebyintroduc- positionalnon-uniformityfortheglassstability. Herewe
ing1-5“hot”moleculesaboveasubstrateandthen(arti- createglassesbycoolingfromthesupercooledliquidand
ficially)coolingthemtothesubstratetemperature. This weassesstheirkineticandmechanicalstabilitybyexam-
procedure resulted in a simulated glass with a higher ki- ining their kinetic stability against reheating, properties
neticstability,butalsopronouncedstructuralanisotropy oftheirpotentialenergylandscape,andtheirshearmod-
in the direction perpendicular to the substrate. ulus.
More recent studies used well-known glass-forming bi- We simulated the Kob-Andersen (KA) [5] 80:20
nary Lennard-Jones mixtures [5]. First, Singh, Ediger, binary Lennard-Jones mixture in three dimensions.
2
The particles interact via the potential V (r) = 1
αβ
4(cid:15) (cid:104)(cid:0)σαβ(cid:1)12−(cid:0)σαβ(cid:1)6(cid:105), and the parameters are (cid:15) =
αβ r r BB 0.8
0.5(cid:15) , (cid:15) = 1.5(cid:15) , σ = 0.88σ , and σ =
AA AB AA BB AA AB
0.8σ . The majority species are of type A and both
AA 0) 0.6
masses are equal. We cut the potential at 2.5σαβ. The =w
raapentasudrurtlei(cid:112)tc,sleamasnrAaedtσptAa2riemAns/uee(cid:15)n,mAtrbeAedesrpbineedcierntneigvdsieuttlhycyee.NduW/unVneitits=ssifmwo1r.ui2tl0hlaetnuσegsdAitAnhNg,,(cid:15)ptA=eeAmri8/op0kde0Bri0c-, q(t,ts00..24 2〈δ〉r(t)11110000-0211
10-2
boundary conditions. 100 102 104
t
0
WeranNVTsimulationsusingaNos´e-Hooverthermo- 0 1 2 3 4
10 10 10 10 10
stat. We started with an equilibrated supercooled liquid t
at the temperature T = 0.5 (for this system the onset
temperature for the slow dynamics is To ≈ 1.0 and the FIG. 1: Average overlap function qs(t,tw = 0) monitored
during the heating trajectories (solid lines) compared to the
mode-coupling temperature is T = 0.435). We ran 80
c
cooling runs at cooling rates T˙ =∆T/∆t of 3.33×10−n averageoverlapfunctionfortheequilibriumliquidatT =0.5
(dashedline). Solidlinescorrespondtoglassesobtainedwith
where n = 3, 4, 5, 6, and 7. Each run started at an cooling rates T˙ = 3.33 × 10−n, where n = 3, 4, 5, 6, 7,
equilibrium configuration at T = 0.5 and was cooled to
and 8 from left to right. The inset shows the mean square
T = 0.3. We also ran 4 cooling runs at a cooling rate displacement versus time for the same cooling rates (solid
of 3.33×10−8. For T˙ > 3.33×10−8, we started from lines)andfortheequilibriumliquidatT =0.5(dashedline).
the final configuration of each of the 80 cooling runs and
annealedthesystemforatleast100timesτ determined
α
from T = 0.5 (we define τ through the average over- itythataparticlemovedoverdistanceaduringthetime
α
lap function discussed later). For the cooling rate of between t and t+t , where t is the waiting time.
w w w
3.33×10−8, we performed 5 different annealing runs for
(cid:42) (cid:43)
each of the 4 configurations, each beginning with a dif- 1 (cid:88)
q (t,t )= q (t,t ) , (1)
ferent initial random velocity. In addition, we heated s w N m w
m
the final configurations of the cooling runs by ramping
up the temperature from T = 0.3 to T = 0.5 at a con- where q (t,t ) = Θ(a−|r (t+t )−r (t )|), Θ is
m w m w m w
stant rate over a time t = 10 (we tried increasing the Heaviside’s step function, and r (t) is the position of
m
temperature instantaneously from T = 0.3 to T = 0.5 a particle m at a time t. t , the waiting time, is mea-
w
but found that this resulted in large, non-physical, os- sured from the start of the heating trajectory. We chose
cillations of the kinetic and potential energies). After a = 0.25 to be consistent with previous work [16]. We
ramping up the temperature to T = 0.5, we continued define τ as when q (τ ,t ) = e−1 for an equilibrium
α s α w
constant temperature simulations (at T = 0.5) until the system (note that for an equilibrium system q (t,t )
s w
mean-square displacement started growing linearly with does not depend on the waiting time, i.e. it is time-
time. We refer to this ramping up of temperature and translationally invariant).
the subsequent run at T = 0.5 as a heating trajectory.
Fig. 1 shows average overlap function q (t,t =0) ob-
For T˙ > 3.33×10−8, we ran 80 heating trajectories, in s w
tainedwhileheatingglassespreparedatdifferentcooling
eachcasestartingfromthefinalconfigurationofadiffer-
rates (solid lines) and obtained from the equilibrium run
entcoolingrun. ForT˙ =3.33×10−8,werantheheating
at T = 0.5 (dashed line). The feature at t = 10 is a
trajectories15timesusingdifferentinitialrandomveloc-
consequence of changing the thermostat from a ramp-up
ities for each of the configurations obtained from the 4
of the temperature at a constant rate from T = 0.3 to
cooling runs. Finally, we also performed an equilibrium
T = 0.5 to maintaining constant temperature T = 0.5.
run at T = 0.5, for at least 100τ . We used HOOMD-
α For t > 10, q (t,t = 0) exhibits prolonged plateaus for
s w
blue [9, 10] for equilibration at T = 0.5 and the cooling
the three slowest cooling rates followed by a very rapid
simulations. The heating trajectories, the T = 0.3 runs,
decayatlatertimes(theglassgeneratedusingtheslowest
andthesinglerunatT =0.5usedLAMMPS[11,12]run cooling rate,3.33×10−8, relaxes faster than for ballistic
on a GPU [13–15].
motion). With decreasing cooling rate both the plateau
We start with an assessment of the kinetic stability of height and its extent increase. This indicates that with
the simulated glasses. First, we we examine the time it decreasing cooling rate the cage diameter decreases and
takes for the particles to rearrange after the sudden in- ittakeslongerfortheparticlestomoveoutoftheircages.
creaseintemperature,i.e. alongtheheatingtrajectories. This is the first indication of increasing kinetic stability
We quantify these rearrangements through the average against re-heating.
overlap function q (t,t ), which measures the probabil- The inset in Fig. 1 shows the mean square displace-
s w
3
2 -7.78
10
-8.22
-7.8 S-8.26
EI
-8.3
1
-7.82
S101 0.8 〉U 10-8 10-6 • 10-4
)w0.6 〈 T
q(t,ts0.4 -7.84
0.2
0
100 101 102 103 104 -7.86
t
0
1010-8 10-7 10-6 10-5 10-4 10-3 10-2 10-8 10-7 10-6 10-5 10-4 10-3
• •
T T
FIG.2: DependenceofstabilityratioS =t /τ oncooling FIG. 3: Average potential per particle, (cid:104)U(cid:105), as a function of
trans α
rate T˙. The inset shows the average overlap function versus cooling rate, T˙. The inset shows average inherent structure
time for the cooling rate T˙ =3.33×10−6 and waiting times, energy per particle, EIS, as a function of cooling rate. The
from top to bottom, 0, 10, 500, 1000, 3000, 4250, 5000. calculationsweredoneatcoolingrates3.33×10−n,withn=
8,7,6,5,4,3, the lines are to guide the eye.
ment (cid:10)δr2(t,t )(cid:11) = N−1(cid:68)(cid:80) [r (t+t )−r (t )]2(cid:69)
w n n w n w
bility ratios of the most stable simulated glasses, which
for a waiting time t =0 obtained while heating glasses
w
approach approximately 400 [8], we would need to de-
prepared at different cooling rates (solid lines) and ob-
crease the cooling rate by about 3 decades (to achieve
tainedfromtheequilibriumrunatT =0.5(dashedline).
We find that (cid:10)δr2(t,t )(cid:11) encodes similar information as stabilityratiosofexperimentalultrastableglasses,which
w
approach approximately 103.5 [17] we would need to de-
q (t,t = 0). In particular, we note that with decreas-
s w
crease the cooling rate by about 22 decades). We con-
ing cooling rate the time required for the mean-square
cludethatevensimulatedstableglasseswouldbeimpos-
displacement to reach the equilibrium curve increases.
sibletoobtainbycooling,atleastwiththecomputational
This is another indication of increasing kinetic stability
resources available to us at present.
against re-heating.
To quantify kinetic stability we follow earlier experi- Next, we turn to an assessment of the mechanical sta-
mental [17] and simulational [8] studies and we evaluate bility of the simulated glasses. The simplest, albeit in-
the transformation time t and the stability ratio S. direct, measure of the mechanical stability is provided
trans
The transformation time is defined as the time it takes by the inherent structure energy [18], since in the en-
for the liquid obtained by heating the glass to return ergy landscape picture of the glass transition, more sta-
to equilibrium. In the experimental study of Ref. [17] ble glasses are in a lower energy basin of the landscape
the transformation time was obtained by monitoring the and have to overcome a larger energy barrier to deform.
time evolution of the dielectric response. We follow the We show the dependence of the average inherent struc-
simulationalstudyofRef. [8]andobtainthetransforma- ture energy (obtained using the Fire algorithm [19] in
tion time from monitoring the waiting time dependence HOOMD-blue) at T = 0.3 on the cooling rate used to
of the average overlap function. Specifically, we calcu- reachthistemperatureintheinsetinFig.3. Weseethat
late q (t,t ) for a range of t (see the inset to Fig. 2 for the average inherent structure energy decreases with de-
s w w
q (t,t ) for various waiting times for T˙ = 3.33×10−6) creasing cooling rate. In Fig. 3 we show the dependence
s w
andwedefinethewaitingtime-dependentrelaxationtime oftheaveragepotentialenergyperparticleatT =0.3on
τ (t ) through the equation q (τ ,t ) = e−1. We define the cooling rate. It also decreases with decreasing cool-
s w s s w
the transformation time, t , as the waiting time such ing rate. We note that for their vapor deposited glass at
trans
that τ (t ) = τ , where τ is the equilibrium relaxation T = 0.3, Ediger and de Pablo got an inherent structure
s w α α
time. Finally, we follow Refs. [8, 17] and we define a energy of -8.35 [7] and a potential energy of -7.8 (esti-
stability ratio S =t /τ . mated from [7]). These values are close to our values of
trans α
In Fig. 2 we show the cooling rate dependence of the -8.30 and -7.86. We note, however, that the vapor de-
stability ratio. At larger cooling rates S increases rather position inspired algorithm used in Ref. [7] resulted in
quicklywithdecreasingT˙. However,forcoolingratesless glasses that were non-uniform in contrast to our simu-
than 3.33×10−5 we get a slower increase of the stability lated glasses, and thus, a literal comparison of the two
ratio with cooling rate. For the slowest cooling rates, an studies is impossible.
order of magnitude decrease in the cooling rate results A more direct measure of mechanical stability is pro-
in a doubling of the stability ratio. Thus, to achieve sta- vided by elastic constants. Ultrastable glasses prepared
4
in experiments [2] and stable glasses obtained in simula- 19
tions using vapor deposition-inspired algorithms [6] have
largerelasticconstantsthanglassespreparedbyordinary 18
cooling. Weexamineherethecoolingratedependenceof
twheeussheeaanr mapopdruolaucshorfecoeunrtlsyiminutlraotdeducgeldasbsyest.wTooofthuiss[e2n0d]. µ 17 -1q;t)]
This method is based on the relationship between corre- 16 S(4
lations of transverse particle displacements in the glass 2q
ρ[ 101
and the glass’s shear modulus. The final result is that T
15 kB
the shear modulus of a glass can be measured by study- 2 100
q
ingthesmallwave-vectorq behavioratlongtimesofthe
14
four-point structure factor -8 -7 -6 -5 -4
10 10 10 10 10
•
T
(cid:42) (cid:43)
1 (cid:88)
S4(q;t)= N g[δrn(t)]g∗[δrm(t)]eiq·[rn(0)−rm(0)] , FIG. 4: Dependence of the shear modulus, µ, on the cool-
n,m ing rate T˙ for glasses at T = 0.3. The inset shows
(2) 2k Tρ(cid:2)q2S (q;t)(cid:3)−1 versus q plotted at a time t = 100 for
B 4
where the weighting function g[δrn(t)] = δr⊥n(t) = T =0.3,forthefourslowestcoolingrates: 3.33×10−n where
(r (t)−r (0))⊥ and δr⊥(t) is the component of n=5,6,7,8, from bottom to top.
n n n
the displacement of particle n perpendicular to the
wave-vector q, δr⊥(t) · q = 0. As argued in
n
Ref. [20], the shear modulus for a glass can be ob- relative increases of the shear modulus achieved in our
tained from the small q, long t behavior of S4(q;t), and Ashwin et al. simulational studies are comparable
µ = lim lim 2k Tρ(cid:2)q2S (q;t)(cid:3)−1. In practice to the relative differences between the shear modulus of
q→0 t→∞ B 4
one can examine the small wave-vector behavior of experimentalultrastableglassesandexperiemntalglasses
2k Tρ(cid:2)q2S (q;t)(cid:3)−1 for times longer than the initial generatedbycooling(19%forindomethacinand15%for
B 4
transient dynamics but shorter than time scale for the tris-naphtylbenzene [2]). It is unclear to us whether this
aging process of a glass. Also, in order to get good fact has some deeper meaning.
statistics one either needs many more independent con- In summary, we showed here that by decreasing the
figurationsattheendofeachcoolingprocessthanthe80 cooling rate in simulations, we generated simulated
that we generated or one needs to continue running at glasses that have increasing kinetic and mechanical sta-
T =0.3 and average over different time origins. In order bility. The highest kinetic stability that we achieved is
to measure properties of the glass as prepared using a about 66. Decreasing the cooling rate by one decade,
specific cooling rate, one can only average over different which would result in the kinetic stability of about 130,
time origins before the start of the aging process. For would require approximately 1 year on a single GPU us-
faster cooling rates the aging starts very quickly and the ing HOOMD-blue per one cooling trajectory. This re-
range of available time origins is rather small. Fortu- sult establishes a baseline against which specialized al-
nately, the present method to calculate the shear mod- gorithms that are inspired by vapor deposition or use
ulus does not require very long trajectories [20]. We random pinning should be measured.
show 2kBTρ(cid:2)q2S4(q;t)(cid:3)−1 versus q at a time t = 100 We believe that the kinetic stability, quantified by the
for T =0.3 in the inset to Fig. 4. We averaged over time stability ratio S =t /τ is the best quantity to com-
trans α
origins for t≤100 for the 3 faster cooling rates, and for pare stability of glasses generated using different proto-
t≤400 for the slowest cooling rate. cols and algorithms. We note that the average potential
In Fig. 4 we show the cooling rate dependence of the energy is very sensitive to the average density, compo-
shear modulus. We find that the modulus increases with sitional nonuniformities and possibly to anisotropy. In
decreasing cooling rate, although the rate of increase fact, our glasses generated by cooling at constant rate
seemstobedecreasingwithdecreasingT˙. Attheslowest have lower average potential energy than stable glasses
cooling rates an order of magnitude decrease of the cool- generated by Singh et al. [6] using a vapor deposition
ing rate results in an increase of the modulus by 4.6%. inspired algorithm. We note the average inherent struc-
We note that Ashwin, Bouchbinder, and Procaccia [21] ture energy seems to correlate better with the stability
alsoinvestigatedthecoolingratedependenceoftheshear but most likely it is also sensitive to the density and
modulus. Theyusedadifferentsystem,intwospatialdi- compositional nonuniformities. Finally, we note that the
mensions and thus their results cannot be literally com- relativechangeoftheshearmodulusofsimulatedglasses
pared to ours. However, we do note that in a similar formed by cooling at different cooling rates is compara-
range of cooling rates our shear modulus increases by ble to the relative difference between the shear modulus
15% and theirs increases by 12%. Interestingly, these ofexperimentalglassescreatedthroughvapordeposition
5
and by cooling. This surprising finding, which deserves [7] I.Lyubimov,M.D.Ediger,andJ.J.dePablo,J.Chem.
further investigation, suggests that that the kinetic and Phys. 139, 144505 (2013).
mechanical stability are not necessarily correlated. [8] G.M.Hocky,L.Berthier,andD.R.Reichman,J.Chem.
Phys. 141, 224503 (2014).
The methods presented here can be used to assess ki-
[9] http://codeblue.umich.edu/hoomd-blue.
netic and mechanical stability of simulated glasses gen-
[10] J.A.Anderson,C.D.Lorenz,andA.Travesset,J.Com-
erated by using specialized algorithms.
put. Phys. 227, 5342 (2008).
We gratefully acknowledge the support of NSF grant [11] http://lammps.sandia.gov.
CHE 1213401. [12] S. Plimpton, J. Comput. Phys. 117, 1 (1995).
[13] The GPU package was created by Mike Brown.
[14] W.M.Brown,P.Wang,S.J.Plimpton,andA.N.Thar-
rington, Comput. Phys. Commun. 182, 898 (2011).
[15] W. M. Brown, A. Kohlmeyer, S. J. Plimpton, and A. N.
[1] S.F.Swallen,K.L.Kearns,M.K.Mapes,Y.S.Kim,R. Tharrington,Comput.Phys.Commun.183,449(2012).
J.McMahon,M.D.Ediger,T.Wu,L.Yu,andS.Satija, [16] E. Flenner and G. Szamel, J. Phys.-Condens. Mat. 19,
Science 315, 353 (2007). 205125 (2007).
[2] K.L.Kearns,T.Still,G.Fytas,andM.D.Ediger,Adv. [17] A.Sepu´lveda,M.Tylinski,A.Guiseppi-Elie,R.Richert,
Mater. 22, 39 (2010). and M.D. Ediger, Phys. Rev. Lett. 113, 045901 (2014)
[3] S.L´eonardandP.Harrowell,J.Chem.Phys.133,244502 [18] F. Sciortino, J. Stat. Mech. P05015 (2005).
(2010). [19] E. Bitzek, P. Koskinen, F. Gahler, M. Moseler, and P.
[4] S.SinghandJ.J.dePablo,J.Chem.Phys.134,194903 Gumbsch, Phys. Rev. Lett. 97, 170201 (2006).
(2011). [20] E.FlennerandG.Szamel,Phys.Rev.Lett.114,025501
[5] W. Kob and H. C. Andersen, Phys. Rev. E 51, 4626 (2015).
(1995). [21] AshwinJ.,E.Bouchbinder,andI.Procaccia,Phys.Rev.
[6] S. Singh, M. D. Ediger, and J. J. de Pablo, Nat. Mater. E 87, 042310 (2013).
12, 139 (2013).