Introduction

Orientational correlation function

A. Orientational motion of a rigid molecule

For a general case, the principal moments of inertia of a rigid, anisotropic molecule are given by three different values Ix≠Iy≠Iz. The three unit vectors

eˆ ,

^x

eˆ and

^y

eˆ

z

along each principal axis of the molecule are mutually orthogonal and form a body-frame coordinate system of the molecule. In the space frame, the evolution of the three unit vectors is governed by the equations of motion

)

where the subscript index

μ

can be x, y or z. The instantaneous angular velocity of the molecule,

ω

^r(t), has three components

ω

_x(t)_,

ω

_y(t) _and

ω

_z

_{(t )}

at an instant t with their rotation axes lying along the principal x-, y- and z-axes in the body frame, respectively and, observed in the space frame,

ω

^r(t) is given as

as

B. Orientational correlation function of a rigid-molecule liquid

We consider the orientational dynamics in a liquid of such anisotropic molecules.

The single-molecular orientational dynamics of the liquid is characterized by the l-th order orientational correlation function (OCF) ^Cμ^(l)

⁽

^t

⁾

, with

μ = x, y or z and l being

an integer, which is defined as

[

^{ˆ ( ) ˆ (0)}

]

where Pl

(x) is the l-th order Legendre polynomial and the angular brackets denote an

ensemble average for the liquid at temperature T. For l=1, with P1

(x)=x,

)

For liquids at thermal equilibrium, the first time derivative ^C

&

μ⁽¹⁾

⁽

^t

⁾

is always zero.

where an integration by part in time has been used in the second equality. By using eq.

(II.2) for

e

ˆ t&_μ( ) and

e&

ˆ_μ(0) individually, ^C

&&

μ⁽¹⁾

⁽

^t

⁾

can be written as

where the ensemble averages for the correlations between unit vectors and those between the angular velocity components have been separated out due to the independence of the two kinds of correlations.

In the short-time regime, the correlation functions of the unit vectors and those of the angular velocity components, respectively, satisfy the following approximations

)

For the first approximation in eq. (II.4), if t is very small, the directions of the three unit vectors are not deviated much from their initial directions, so that

e

ˆ t_ν( ) at small

t is still almost orthogonal to other two unit vectors at t=0. The second approximation

in eq. (II.5) is resulted from the statistical independence among the three angular velocity components of molecular liquids at thermal equilibrium. Thus, with eqs. (II.4) and (II.5), ^C

&&

μ⁽¹⁾

⁽

^t

⁾

can be approximated as

where the subscript indices μ

≠

ν

≠

η are any combinations of

x ,

,

y z

.

By substituting the approximated ^C

&&

μ⁽¹⁾

⁽

^t

⁾

in eq. (II.6) into eq. (II.3),

C

_1,μ(

t

), with

z

y x ,

,

μ

= , for the three axes satisfy the combined integral equations

[ ^{( ) ( 0 ) ( ) ( ) ( 0 ) ( )} ]

Further, in the second cumulant approximation, we have

( ( ) )

With the similar method, we can prove that in the second cumulant approximation,

( ^{3 ( )} )

We define the normalized angular velocity autocorrelation function (AVAF) as

)

Using the equipartition theorem for each rotational degree of freedom of a molecule, we have

ω

η²(0) =

k

_B

T

/

I

η and, therefore, By substituting eq. (II.9) into the integrand in eq. (II.8),

⎥⎥ Fourier cosine transformation

∫

^∞Ω

(II.11) into eq. (II.10) and exchanging the integration order between t and

ω

, we obtain

C_μ in the short-time regime. This is a surprised result that, in the second

cumulant approximation, ^Cμ⁽¹⁾

⁽

^t

⁾

and ^Cμ⁽²⁾

⁽

^t

⁾

of the

μ

-th principal axis has nothing to with the power spectrum of the AVAF of the

μ

-th axis but is related with the power spectra of the AVAFs of the other two axes.

By expanding cos( t

ω

) in eq. (II.12) into Taylor series in the short-time limit such

and using the normalization of

D

_η(

ω

) and

D

_ν(

ω

),

ψ

_μ(t) in the lowest order in t is

This result means that although they interact with each other the interacting molecules at sufficiently short times rotate as free rotators so that the short-time behavior of

)

)

