Classification of local structures: Voronoi groups

D. Classification of local structures: Voronoi groups

The Voronoi polyhedral (VP) analysis has been applied for the computer simulations of liquid water [7][8][9]. Considering a simple liquid system, the Voronoi cell of particle i is defined as the cell where any point inside it is closer to that particle than to any other particles. The geometry of the Voronoi polyhedron is, therefore, a generalization of the Wigner-Seitz unit cell in a crystal. However, in liquid water, there are three atoms in a water molecule so that it is hard to define Voronoi unit cell. Nevertheless, we can apply this idea to either the centers of mass of molecules or the oxygen atoms of every molecules. Here, we choose to apply for the oxygen atoms. The algorithm for determining the vertexes of a Voronoi unit cell is given in the program F. 35 of Ref. [4].

For a Voronoi unit cell, the asphericity

ξ

is defined as

2 3

36 V

A ξ

=

π

,

where A and V are the total surface area and the volume of the Voronoi cell [10]. The

smallest

ξ

equals one for a sphere. As

ξ

increases, the geometry of the polyhedron becomes more aspherical. For the ice

I structure, in which most molecules are in

^h the tetrahedral arrangement,

ξ

is equal to 2.25.

We have performed the VP analysis for the oxygen atoms in our simulated liquid water and the asphericity distribution of the VP has a shape similar as the one of the TIP4P water model at the same temperature [11]. By dividing the whole distribution into four adjacent sections, with their ranges given in Table II, we separated the molecules of a configuration into four groups, indexed as Voronoi group (VG) I to IV with increasing asphericity. Most of the molecules are in VG II and III; VG I and IV contain only few percentages. The selection operator

Θ

_ξ

(i )

for molecule i in VG

ξ

is defined as

To examine the local structures of the molecules in each VG, we present in Fig. 3 the oxygen-oxygen, oxygen-hydrogen and hydrogen-hydrogen radial distribution functions for each of the four VGs, which are defined as

∑

^{Θ >}

∑

^{≠ − Θ} local structures of liquid water can be distinguished into two extremes: one is the local structure of the molecules in the VG I, which has a randomly structure indicated by

⎪⎩

the plateau from r = 2.76 Ǻ to 3.81 Ǻ in

g

_OO^I (r) and the other is the similar ordered

structure represented by the VG II, III and IV. The plateau in

g

_OO^I (r) shows that the local structures of water molecules in the VG I have much more random in the radial direction and more weakly connected through H-bonds with their neighbors than the other three VG groups. Thus, from the VP analysis, a structural description for the two-state mixture model of liquid water has been proposed [11].

We also calculated the length distributions of the H bonds attaching to the molecules in each VG and the results, also shown in Fig. 3, obviously point out that the distribution shifts toward shorter length as the asphericity increases.

For the dynamical properties of the water molecules in various VGs, we consider the AVAF of the VG

ξ

defined in the following

By defining the projection operator for the VG

ξ

_as

∑

=

the total INM DOS can be separated into components as

( )

( ^{( )} )

or obtained by the INM approximation given in eq. (III.3). These results will be shown in next chapter.

Chapter IV

Results and Discussions

In this chapter, we will focus on the INM DOS, the angular velocity

autocorrelation function and the behaviors of the orientational correlation function at the short-time, long time and intermediate-time scales.

A. The instantaneous normal mode DOS

Before discussing the INMs, we schematically describe in Fig. 1 the body frame of a SPC/E water molecule. In the body frame, a water molecule is placed on the x-y plane and the center of mass of the molecule locates at the origin. The oxygen atom is on the positive x-axis and the line jointing the two hydrogen atoms are parallel to the y-axis. For such a molecule, the moments of inertia along the x-, y-, and z-axis are IX=7.41×10⁻², IY=3.29×10⁻², and IZ=1.07×10^{−1 m}Ų, respectively. Here, ^m is the mass of a water molecule.

In Fig. 4, we present the INM DOS and its translational and the rotational components. A distinct difference in the imaginary spectrum of the DOS between liquid water and the Lennard-Jones liquids is found. In liquid water, the fraction of the imaginary spectrum is only about 12 percents of the INM DOS, while in Lennard-Jones liquids, more than 25 percents (depending on the density of the liquid) of the DOS are contributed from the imaginary spectrum [12]. This suggests that the

