Transitions between Inherent Structures in Water

8 Pages·0.33 MB·English

by  Nicolas Giovambattista

Transitions between Inherent Structures in Water
Nicolas Giovambattista1, Francis W. Starr2, Francesco Sciortino3,
Sergey V. Buldyrev1, H. Eugene Stanley1
1Center for Polymer Studies and Department of Physics, Boston University,
Boston, Massachusetts 02215.
2Polymers Division and Center for Theoretical and Computational Materials Science, National Institute of Standards and
Technology, Gaithersburg, MD, 20899 USA
2 3 Dipartimento di Fisica, Istituto Nazionale per la Fisica della Materia,
0 and I.N.F.M. Center for Statistical Mechanics and Complexity,
0 Universit`a di Roma La Sapienza, P.le A. Moro 2, I-00185 Roma, ITALY
2 (EH8157. Time scale corrected in Figs 1,2,8: 31 Dec2001)

n
a
J
While such a picture of liquid dynamics is difficult to

3 The energy landscape approach has been useful to help verify experimentally, computer simulation offers an ex-
understandthedynamicpropertiesofsupercooledliquidsand
cellent opportunity to explore these ideas. For a pre-
] the connection between these properties and thermodynam-
t defined liquid potential, a liquid trajectory can be gen-
f ics. The analysis in numerical models of the inherent struc-
o erated via molecular dynamics simulation and the lo-
ture (IS) trajectories — the set of local minima visited by
s cal potential energy minima can be evaluated by an en-
the liquid — offers the possibility of filtering out the vibra-
at. tional component of the motion of the system on the poten- ergy minimization method [2]. With this procedure, the
motion in phase space is converted into a minimum-to-
m tial energy surface and thereby resolving the slow structural
minimum trajectory, or IS trajectory. A general picture
component more efficiently. Herewe report an analysis of an
- of the system moving among a set of basins surrounding
d IS trajectory for a widely-studied water model, focusing on
the multitude oflocalminima has evolved. More specifi-
n the changes in hydrogen bond connectivity that give rise to
o many IS separated by relatively small energy barriers. We cally, simulations have shownthat both the depth of the
c findthatwhilethesystemtravels throughtheseIS,thestruc- minima sampled by the system, as well as the number
[ ture of the bond network continuously modifies, exchanging of these minima, decrease on cooling [5,10,11]. Simula-
1 linearbondsforbifurcatedbondsandusuallyreversingtheex- tionshavealsoshownthatbelowacrossovertemperature
v changetoreturntonearlythesameinitialconfiguration. For T×, only rare fluctuations bring the system to a saddle
8 the216moleculesystemweinvestigate,thetimescaleofthese point and hence activated processes become important
2 transitions is as small as the simulation time scale (≈ 1 fs). for relaxation of the liquid [7,9,12–14]. For the system
0 Hence for water, the transitions between each of these IS is we study, T× coincides numerically [9,14–17] with the
1
relatively small and eventual relaxation of the system occurs critical temperature identified by mode coupling theory
0
only by many of these transitions. We find that during IS (MCT), which has been widely used to understand the
2
0 changes, themolecules with thegreatest displacements move dynamicsofliquidson“weak”supercooling(theT range
/ in small “clusters” of 1-10 molecules with displacements of wherecharacteristicrelaxationtimesapproach≈10−8 s)
t ≈ 0.02 −0.2 nm, not unlike simpler liquids. However, for
a [18].
m water these clusters appear to be somewhat more branched The description of the real motion of the system as
than thelinear “string-like” clusters formed in a supercooled
- an IS trajectory becomes a powerful way of separating
d Lennard-Jones system found by Glotzer and her collabora- the vibrational contribution, responsible for the thermal
n tors.
broadeningofinstantaneousmeasurementsfromtheslow
o
structural component[19]. Such an approach becomes
c
: even more powerful below T×, since most of the instan-
v
taneous configurations are far from saddles, making cor-
i
X relation functions calculated from the IS trajectory fully
I. INTRODUCTION
account for the α-relaxation dynamics[7].
r
a Here we study the IS trajectory below T× for the ex-
It is believed that the properties of liquids can be un- tendedsimplepointcharge(SPC/E)potential[20],awell-
derstood as motion of the system in a high-dimensional studied model for water. The dynamics of the SPC/E
complexpotentialenergysurface(PES)[1–6]. Asaliquid modelhavebeenshownto be consistentwith the predic-
is cooled toward the glassy state, the system is increas- tions of MCT [21,22]. Additionally, the PES of SPC/E
ingly found near local potential energy minima, called has been studied in detail above T× and a thermody-
inherent structure (IS) configurations [2]. As tempera- namic descriptionof the supercooledstates basedon the
ture decreases, the description of the dynamics in terms PEShasbeenpresented[11,23]. Herewefocusonthege-
of motion on the PES becomes increasingly appropriate. ometrical properties of the motion, once the vibrational
In this description, in the glassy state, the system is lo- component is subtracted. The possibility of performing
calized in one of the potential energy basins [7–10]. suchastudybelowT×,withaveryfinetimecoarsegrain-
1
ing, allows us to examine the structural changes that characterized by an energy change of ≈ 10−20 kJ/mol
accompany the basin transitions and to describe an ele- and an oxygen atom square displacement of the order of
mentarystepofthediffusiveprocessintermsofhydrogen 0.01 ˚A2; they appear to constitute the elementary step
bond network rearrangement. underlying the diffusional process in the system.
This work is organized as follows. In Sec. II we pro- We note that it is impossible to identify the transi-
vide the simulation details. In Sec. III we analyze the tion unless the quenching time step is of the order of
displacementofthemoleculesinthese IS-transitionsand the simulation time step. This is in distinct contrast
in Sec. IV we study the corresponding changes in the to a Lennard-Jones liquid, where such small continuous
hydrogen-bond (HB) network. Finally, in Sec. V we changesarenotpresent[7]. Thedifference,wewillshow,
present a brief summary. isattributabletothehydrogenbonds. Thistimescalefor
IS transitions is in accordwith other workon the TIPS2
water model, at T =298 K [19].
II. SIMULATION In Fig. 2, we see the sharp changes in EIS(t) coincide
withthesharpchangesinhr2(t)i. Thisconfirmsthatthe
Our results are based on molecular dynamics simula- systemisrepeatedlyvisitingspecificconfigurations,since
tionsoftheSPC/Emodel[20]ofwaterfor216molecules, hr2(t)i and EIS(t) take on discrete values. The results
at fixed density ρ=1 g/cm3 . The numerical procedure shown in Fig. 2 imply that the system often returns to
is described in Ref.[22]. The integration time step δt theoriginalbasinbecauseofboththedifferenceinenergy
is 2 fs. The mode coupling temperature for this den- andthedisplacementapproachzeroattheendofthetime
sity is TMCT ≈ 194K [22]. We analyze trajectories at interval.
T =180K,sothatthesystemisinthedeepsupercooled To aid in understanding the distribution of the dis-
liquid state. At this temperature, the diffusion coeffi- placementsduringtheISchanges(suchasthosebetween
cientis morethanfourordersofmagnitudesmallerthan the two IS labeled A and B in Fig. 2(b)), Fig. 3 shows
its value at T = 300 K and only a few molecules move the displacements u of all 216 individual molecules from
significantly (with displacements larger than 0.025 nm) IS-AtoIS-B.Weseethatthereisarelativelysmallsetof
at each simulation time step[24]. Our system is started molecules with a large displacement. A snapshot of the
from equilibrated configurations at 190 K, which relax eightmoleculeswiththelargestdisplacementisshownin
for nearly 920 ns at 180 K before we record and analyze Fig. 4. Interestingly, we find that this set of molecules
the trajectory. At such a low temperature, a slow aging formsaclusterofbondedmolecules. Indeed,forallcases
in the trajectory could be present, however, the aging studied, we found that the set of molecules which dis-
should not affect the qualitative picture we present. place most form a cluster of bonded molecules. The ob-
We have generated one trajectory of 30 ns, sampling served clustering phenomenon characterizes the IS tran-
configurations at each 1 ps. For each configuration we sitions in water and can be interpreted as the analog of
find the corresponding IS using conjugate gradient min- the string-like motion observed in simple atomistic liq-
imization. In this way, we obtain 30,000 configurations uids [7] and connected to the presence of dynamical het-
with the corresponding IS. Since we could miss some IS erogeneities[25]. Similar results were found by Ohmine
transitionswith1pssampling,wealsoranfourindepen- et al. using the TIP4P andTIPS2 models for water [26].
dent 20 ps simulations sampling the IS at 4 fs. In this To characterize the distribution of individual molecu-
way, we obtain another 20,000 configurations with the lardisplacementsbetweendifferentISmorecarefully,we
corresponding IS. show the distribution of displacements u of the oxygen
atoms P(u,δt = 4 fs) in Fig. 5. The time difference δt
is two times the integration time step[27]. We collected
III. INHERENT STRUCTURE TRAJECTORIES datafromalltransitionsinfourindependent10pstrajec-
tories for the results reported in Fig. 5. Note that P(u)
was previously studied by Schrøder et al. for a binary
Figure 1 shows an example of the potential energy
Lennard-Jones (LJ) mixture [7].
of the IS, EIS(t), and the mean square displacement of
the oxygenatomshr2(t)istarting froma singlearbitrary Wefindthatthedistributionofdisplacementsismono-
tonic and does not allow us to detect a characteristic
starting time. At T = 180 K, the slowest collective re-
length which could help in distinguishing between diffu-
laxation time τα > 200 ns [22]. The IS trajectory in
sive and non-diffusive basin changes. P(u,t) for t = δt
Fig. 1 has a mesh of 1 ps and covers a total time of
30 ns. In this time interval, hr2(t)i is about 1˚A2, i.e. shows an apparent power-law with a negative slope of
about −2.8± 0.2, followed by an exponential tail. To
much less than the corresponding value of the average
nearestneighbordistanceof2.8˚A. Figure2showsanen- highlight the exponential tail, we present the data in a
linear-logplotFig.5(upperpanel),showingthatthetail
largementofthe IStrajectoryusingamuchsmallertime
in P(u,δt = 4 fs) is mostly due to these highly “mobile
mesh (4 fs, two times the simulation time step). Fig. 2
molecules”. The characteristic length of the exponen-
shows that changes occur via discrete transitions, with
tial tail is 0.2˚A, about 15 times smaller than the nearest
an average duration of ≈ 0.2 ps. The transitions are
neighbor oxygendistance of 2.8˚Aand in agreementwith
2
the finding of Schrøder et al. for the LJ case, where the whose fraction is both temperature and density depen-
characteristiclength was about1/8of the nearestneigh- dent. The HB network tends toward a perfect random
bor distance. tetrahedral network on cooling, or on lowering the pres-
Schrøder et al. interpreted the power-law contribu- sure[23]. Here we emphasize the changes in the network
tionto P(u,δt)asaresponseofthe systemto the diffus- associated with IS-transitions.
ing molecules. In elasticity theory[28], a local displace- To quantify the bonding changes, it is necessary to
ment at the origin produces a distribution of displace- employ an arbitrary bond definition. Previous works
ments which scales with the distance R from the origin haveindicatedthat severaldefinitions providephysically
as R−2. Hence the distributionofdisplacements uscales reasonable results. Here we use the definition that two
as u−2.5[28]. In contrast with the LJ case, the exponent molecules are bonded if their oxygen-oxygen distance is
we find is slightly larger than 2.5, a discrepancy which lessthan3.5˚Aandtheirmutualpotentialenergyisnega-
may be related to the application of a continuum theory tive[31]. Usingthisdefinition,weobtaintheequilibrium
at a length scale where molecular details are still rele- distribution of the potential energy for the LB and BB,
vant. In the spirit of elasticity theory, the deviation of shown in Fig. 6.
the distribution at small u values from the power law These results agree with previous findings, based on
arises from the cutoff introduced by the finite size of the different models[32–34]. The distribution of BB ener-
simulation. Indeed, small u values are produced at large gies is bimodal with peaks at roughly −6 kJ/mol or
R and hence are missing in a finite simulation box. −22.5 kJ/mol, while the LB energy distribution is uni-
Althoughasubsetof“highlymobile”moleculesisiden- modal with a peak at roughly −24 kJ/mol. Therefore,
tifiable using a pre-defined threshold value in a single the energy associated with a change in the HB network
basin change, there is no unambiguous general criterion duetolosingoneLBandcreatingtwoBBrangesroughly
foridentifyingthemoleculesresponsibleforasinglebasin from −21 kJ/mol to 12 kJ/mol, depending on which of
transition. the two possible BB are created. The relative intensities
An important open question is the relation between of the peaks of the BB energy distribution suggest that
thedisplacementdistributionfunctionsfort=δt(Fig.5) such a mechanism would more likely lead to an increase
andthesamedistributionevaluatedforanarbitrarytime in the overall energy. The comparison of the HB-energy
interval t. A full knowledge of such relation may shed changes with those found in Fig. 1 suggest that inter-
light on the elementary stochastic process which better changesbetweenLBandBBcanexplaintheexistenceof
describes the dynamics of the deep supercooled state. amultitudeofISseparatedbyverysmallenergybarriers.
This calculation requires a significant computational in- More specifically, we hypothetize that changes increas-
vestment due to the very large number of IS configura- ing EIS found in Fig. 1 are due to processes LB→BB,
tions which need to be evaluated. Despite such difficul- while the changes decreasing EIS are due to the pro-
ties, it is a promising researchline for the near future. cesses BB→LB. To confirm or reject this hypothesis, we
study the distributions of LB and BB for IS just before
and just after positive and negative “jumps” in the po-
IV. HB NETWORK RESTRUCTURING tentialenergy(Fig.1b)withenergylargerthan9kJ/mol
[35]. These IS are schematically shownin Fig. 7(a). The
Liquid water is an interesting system for studying the correspondingdistributions for LB and BB are shown in
physicalprocessesthataccompanybasintransitions. The Fig.7(b) and 7(c). We see that during PE increases,the
unusualdynamicandthermodynamicpropertiesofwater averagenumberofLBdecreaseswhile the number ofBB
are believed to be connected to the microscopic behav- increases; the opposite situation occurs when the poten-
iorofhydrogenbonding. Manyexperimentssuggestthat tial energy decreases. The distribution for BB has peaks
this tetrahedral HB network has defects, such as an ex- at even numbers of BB, which we expect since for each
tra (fifth) molecule in the first coordination shell[29,30]. LB lost, two BB appear, also implying that the distri-
Indeed, such 5-coordinated molecules have been directly bution should be zero for odd numbers of BB [36]. The
identifiedinsimulationsandthedefectswerefoundtobe anti-correlation of the BB and LB changes is a strong
a catalyst for motion in the system[31], making an ob- indicatorthata mechanismwhereby the systemaccesses
vious possible connection between network defects pro- higher energy states is via LB→BB transitions.
moting diffusion and the basin transitions that give rise To reinforce the hypothesis that the basin change is
to diffusive motion of that system. associated with a restructuring of the local connectivity,
To obtain a physical picture of the IS transitions for Fig. 8 shows the number of molecules with a coordina-
water and to hopefully better understand the source of tionnumber equalto 3,4 or 5 asa function oftime for a
these transitions, we focus on the small changes in the characteristic time interval and contrasts this data with
HB network along a minimum to minimum trajectory. the time dependence of hr2(t)i. A clear anti-correlation
Analysis of the HB network based on IS configurations is observed between the time dependence of the number
has shown that local PES minima are characterized by of 3 and 5-coordinated molecules compared to the time
both linear bonds (LB) and bifurcated bonds (BB) [31] dependence of the 4-coordinated molecules, supporting
the proposed interchange mechanism. The fact that in-
3
creases in the fraction of BB coincides with changes in [3] F.H. Stillinger, Science 267, 1935 (1995).
hr2(t)i supports the expectation that motion is a result [4] C.A. Angell, Science 267, 1924 (1995).
of network imperfections. [5] S. Sastry, P.G. Debenedetti, and F.H. Stillinger, Nature
393, 554 (1998).
[6] P.G. Debenedetti, and F.H. Stillinger, Nature 410, 259
V. SUMMARY (2001).
[7] T.B. Schrøder, S. Sastry, J.C. Dyre,and S.C. Glotzer, J.
Chem. Phys.112, 9834 (2000).
We have presented a detailed analysis of the motion
[8] A. Heuer,Phys.Rev.Lett. 78, 4051 (1997); S.Buchner,
on the PES of a 216 water molecule system interacting
and A.Heuer, Phys.Rev. E 60, 6507 (1999).
via the SPC/E potential in a deeply supercooled state,
[9] L.Angelani,G.Parisi,G.Ruocco,G.Viliani,Phys.Rev.
below the MCT temperature. At this temperature, the
Lett. 81, 4648 (1998).
system populates basins of the local minima. At these
[10] F.Sciortino, W.Kob,andP.Tartaglia, Phys.Rev.Lett.
conditions, the analysis of the IS trajectories provides a 83, 3214 (1999).
very clean description of the slow alpha relaxation pro-
[11] A.Scala,F.W.Starr,E.LaNave,F.Sciortino,andH.E.
cess filtering nearly all vibrational motions. Stanley, Nature406, 166 (2000).
We have shown that an inherent structure transition [12] J.P.K. Doye, and D. Wales, Phys. Rev. Lett. 86, 5719
is observed about every 0.2 ps. It is the collection of (2001).
these numerous small transitions that gives rise to the [13] T. Keyes,J. Phys.Chem. A 101, 2921 (1997).
structural relaxation of the system. These fast transi- [14] E.LaNave,A.Scala,F.W.Starr,F.Sciortino,andH.E.
tions are characterized by a broad distribution of indi- Stanley, Phys. Rev. Lett. 84, 4605 (2000); E. La Nave,
vidual molecule displacements, without a clear charac- A. Scala, F. W. Starr, H.E. Stanley, and F. Sciortino,
teristic length. Future work must address the issue of Phys. Rev.E 64, 36102 (2001).
the prediction of P(u,t) from the knowledge of P(u,δt). [15] F. Sciortino, and P. Tartaglia, Phys. Rev. Lett. 78 2385
We perform an analysis of the geometry of the indi- (1997).
vidualeventandfindthatthemostmobilemoleculesare [16] C. Donati, F. Sciortino, and P. Tartaglia, Phys. Rev.
clustered. The analysis of the changes in HB connectiv- Lett. 85, 1464 (2000).
ity associated with IS changes reveals that these transi- [17] E. La Nave, F. Sciortino, and H.E. Stanley, Phys. Rev.
Lett. 88, 035501 (2002).
tions are associated with the breaking and reformation
[18] W. G¨otze, and L. Sjøgren, Rep. Prog. Phys. 55, 241
of HB. This result is in accord with work of Ohmine
(1992).
and colleagues [19,26] for TIP4P. We have shown that
[19] I. Ohmine, and H. Tanaka, Chem. Rev. 93, 2545 (1993)
the transitions associatedwith an increase in the energy
and references therein.
correspond to the breaking of linear bonds and to the
[20] H.J. Berendsen, J.R. Grigera, and T.P. Stroatsma, J.
simultaneous formation of bifurcated bonds. Similarly,
Phys. Chem. 91, 6269 (1987).
the transitions associated with a decrease in the energy
[21] F. Sciortino, P. Gallo, P. Tartaglia, and S-H. Chen,
correspond to the breaking of bifurcated bonds and to Phys. Rev. E 54, 6331 (1996); P. Gallo, F. Sciortino,
the simultaneous formation of linear bonds. P. Tartaglia, and S.H. Chen Phys. Rev. Letts 76 2730
This result supports the hypothesis that the linear to
(1996); F. Sciortino, L. Fabbian, S-H. Chen, and P.
bifurcatedtransitioncanbe consideredas anelementary Tartaglia, Phys. Rev. E 56, 5397 (1997); L. Fabbian,
step in the rearrangementof the HB network. A. Latz, R. Schilling, F. Sciortino, P. Tartaglia, and C.
Theis Phys. Rev. E. 60 5768 (1999); ibid. Phys. Rev.
E. 62, 2388 (2000); C. Theis, A. Latz, R. Schilling, F.
VI. ACKNOWLEDGMENTS Sciortino, P. Tartaglia Phys.Rev.E. 62, 1856 (2000).
[22] F.W.Starr,F.Sciortino,H.E.Stanley,Phys.Rev.E60,
WethankT.B.Schrøder,S.SastryandS.C.Glotzerfor 6757 (1999).
interesting discussions. N.G. wants to thanks I. Ohmine [23] F.W.Starr,S.Sastry,E.LaNave,A.Scala,H.E.Stan-
for fruitful discussions and suggestions. F.S. acknowl- ley, and F. Sciortino, Phys. Rev.E 63, 041201 (2001).
[24] The Adam and Gibbs hypothesis [37] states that flow
edgessupportfromINFM-PRA-HOPandInitiazivaCal-
in a supercooled fluid involves cooperative motion of
colo Parallelo and MURST COFIN2000. This work was
molecules. For recent experimental work on cooperative
supported by the NSF Chemistry Program.
motions see Ref.[38].
[25] M. Hurley, and P. Harrowell, Phys. Rev. E 52, 1694
(1995); W.Kob, C. Donati, S.J. Plimpton, P.H. Poole,
and S.C. Glotzer, Phys. Rev. Lett. 79, 2827 (1997); C.
Donati,J.F.Douglas,W.Kob,S.J.Plimpton,P.H.Poole,
and S.C. Glotzer, Phys. Rev. Lett. 80, 2338 (1998); B.
Doliwa, andA.Heuer,Phys.Rev.Lett.80,4915(1998).
[1] M. Goldstein, J. Chem. Phys. 51, 3728 (1969). [26] M.MatsumotoandI.Ohmine,J.Chem.Phys.104,2705
[2] F.H. Stillinger, and T.A. Weber, Phys. Rev. A 28, 2408 (1996); I.OhmineandS.Saito,Acc.Chem.Res.32,741
(1983).
4
(1999), and references therein.
[27] Ifduringthetimeintervalδtthesystemdoesnotchange 0.015
IS, the contribution to P(u,δt) appears only for u = 0 (a)
andhencedoesnotmodifytheshapeofthedistribution. 2]
[28] J.C. Dyre, Phys. Rev. E 59, 2458 (1999); ibid. 59, 7243 m 0.010
n
(1999). > [
[29] J. Dore and J. Texeira (eds.), Hydrogen Bonded Liquids (t) 0.005
(Kluwer, Dordrecht),(1991), pp 171-183. 2
r
[30] E.Gru¨nwald,J.Am.Chem.Soc.108,5719(1986); A.H. <
Narten and H.A. Levy, Science 165, 447 (1969); P.A. 0.000
Giguere,J.Chem.Phys87,4835(1987); G.E.Walrafen,
M.S. Hokmabadi, W.H. Yang,Y.C. Chu and B. Mono- ol] −129.8 (b)
smith, J. Phys. Chem. 93, 2909 (1989). m −130.0
[31] F. Sciortino, A. Geiger, and H.E. Stanley, Phys. Rev. kJ/
Lett.56,3452(1990); Nature354,218(1991);J.Chem.20 −130.2
1
[32] FP.hHy.s.S9ti6ll,in3g8e5r7a(n1d99A2.).Rahman,J.Chem.Phys.60,1545 (t) [ −130.4
S
(1974). EI −130.6
[33] H.J.C.Berendsen,J.C.M.Postma,W.F.VonGunsteren, 0 10 20 30
and J. Hermans in Intermolecular Forces edited by B.
t [ns]
Pullman (Reidel, Dordrecht,1981). p.331.
[34] W.L. Jorgensen, J.Chandrasekhar, J.D. Madura, R.W.
Impey,and M.L. Klein, J. Chem. Phys. 79, 926, (1983).
FIG. 1. (a) Mean square displacement and (b) IS energy
[35] We find from Fig. 1 that most of the “jumps” in poten-
fortheinherentstructuresas afunction oftimefor thestud-
tial energy correspond to changes in |EIS| larger than ied216-molecule system. Thetimeintervalbetweenadjacent
9kJ/mol. We found similar results for a lower cutoff in
IS in both figures is 1 ps. While it is possible to track IS
energy of 1kJ/mol.
transitions from the potential energy, it is not the case for
[36] Wefind that thedistribution is not exactly zero for odd
the mean square displacement. Note the amplitude of the
numbersofBB,whichisaconsequenceofthebonddefini-
peaks of the potential energy is ≈10−20 kJ/mol, the same
tion.Usingthedefinitionthattwomolecules arebonded
order of magnitude of a HB energy.
if their oxygen-oxygen distance is less than 0.35 nm
and their interaction energy is less than a fixed value
EHB = 0, we find that there exist few cases in which
one hydrogen atom is “shared” by 3 oxygens (instead of
2 as in the case of a BB). We try different values for
EHB and we find that these “tri-furcated” bonds (TB) −130.1
disappearifEHB isreducedto-10kJ/mol.TheseTBcan
alsocontributetothepresenceofanevennumberof BBol] (a)
m −130.2
appearing in the distribution. The ratio of BB:TB for
J/
tEhHaBt t=he0d,is−tr5ibkuJt/iomnofloirsB2B0:1atanEdH2B4:=1. −W1e0aklsJo/mchoelciks20 k −130.3
zero for odd numbersof molecules. 1
[[3378]] GM..TA.daCmicearnodneJ,.HF..GRi.bbBsl,aJc.kCbuhrenm.aPnhdysM.4.D3.,1E3d9ig(1e9r,65J)..E(t) [IS −130.4
Chem.Phys.102 471 (1995); E.R. Weeks,J.C. Crocker, −130.5
A.C. Levitt, A. Schofield, D.A. Weitz, Science 287 627
(2000); K. Schmidt-Rohr and H.W. Spiess, Phys. Rev. B
Lett.663020(1991); A.Heuer,M.Wihelm,H.Zimmer- 0.0002 (b)
mannandH.W.Spiess,Phys.Rev.Lett.752851(1995).2m]
n
> [
) 0.0001
(t
2
r
<
A
0.0000
0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4
t [ps]
5
FIG.2. (a)ISenergyand(b)meansquaredisplacementfor
the IS obtained using a sampling interval of 4 fs, two times
the simulation time step. The correlation between EIS and
hr2(t)i is evident. Also, we see that it is necessary to sample
IS with a mesh of the order of the simulation time step to
detect all theIS visited by thesystem.
0.08
0.06
]
m
n 0.04
[

u
0.02
FIG. 4. Snapshot of the system in the IS labeled A in
0 Fig. 2(b). Only the eight molecules with displacement larger
0 30 60 90 120 150 180 210 than 0.025 nm are shown here. Hydrogen-bonded molecules
areconnectedbytubes. Notethatall8molecules arenearby
molecule number
andform acluster,whichunliketheLJcase,arebondedand
less string-like.
FIG. 3. Displacement of each molecule in the transition
from IS-A to IS-Bshown in Fig. 2(b).
103
(a)
102
~ exp(−u/0.021)
101 ~ u−2.8
) ~ u−2.5
u, t 100
(
P
10−1
10−2
10−3
0.00 0.02 0.04 0.06 0.08 0.10
u [nm]
6
FIG. 6. Distribution function for the pair interaction en-
103 ergy Vij for the LB and BB in the IS. Note the bimodal dis-
tribution for BB and unimodal distribution for LB.
(b)
102
101
)
u, t 100
(
P
10−1
~ exp(−u/0.021)
10−2 ~ u−2.8
~ u−2.5
10−3
0.001 0.010 0.100
u [nm]
FIG. 5. Distribution of displacements u of the oxygen
atoms between ISchanges, P(u,t),evaluated at t=4 fs. We
showthedistributionofallthedisplacementsofthemolecules
(squares) and the distribution for the largest displacement
in an IS-transition (circles). For comparison, we divided the
largestdisplacementdistributionbythenumberofmolecules.
Wepresentboth(a)log-logand(b)linear-loggraphstoshow
the power law behavior for u ≈ 0.01 nm and the exponen-
tial tail of the distributions. We also show the prediction of
elasticity theory P(u)∼u−2.5.
0.15
LB
0.10 BB
)
Vij
(
P
0.05
0.00
0 10 20 30
− V [KJ/mol]
ij
7
FIG.7. Panel(a) showsschematically theseparation of IS
usedtomeassuretheprobabilitydistributionsfor(b)LBand
(c)BB.Circlescorrespondtothedistributioncalculatedover
IS just before an increase in potential energy larger than 9
kJ/mol, i.e. when the system is in an IS A shown in (a).
(a) SquarescorrespondtothedistributioncalculatedoverISjust
afteranincreaseinpotentialenergylargerthan9kJ/mol,i.e.
when the system is in an IS B shown in (a). Similarly, dia-
y
g monds correspond to the distribution calculated over IS just
r
e before a decrease in potential energy larger than 9 kJ/mol,
n
E B C when the system is in an IS C shown in (a). Triangles cor-
S respond to the distribution calculated over IS just after a
I decrease in potential energy larger than 9 kJ/mol, when the
>9kJ/mol system is in an IS D shown in (a).
A D
Time 30
CN=4
(b) 20
n 0.2 N CN=5
o C
ti
u 10
b A
ri
t
s B
Di CN=3
C 0
y
t
uli 0.1 D 2m] 0.0002
b
a n
b [
o > 0.0001
r )
P (t
2
r
< 0
0 0.2 0.4 0.6 0.8 1 1.2 1.4
0.0
t [ps]
415 425 435
number of LB
FIG. 8. Number of molecules with coordination number
0.2
CN equals 3,4 and 5 versustime, for the IScorresponding to
Fig. 2. The plot for CN = 4 is shifted down by 400 units
(c) forbettercomparison. Alsoshownishr2(t)i. Weseehowthe
n tetrahedralnetworkacquiresbothtypesofdefects(CNequals
o
ti to 3 and 5) while the system explores different IS.
u
b A
ri
t B
s
Di 0.1 C

y
t D
bili
a
b
o
r
P
0.0
0 10 20 30
number of BB
8