(

1

( t

C_μ is a Gaussian decay, with a decay time constant

τ

_μ only associated with the temperature of the molecular liquid. We can get

τ

_x= 0.06026ps,

τ

_y= 0.07949ps and

τ

z= 0.05734ps.

C.1 Instantaneous normal mode (INM) approximation of D

_η(

ω

)

Recently, the INM theory of liquids has been well developed and the short-time dynamics of liquids are accurately described by the theory. In this subsection, we consider the INM approximation of ^Cμ⁽¹⁾

⁽

^t

⁾

and ^Cμ⁽²⁾

⁽

^t

⁾

in the short-time limit.

The INMs of a liquid are referred to the eigenmodes of dynamic matrices of the liquid at an instant. For a liquid of N rigid molecules, the i-th molecule in a configuration is specified by

r

r_i ≡(

x

_i,

y

_i,

z

_i,

ξ

_ix,

ξ

_iy,

ξ

_iz)

, where the first three are the

center-of-mass coordinates of the molecule.

ξ

ix_,

ξ

_{iy and}

ξ

_iz are the rotational angles with respect to the instantaneous principal axes of the molecule in direction

) ˆ t(

e

_ix ,

e

ˆ t_iy( ) and

e

ˆ t_iz( ), respectively, A configuration of the molecular liquid is specified by the 6N-dimensional mass-weighted coordinates

R≡ { z

r₁,⋅ ⋅⋅,

z

r_N

}

, with

z

r≡(

m x

_i,

m y

_i,

m z

_i,

I

_x

ξ

_ix,

I

_y

ξ

_iy,

I

_z

ξ

_iz) ,

where m is the molecular mass and Ix, Iy and Iz are the three principal moments of inertia of the molecule.

At a configuration R0, the 6N×6N dynamic matrix D(R0) is the second derivatives of the total potential energy V(R) with respective to the mass-weighted coordinates given by the formula

0

where the particle indices, i and j, run from 1 to N and the indices for the degrees of freedom of a molecule, μ and ν, run from 1 to 6, with 1, 2 and 3 for the translation of the center of mass and 4, 5 and 6 for the rotations along the three principal axes. Each dynamic matrix can be diagonalized by an orthogonal transformation with the 6N×6N matrix U(R0) and has 6N eigenvalues given as

)

direction of the mass-weighted coordinate

z . Because of the normalization of each

^iμ eigenvector and the completeness of the 6N orthogonal eigenvectors, we have the following sum rules of the projection components

1

With the INM frequencies obtained at each configuration R0, the INM density of states (DOS) is defined as

( )

where the angular brackets represents an ensemble average of the liquid configurations. In principle,

D

^INM(

ω

) includes the contributions from both the translational and rotational degrees of freedom of all molecules. In terms of the eigenvector components

U

_α,iμ(

R

₀), we define the projection operators as

By the projection operators, the INM DOS contributed by the translational degrees of freedom is defined as

Similarly, the INM DOS contributed by the rotational degrees of freedom is given as

( )

the rotational INM DOS, contributed by the rotations about the x, y and z principal axes of a molecule, respectively. With the normalization conditions in eq. (II.14), we have

)

Since the configurations at which the dynamic matrices are evaluated may not correspond, in general, to the local minima of V(R), the eigenvalues

ω

α² _{are not}

guaranteed positive definite and the INM frequencies

ω

^α may thus be either real or imaginary. Therefore, the INM DOS

D

^INM(

ω

) generally consist of two lobes,

)

,^s(

ω

D

INM and

D

^INM,u(

λ

), which are the DOS of the stable-INM (denoted by the superscript s) with real frequencies

ω

^α and the unstable-INM (denoted by the superscript u) with imaginary frequencies

ω

_α =

i λ

_α, respectively. The normalization of

D

^INM(

ω

) yields a unit area under the curves of

D

^INM,s(

ω

)

and D

^INM,u(

λ

)

,

which are usually plotted along the positive and negative axes of ω respectively.

That is, has a bad divergence in the long-time limit due to the imaginary-frequency term in eq.

(II.16). In the stable-INM approximation in which the contribution due to the

imaginary-frequency INMs is neglected, Ω_η(t) is further simplified as imaginary-frequencey INMs are a small fraction in the total number of the INMs so that the difference between

D

_η(

ω

) and

D

~_η^INM_,^s(

ω

) is expected to be small.

INM-approximated effective spectrum given in eq. (II.17),

ψ

_μ(t) is approximated as

)

INM