water system spends most of time around the local minima of the potential energy surface and the water molecules are hard to diffuse. We find that the maximum peak in the real-frequency spectrum occurs at about 70

cm . This peak can be explained by

⁻¹ the oxygen-oxygen bending motion [13] and the peak position is consistent with the prediction of the random network model [14]. Besides the translational spectrum, the peaks of the rotational INM spectrum of the x-, y- and z-axis occur at about 500, 830, and 530

cm , respectively. The peak positions of the rotational INM DOS along the

⁻¹ three axes are related to the moments of inertia associated with the corresponding axis.

For the y-axis, the moment of inertia is the smallest, which causes the largest in the mean value of the rotational DOS. It is expected that the mean value of the rotational DOS along the x-axis is larger than that along the z-axis, since the moment of inertia along the x-axis is smaller than that along the z-axis. However, the mean values of the rotational DOS along the x- and z-axis are almost the same. This is understood for the involved intermolecular interactions with surrounding molecules for rotating along the z-axis is stronger than that along the x-axis [15].

By the VP analysis in eq. (III.4), we can decompose the rotational INM DOS to investigate the effect of local structures. In Fig. 5, we normalize all rotational INM spectra along the three axes for various VGs. There are two interesting features. For all rotational axes, with the increase of the asphericity, the mean frequency of the real rotational INM DOS shifts toward the high-frequency end. Also, the fraction of the imaginary rotational DOS of a VG is found to decrease with the increase of the asphericity. The blue shift in the rotational DOSs for different VGs can be explained by the fact that the number of the hydrogen bonds (HBs) increases with the asphericity. In the VG IV, where most water molecules are bonded by four HBs, the mean frequency in rotational DOS is higher than that in the VG I, where a water molecule has fewer HBs [5].

B. Angular velocity autocorrelation function

Shown in Fig. 6 are the AVAFs along the x-, y-, and z-axis calculated by MD simulations. At the short-time scale, the AVAF along the y-axis decays faster than those along the x- and z-axis. Less attenuation is found in the AVAF along the z-axis and the AVAFs along the x- and z-axis are quite similar.

We also use the VP analysis in eq. (III.2) to discuss the AVAFs. In Figs. 8(a), 8(b), and 8(c), the simulated results of the AVAFs along x-, y-, and z-axis for the four VGs are shown. It is found that at the short-time scale, the AVAF becomes faster decaying with increasing the asphericity. This is consistent with the blue shift in the INM rotational DOSs for different VGs.

C. Short-time behavior of

C_μ⁽¹⁾

(

t

) and

^Cμ(2)

⁽

^t

⁾

In this section, we focus on the short-time behaviors of the first and the second order orientational correlation functions (^Cμ^(l)

⁽

^t

⁾

, l = 1 or 2,

μ

= x, y, or z) for t smaller than 0.2ps. For these OCFs under the second cumulant approximation given in eq. (II.8),

ψ

(t) can be calculated by three different ways. Two of these are calculated by the total INM DOS and the real-frequency INM DOS; the other is by the power spectrum of the AVAF.

Figs. 9(a), 9(b), and 9(c) show the ^D_μeff

( ω )

^kb^T , with μ = x, y, z, of the AVAF power spectrum and those of the stable INMs for the four VGs. The shape and the area for the ^{D k}b^T

eff

( ω )

μ spectrum of the z-axis are found to be similar to those for

the spectrum of the x-axis and are larger than those for the spectrum of the y-axis. In Fig. 4, Fig. 5, the DOS of the x-axis is similar to that of the z-axis. D_μ^eff

⁽ ω ⁾

is a combination of the DOSs of the others two axes divided by the corresponding moment of inertia. This is the reason why the shapes of

D

^eff_X (

ω

)

k

_b

T

_and

T k

D

_Z^eff(

ω

) _b are similar and the area of

D

_Y^eff(

ω

)

k

_b

T

is the smallest. In addition to

the spectrum shape and area, ^{D k}b^T eff

( ω )

μ is 0 at ω = 0, but the power spectrum of AVAF is not. Also, as the asphericity of the VG increases, the value of ^{D k}b^T

eff

( ω )

μ

at ω = 0 is decreased. This difference in the spectrum values at ω = 0 makes the OCF

)

)

(

( t

C_μ^l calculated by the real-frequency INM and the power spectrum of AVAF different.

In Figs. 10 and 11, the three fitting lines in the short-time scale (within about 20 fs) can exactly fit the MD simulation results. Among these fittings, the best one for the longest time is the ^Cμ^(l)

⁽

^t

⁾

calculated by the power spectrum of AVAF, because the AVAF is calculated by the MD simulation and

D

_η(

ω

) in eq. (II.12) is obtained from the power spectrum of Ω_η(t) by the inverse Fourier cosine transformation. For this reason, the ^Cμ^(l)

⁽

^t

⁾

calculated by the power spectrum of AVAF has the longest

accurate fitting time than other. The ^Cμ^(l)

⁽

^t

⁾

calculated by the total INM DOS provides an accurate estimation within t < 20 fs, but fails at longer times. Eq. (II.15) point out that the total INM DOS is larger than that of the real-frequency INMs such that D_μ^eff

⁽ ω ⁾

in eq. (II.8) will decay faster, and the ^Cμ^(l)

⁽

^t

⁾

calculated by the real-frequency INMs decays slowly than the results of MD simulation. On the other hand, comparing the results of the x-, y-, and z-axis, we can see the four lines of

y-axis decay slowly than others, and the z-axis is fastest. Because the area of

T

k

D

_Y^eff(

ω

) _b is small than others, in eq. (II.8), (II.8-2), and (II.12), tell us when the

area of D_μ^eff

⁽ ω ⁾

increase, the decay speed of ^Cμ^(l)

⁽

^t

⁾

will increase too.

As t < 0.2ps, both the ^Cμ^(l)

⁽

^t

⁾

of the real-frequency INMs and the power spectrum of AVAF have good predictions, but after t = 0.2ps, they have different behaviors. We have discussed previously that the ^{D k}b^T

eff

( ω )

μ at ω = 0 for the real-frequency INMs and the power spectrum of AVAF are distinct. Because, when t > 0, ^Dμ^eff

^{( 0 )}

will be more important influence on ^Cμ^(l)

⁽

^t

⁾

, so as time increase, the ^Cμ^(l)

⁽

^t

⁾

of the power spectrum of AVAF will decay faster than the real-frequency INMs. Besides, comparing ^Cμ⁽¹⁾

⁽

^t

⁾

and ^Cμ⁽²⁾

⁽

^t

⁾

, the real-frequency INMs have a better qualitatively correct behavior than the power spectrum of AVAF in ^Cμ⁽²⁾

⁽

^t

⁾

, but when t < 35fs, the

power spectrum of AVAF have longer accurate fitting time either ^Cμ⁽¹⁾

⁽

^t

⁾

or ^Cμ⁽²⁾

⁽

^t

⁾

. By VP analysis, in Figs. 12(a), (b), (c) and Figs. 13(a), (b), (c), we can observe that as VG

ξ

increases, the number of bumps increases. This is very reasonable, because we know that as

ξ

increases, the average value of H-bond length decreases. The H-bond will restrict water molecular motions, like inside a cage. When the H-bond length decreases, the restricting force will be stronger and make the rebound in molecular motion occur. So, the bumps come from the reversion or palliation of rotational motion in water molecule. That is the reason why the bumps increase with

ξ

increase. From VP analysis, there is another appearance. When

ξ

increasing, the accurate fitting time will be extended in the ^Cμ^(l)

⁽

^t

⁾

of the real-frequency INMs and the power spectrum of AVAF. This may reveal that the numbers of H-bond affect the

accuracy of prediction by the real-frequency INMs.

From the above discussions, we find that the short-time behavior of ^Cμξ^(l)

⁽

^t

⁾

is a Gaussian decay within t < 20fs and, even for different VGs, can be accurately fit by the three kinds of method. When t > 0.2ps, these methods do not provide correct predictions. Therefore, we have to use another theory to describe their behaviors at the long-time scales.

D. Long-time behavior of

C_μ⁽¹⁾

(

t

) and

C_μ⁽²⁾

(

t

)

In the previous section, we use the total, real-frequency INMs and the power spectrum of AVAF to fit the short-time behavior of the OCF. But, for the long-time scales, these methods do not provide correct predictions. By eq. (II.18), the long-time behavior of ^Cμ^(l)

⁽

^t

⁾

is expected to be a single-exponential decay. Thus, we change to fitting parameters for different rotational axes and VGs in Table III and Table IV. The linear fitting function f_μ^(l)

(

t

)

for the

μ

-th principal axis and for VG

ξ

are written,

where a_μ^(l) is a constant,

τ

^D(lμ⁾ is the rotational-diffusion decay time, l is the index for the l-th order Legendre polynomial, and

D is the rotational diffusion coefficient for

_μ the

μ

-th principal axis.

In Fig. 14 and Fig. 15,

^ln

^Cμ^(l)

⁽

^t

⁾

is well fit by the linear function ^fμ^(l)

⁽

^t

⁾

so that the long-time behavior of

^ln

^Cμ^(l)

⁽

^t

⁾

is linear in time for each axis. For ^fμ⁽²⁾

⁽

^t

⁾

, although the

τ

^D(lμ⁾ values of the three axis are different but the three lines are very

close to parallel. For ^fμ⁽¹⁾

⁽

^t

⁾

, the y-axis decays faster than the x-axis instead of parallel to each other.

By VP analysis, in Figs. 16(a), (b), (c) and Figs. 17(a), (b), (c), the

^ln

^Cμξ⁽¹⁾

⁽

^t

⁾

and

)

(

ln

C_μξ⁽²⁾ t , with

μ

= x-, y-, or z-axis, also can be well fit by the linear function

)

)

(

( t

f_μξ^l for each VG

ξ

. In Table III, and IV, we can find that the changes of

ξ

_and

) (l D_μξ

τ

are not related and, thus, these lines are almost parallel. According to these

results, the long-time behavior of ^Cμξ^(l)

⁽

^t

⁾

calculated by MD simulation is an exponential decay. At the short-time scales, we know the bumps of water molecule motion increase when VG

ξ

increases. These bumps will decrease the decay speed.

Therefore, in the long-time scale, the sequence of the ^Cμξ^(l)

⁽

^t

⁾

from MD simulation in different VG

ξ

is reasonable. Moreover, these nearly parallel lines show that the decay rates in long-time scale are independent of their VGs.

In this section, we have already discussed the long-time behavior for different axes and the four VGs. We know that the long-time behavior is an exponential decay when t > 2ps. But, we find that the intermediate time scale behavior (0.2ps. < t < 2ps) is

fitting function to fit the intermediate time-scale behavior of ^Cμ^(l)

⁽

^t

⁾

and the discussions are given in the next section.

E. Intermediate time scale of

^Cμ⁽¹⁾

⁽

^t

⁾ and

^Cμ⁽²⁾

⁽

^t

⁾

Although the intermediate time behavior is not described by a single Gaussian decay or a single exponential decay, we think that the transition process from the short-time scale to the long-time scale must have parts of the residual memory.

Therefore, we use a fitting function that is composed of two kinds in mathematical formula. In Ref. [16], they use a combination of two terms of power law and multiplied by exponential function to fit the intermediate time-scale behavior of the orientational correlation function of water in the regime of supercooled liquids. Thus we assume the transition process is a power-law-like behavior.

Before constructing the fitting function in the intermediate time scale, we make

)

To analyze the intermediate-time behavior (within 0.2ps < t < 2ps), we first extend the fitting function ^Fμ^(l)

⁽

^t

⁾

for the long time behavior from the range t > 2ps to t > 0.2ps. The results for their comparison with the simulated ^Cμ^(l)

⁽

^t

⁾

are shown in the left-hand panels in Fig. 18, and Fig. 19. We find that the single exponential decay

for the long time behavior is not fit well within 0.2ps < t < 2ps. So, we define

which is a power-law function multiplied by a Gaussian function. In the fitting function,

s and

^(l)

p are parameters of the power-law function and

^(l)

τ

_i^(l) is the rotational-diffusion decay time in the intermediate-time scale.

For each VG, the corresponding fitting function is given as

].

axes. In the previous section, we know that the bumps of water molecule motion will decrease the decay rate. For different VGs, we can clearly observe that the relation between

ξ

and the decay rate in Figs. 20(a), (b), (c) and Figs. 21(a), (b), (c).

From these results, we also tabulate these fitting parameters for different axes and VGs in Table V and Table VI, respectively. We find that the appearance time of bumps will be extended when

ξ

increase. This result is consistent with the previous one, especially for the y-axis, the geometry make IY smaller than other two axes, therefore, the initial motion memory keeps in a longer period. Besides, ^Cμ⁽¹⁾

⁽

^t

⁾

decays more slowly than ^Cμ⁽²⁾

⁽

^t

⁾

. Obviously, in Fig. 20(b),

Δ

^Cμξ⁽¹⁾

⁽

^t

⁾

in VG III and

IV have inflection points.

The fitting functions provide good predictions to us in the intermediate-time behavior. This result indicates that the OCF after getting rid of the H-bond restrictions in the intermediate-time scale is a combination of the residual memory of the Gaussian decay and the exponential decay extending into a long-time scale.

Chapter V Conclusions

By MD simulations with the SPC/E model, we have investigated the behaviors of the OCF of water molecules at different time scales. For examining the effects of interaction force, especially the H-bonds between molecules, we use the VP analysis to distinguish the local structures of molecules in a configuration.

In chapter II, we have derived the short-time behavior of the OCF, which is a Gaussian decay. We show that the behavior of the OCF within t < 0.02ps can be well described with three kinds of approximation: the total and the real-frequency INM spectra and the power spectrum of the AVAF. For t > 0.02ps, some inflection points in the OCF occur, with the number of inflection points determined by the strengths of the H-bonds attaching to a molecule and, therefore, the H-bond strengths determined by the local structures of the molecule. Thus, the local structures of water molecules have an enormous impact in the short-time behavior of the OCF.

According to the Debye theory for the rotational diffusion in molecular liquids, the long-time behavior of the OCF is a single-exponential decay with the rotational-diffusion decay time given as ⁽⁾ =

[

( +1) _μ

]

⁻¹

τ

_D^lμ

l l D

, where l = 1, 2 and Dμ. is the rotational-diffusion coefficient. Thus, the ratio

τ

_D⁽¹_μ⁾/

τ

⁽_D²_μ⁾ of the decay times for the OCFs with l=1 and 2 in the long time limit is expected to be three. For my MD simulation results, the logarithm of the OCFs with l=1 and 2 in long-time limit are

indeed well. Fit with linear functions; however, given in Table III and Table IV, the ratio of the slopes of the linear functions for l=1 and 2 is about 1.6 ~ 2.2. This phenomenon has been discussed in paper [21], in which similar MD simulation results of the SPC/E model are reported. In this paper, the difference between the simulated ratio and the prediction by the Debye theory is explained by the rearrangements of the H-bond network by constantly breaking and forming H-bonds through the reorientation of water molecules. Therefore, the Debye diffusive model is unsuitable for describing the reorientational dynamics in water and the authors of the paper introduced the extended jump model (EJM) for describing the reorientations in liquid water. Although the Debye theory can not predict accurately the rotational-diffusion decay time, the long-time behaviors of the OCFs still follow a single-exponential decay.

For investigating the intermediate-time behavior of the OCF, I introduce

Δ

^Cμ^(l)

⁽

^t

⁾

, which is the difference between the OCF obtained by the MD simulation and the single-exponential decay in the long time limit. In general, the behaviors of

Δ

^Cμ^(l)

⁽

^t

⁾

for molecules with different characteristics of local structures are all well fit by a power-law function multiplied by a Gaussian function, only with different coefficients in the fitting function, which are resulted from the H-bond effect and the three different moments of inertia of the anistopic water molecule. For a special case, the H-bonds attaching to molecules with highly aspherical local structures cause an

for molecules with different characteristics of local structures are all well fit by a power-law function multiplied by a Gaussian function, only with different coefficients in the fitting function, which are resulted from the H-bond effect and the three different moments of inertia of the anistopic water molecule. For a special case, the H-bonds attaching to molecules with highly aspherical local structures cause an

在文檔中 液態水的定向自相關函數：SPC/E模型 (頁 29-0)