(t

ψ

μ . Correspondingly, ^Cμ⁽¹⁾

⁽

^t

⁾

and ^Cμ⁽²⁾

⁽

^t

⁾

in the stable-INM approximation are, respectively, expressed as

[ ^{( )} ]

D. Long-time behavior of

^Cμ⁽¹⁾

⁽

^t

⁾ and

^Cμ⁽²⁾

⁽

^t

⁾

In the long time limit, the direction

e

ˆ t_μ( ) of a molecule in the liquid is generally

independent of its initial direction at t=0 and the OCF ^Cμ^(l)

⁽

^t

⁾

can be described by the Debye theory for the rotational dynamics of molecular liquids.

In the Debye theory for rotational relaxation,

p

_μ(

θ

,

φ

,

t

) is the two-dimensional probability to find the

μ

-th principal axis of a molecule pointing in the direction

( ) θ , φ

at time t , where

θ

and

φ

are the polar and azimuth angles of the spherical polar coordinates, respectively. With this probability distribution, the ensemble averaged ^Cμ^(l)

⁽

^t

⁾

is defined as

The probability distribution

p

_μ(

θ

,

φ

,

t

) satisfies the rotational diffusion equation in angular space

where ∇ is the angular part of the Laplacian operator in spherical coordinates and ²

D is the rotational diffusion coefficient for the

μ

μ

-th principal axis of the molecule.

With the initial condition that the probability distribution is a delta function at

=0

=

φ

θ

, the solution of the rotational diffusion equation is written as

t

With this solution, one can show that

)

a rotational-diffusion decay time ⁽⁾ =

[

( +1) _μ

]

⁻¹

τ

_D^lμ

l l D

. For l=1 and 2,

τ

_D⁽¹_μ⁾ =(2

D

_μ)⁻¹ and

τ

_D⁽²_μ⁾ =(6

D

_μ)⁻¹ so that the ratio

τ

_D⁽¹_μ⁾ /

τ

_D⁽²_μ⁾ of the rotational-diffusion decay times in the long time limit is estimated to be three by the Debye theory.

Chapter III

SPC/E Model, MD Simulation and Local Structures

A. The SPC/E model

Water molecular models have been developed in order to help discovering the structures of water. Using proper model can successfully predict the physical properties of liquid water. For MD simulations, several rigid models of water, for example, the ST2, TIP4P, SPC, SPC/E, etc., have been proposed. Among these rigid models, the TIP4P and the SPC/E models produce good agreements for the oxygen-oxygen radial distribution function [1], although they overestimate its first maximum. At ambient conditions, the internal energy predicted by most rigid models, including the ST3, TIP4P, SPC, and SPC/E, etc., are different from the experimental value within 5% [2][3]. For the diffusion coefficient, the TIP4P and the SPC models predict higher values, while within the simulation error, the SPC/E model gives the same value as experiment at T=298K [3].

In this thesis, the SPC/E model is used in the MD simulation. In this model, the interactions between two water molecules are expressed as

6 .

, 12

,

,j iO jO iO jO

i j i

ij

r

N r

M r

q

V

=

∑∑ q

+ −

α β α β

β α

Here,

r

^iα_,^jβ denotes the distance between the

α

^th _{site of}

i molecule and the

^th

β

^th

site of

j molecule.

^th

q is the charge on the

iα

α

^th site of

i molecule. M and N

^th are model specified parameters which are listed in Table I. With this model, many structural and dynamical properties can be calculated.

B. MD simulations

MD simulations are widely used in studying the structural and dynamical properties of liquid. In this thesis, we have simulated 256 SPC/E water molecules in a cubic box with periodic boundary conditions at density

ρ

=1.0

g cm

³and temperature T=300K. The long-ranged Coulomb interactions between water molecules have been treated by the Edward sum rule. The Lennard-Jones part between two oxygen atoms in the SPC/E model is truncated at half of the box length and shifted upward to make the potential and its first derivative continuous at the cutoff. To integrate the equations of motion, we use the leap-frog algorithm [4] with time step of 1 fs. To make sure the molecules are rigid, the SHAKE algorithm is used.

Before starting MD simulations, we initiated the system at a face-centered cubic lattice structure with randomly distributed velocities which follow the Maxwell-Boltzmann distribution. In order to make sure that the states used in the further calculations equilibrium, we discarded the configurations at the first 300,000 time steps and, then, started to take the data for 1,000,000 configurations at every 1 steps.

C. Radial distribution functions

With the configurations generated by MD simulations, the structures of liquid water are examined by the pair correlation functions

g

_ab(r), where ^a and ^b are

the oxygen or the hydrogen atom. We calculated the oxygen-oxygen, oxygen-hydrogen, and hydrogen-hydrogen radial distribution functions, which are denoted as

g

_OO(r),

g

_OH(r),^gHH

^{(r )}

, respectively.

The results are shown in Fig. 2. In

g

_OO(r), the first maximum occurs at r = 2.75Ǻ with a magnitude 3.06, which well agrees with other results obtained by molecule dynamic simulations [5][6]. In

g

_OH(r), there are two peaks at r = 1.76 Ǻ and 3.26Ǻ with almost equal magnitudes 1.6. In^gHH

^{(r )}

, there are also two strong peaks at r = 2.39Ǻ and 3.87Ǻ, however, the magnitude of the first peak is 1.35 and that of the second peak is 1.17, which is smaller than that of the first peak.

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

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

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