Weinan E. Principles of Multiscale Modeling

Подождите немного. Документ загружается.

4.2. MOLECULAR DYNAMICS 167

4.2.2 Equilibrium states and ensembles

Consider the solutions of Hamilton’s equation (4.2.2) with random initial condition.

Denote by ρ(z

, ··· , z

, t) the probability density of the particles at time t, where z

, p

). It is easy to see that ρ satisﬁes Liouville’s equation:

∂

ρ +

i=1

∇

· (V

ρ) = 0 (4.2.16)

where V

= (∇

H, −∇

H), ∇

= (∇

, ∇

Given a system of interacting particles (atoms or molecules, etc), the ﬁrst question

we should think about is their equilibirium states. These are stationary solutions of the

Liouville equation (4.2.16) that satisfy the condition of detailed balance:

∇

· (V

ρ) = 0

for each i = 1, ··· , N. It is obivous that probability densities of the form:

f(H)

are equilibrium states. Here f is a nonnegative function, and Z

is the normalization

constant, also called the partition function:

f(H)dz

···dz

Of particular importance are the following special cases:

1. The microcanonical ensemble: f (H) = δ(H − H

) for some constant H

2. The canonical ensemble: f(H) = e

−βH

where β = 1/k

T , k

is the Boltzmann

constant and T is the temperature.

The particular physical setting determines the particular equilibrium distribution that

is realized. In the context of molecular dynamics, the microcanonical ensemble is realized

through the NVE ensemble, in which the number of particles N, the volume V accessible

to th e particles and the total energy is ﬁxed. This models a closed system with no

exchange of mass, momentum or energy with the environment. Similarly, the canonical

ensemble is realized through the NVT ensemble in which instead of the total energy, the

168 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

temperature of the system is ﬁxed. Recall that temperature is deﬁned as the average

kinetic energy per degree of freedom due to the ﬂuctuating part of the particle velocity:

Assuming that the average velocity is zero:

j=1

= 0 (4.2.17)

where v

= p

, then

j=1

T (4.2.18)

The NVE ensemble is obtained when one solves directly the Hamilton’s equations

(4.2.2) in a box with, e.g. reﬂection bound ary condition. The value of the Hamiltonian

is conserved in this case. However, to obtain the NVT ensemble in which temperature

rather than the total energy of the system is ﬁxed, we need to modify Hamilton’s equation.

There are several empirical ways of d oing this. The simplest idea is the isokinetic energy

trick (also called velocity rescaling). The idea is very straightforward for a discretized

form of (4.2.2): At each time step, one rescales the velocities v

simultaneously by a

factor such that, after rescaling,

j=1

T (4.2.19)

holds, where T is the desired temperature. This does give the right average temperature,

but it does not capture the temperature ﬂuctuations correctly. Indeed in this case the

quantity

T =

j=1

does not ﬂuctuate. This is not the case for the true

canonical ensemble. In fact, one can show that [18]

T − T )

hT i

(4.2.20)

i.e. the ﬂuctuation in

T is of the order N

−

, which only vanishes in the limit when

N → ∞.

Other popular ways of generating the NVT ensemble include:

1. The Anderson thermostat: At random times, a random particle is picked and its

velocity is changed to a random vector from the Maxwell-Boltzmann distribution:

ρ(v

, v

, ··· , v

) ∼ e

−β

(4.2.21)

4.2. MOLECULAR DYNAMICS 169

2. Replacing Hamilton’s equation by the Langevin equation

= −

∂V

∂y

− γ

+ σ

, j = 1, ··· , N, (4.2.22)

where γ

is the friction coeﬃcient for the j-th particle, σ

is the noise amplitude:

T γ

(4.2.23)

{

, ··· ,

} are independent white noises.

3. The N´ose-Hoover thermostat. This is a way of producing the canonical ensemble by

adding an additional degree of freedom to the system. The eﬀect of this additional

degree of freedom is to supply energy or take energy away from the system according

to the instantaneous thermal energy of the system [18].

In general, these thermostats produce correctly the equilibrium properties of the system.

But it is not clear that dynamic quantities are produced correctly [18].

4.2.3 The elastic continuum limit – the Cauchy-Born rule

From the viewpoint of multiscale, multi-physics modeling, it is natural to ask how

the atomistic models are related to the continuum models discussed in the last section.

We will address this question from three diﬀerent but related angles. In this subsection,

we will discuss the elastic continuum limit at zero temperature. In the next subsection,

we will discuss the relation between the continuum models and the atomistic models at

the level of the conservation laws. Finally, we will discuss how to compute macroscopic

transport coeﬃcients using atomistic models.

At zero temperature T = 0, the canonical ensemble becomes a degenerate distribution

that is concentrated at the minimizers of the total potential energy:

min V (y

, . . . , y

). (4.2.24)

This special case is quite helpful for understanding many issues in solids.

It is well-known that at low temperature, solids are crystals whose atoms reside

on some lattices. Common lattice structures include body-centered cubic (BCC), face-

centered cubic (FCC), diamond lattice, hexagonal closed packed lattice (HCP), etc [32].

For example, at room temperature and atomspheric pressure, iron (Fe) exists in BCC

170 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

lattice, aluminum (Al) exists in FCC lattice, and silicon (Si) exists in diamond lattice.

Lattices are divided into two types: simple lattices and complex lattices.

Simple lattices. Simple lattices are also called Bravais lattices. Their lattice points

take the form:

L({e

}, o) = {x : x =

i=1

+ o ν

, are integers}, (4.2.25)

where {e

}

i=1

is a set of basis vectors, d is the dimension, and o is a particular lattice

site. Note that the basis vectors are not unique.

For the examples listed above, BCC and FCC are simple lattices. For the FCC lattice,

one set of basis vectors are

(0, 1, 1), e

(1, 0, 1), e

(1, 1, 0).

For the BCC lattice, we may choose

(−1, 1, 1), e

(1, −1, 1), e

(1, 1, −1).

as the basis vectors. Here and in what follows, we use a to denote the lattice constant.

Another example of a simple lattice is the two-dimensional triangular lattice. Its basis

vectors can be chosen as:

= a(1, 0), e

= a(1/2,

√

3/2).

Figure 4.4: The graphene sheet as an example of a complex lattice. The right panel

illustrates a unit cell (courtesy of Pingbing Ming).

Complex lattices. Crystal lattices that are not simple lattices are called complex

lattices. It can be shown that every complex lattice is a union of simple lattices, each of

4.2. MOLECULAR DYNAMICS 171

which is a shift of any other: i.e., the lattice points can be expressed in the form:

L = L({e

}, o) ∪ L({e

}, o + p

) ∪ ··· L({e

}, o + p

)

Here p

, ··· , p

are the shift vectors. For example, the two dimensional hexagonal lattice

with lattice constant a can be regarded as the union of two triangular lattices with shift

vector p

= a(−1/2, −

√

3/6). The diamond lattice is made up of two interpenetrating

FCC lattices with shift vector p

= a/4(1, 1, 1). The HCP lattice is obtained by stacking

two simple hexagonal lattices with the shift vector p

= a(1/2,

√

3/6,

√

6/3). Some

solids consist of more than one species of atoms. Sodium-chloride (NaCl), for example,

has equal numbers of sodium and chloride ions placed at alternating sites of a simple

cubic lattice. This can be viewed as the union of two FCC lattices: one for the sodium

ions and one for the chloride ions.

The picture described above is only for single crystals and it neglects lattice defects

such as stacking faults, vacancies and dislocations. Most solids are polycrystals, i.e.

unions of single crystal grains with diﬀerent orientations. Nevertheless, it is certainly

true that the basic building blocks of solids are crystal lattices.

Turning now to the physical models, we have discussed two approaches for analyzing

the deformation of solids: The continuum theory discussed in the last section and the

atomistic model. It is natural to ask how these two approaches are related to each

other. From the viewpoint of atomistic modeling, the relevant issue is the macroscopic

continuum limit of the atomistic models. From the viewpoint of the continuum theory,

the relevant issue is the microscopic foundation of the continuum models.

Recall that in continuum elasticity theory, the potential energy of an elastic solid is

given by

Ω

W (∇u)dx (4.2.26)

where W is the stored energy density. In the last section, we discussed ways of obtaining

empirical expressions for W . We now ask: Can we derive W from an atomistic model?

This question was considered already by Cauchy, who derived atomistic expressions

for the elastic moduli. Cauchy’s work was extend ed by Born who also considered complex

lattices [4]. For crystalline solids, the Cauchy-Born rule establishes an expression for W

in terms of an atomistic model.

Let us consider ﬁrst the case when the crystal is a simple lattice. In this case, given

a value of the deformation gradient A, consider a homogeneously deformed crystal with

172 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

deformation gradient A, i.e. the positions of the atoms in the deformed conﬁguration are

given by y

= x

+ Ax

for all j. Knowing the positions of the atoms in the deformed

crystal, we can compute the energy density of the deformed crystal using the atomistic

model. This gives the value of W

(A).

Figure 4.5: Illustration of the Cauchy-Born rule: W (A) is equal to the energy density of

a homogeneously deformed crystal obtained through the relation y = x + Ax (courtesy

of Pingbing Ming).

Next consider the case when the crystal is a complex lattice. In this case, to deﬁne

the conﬁguration of a homogeneously deformed crystal, we not only have to specify

the deformation gradient A, we also need to specify the shift vectors p

, ··· , p

. Once

they are speciﬁed, we can compute the energy density (energy per unit cell) of the

homogeneously deformed crystal using the atomistic mo d el. We denote this value by

w(A, p

, ···p

). We then deﬁne:

(A) = min

,···,p

w(A, p

, ··· , p

) (4.2.27)

To see where the Cauchy-Born rule comes from, consider the case of a simple lattice

with atoms interacting via a two-body potential:

V =

i6=j

(|y

− y

|) . (4.2.28)

Cauchy’s idea is to assume that the deformation of the atoms follows a smooth displace-

ment vector ﬁeld u, i.e. [16]

= x

+ u(x

), (4.2.29)

4.2. MOLECULAR DYNAMICS 173

Then the total energy of the system is given by:

V =

i6=j

(|(x

+ u(x

)) − (x

+ u(x

))|), (4.2.30)

i6=j

(|x

− x

+ ∇u(x

)(x

− x

) + O(|x

− x

)|).

To leading order, we have

V ≈

i6=j

(|(I + ∇u(x

))(x

− x

)|) =

W (∇u(x

))V

≈

W (∇u(x))dx (4.2.31)

where

W (A) =

i6=0

(|(I + A)(x

)|)

and V

is the volume of the unit cell of the undeformed lattice. This is nothing but the

Cauchy-Born rule. In one dimension with x

= ka, (4.2.31) becomes

V =



1 +

)



(4.2.32)



)



a (4.2.33)

≈



(x)



dx,

where

W (A) =

((1 + A)ka) (4.2.34)

Calculations of this type have been carried out for very general cases in [3]. In

essence, such a calculation demonstrates that to leading order, the atomistic model can

be regarded as the Riemann sum of the nonlinear elasticity model where the stored energy

density is given by the Cauchy-Born rule. The validity of the Cauchy-Born rule at zero

temperature has been studied in great detail in [14] for the static case and in [15] for the

dynamic case. The essence of this analysis is to ﬁnd out under what conditions, Cauchy’s

assumption that the atoms follow a smooth displacement ﬁeld, is valid. These conditions

are shown be in the form of stability conditions for the phonons as well as the elastic

stiﬀness tensor. These conclusions seem to be sharp and quite general.

174 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

We will refer to the nonlinear elasticity model in which W = W

as the Cauchy-Born

elasticity model.

So far we have only discussed the highly idealized situation when the temperature

is zero. The extension of the Cauchy-Born construction to ﬁnite temperature is still

somewhat of an open problem.

4.2.4 Non-equilibrium theory

From a statistical viewpoint, the basic object for the non-equilibrium theory is Liou-

ville’s equation. It contains all the statistical information we need. However, Liouville’s

equation is a very complicated object since it describes the dynamics of a function, the

probability density, of many variables. In the following, we will try to make a more

direct link with non-equilibrium continuum theory. We will do so by deriving, from the

equations of molecular dynamics, conservation laws similar to the ones in continuum the-

ory. We will restrict ourselves to the case when the inter-atomic potential is a two-body

potential:

V (y

, . . . , y

) =

i6=j

(|y

− y

|) (4.2.35)

Given the output from molecular dynamics, {(y

(t), v

(t)), j = 1, ··· , N}, we deﬁne

the empirical mass, momentum and energy densities as follows:

˜ρ(y, t) =

δ (y − y

(t)) , (4.2.36)

m(y, t) =

(t)δ (y − y

(t)) , (4.2.37)

E(y , t) =

(t)|

δ (y − y

(t)) (4.2.38)

i6=j

(|y

(t) − y

(t)|)

δ (y − y

(t)) .

Here δ(·) is the Delta function. Let f

= −∇

(|y

− y

|), we can deﬁne the analog of

4.2. MOLECULAR DYNAMICS 175

the stress tensor and energy ﬂux:

˜σ =

(t) ⊗ v

(t)δ(y − y

(t)) (4.2.39)

j6=i

(t) − y

(t)) ⊗ f

(t)

δ (y − (y

(t) + λ(y

(t) − y

(t)))) dλ,

J =

(t)

j6=i

(|y

(t) − y

(t)|)

δ (y − y

(t)) (4.2.40)

j6=i

(t) + v

(t)) · f

(t) ⊗ (y

(t) − y

(t))

δ (y − (y

(t) + λ(y

(t) − y

(t)))) dλ.

(4.2.39) is the well-known Irving-Kirkwood formula [22]. It is a straightforward but

tedious calculation to show that the following equations hold as a consequence of the

MD equation (4.2.1):

∂

˜ρ + ∇

m = 0, (4.2.41)

∂

m + ∇

· ˜σ = 0, (4.2.42)

∂

E + ∇

J = 0 (4.2.43)

These are the analog of the conservation laws in continuum theory.

Even though molecular dynamics naturally uses the Eulerian coordinates, these for-

mulas also have th eir analogs in Lagrangian coordinates. Denote x

= x

(0), which can

be thought of as being the Lagrangian coordinate of the j-th atom, and let [30]:

˜ρ(x, t) =

δ(x − x

), (4.2.44)

m(x, t) =

(t)δ(x − x

), (4.2.45)

E(x, t) =

(t)|

δ(x − x

) +

i6=j

(|x

− x



x − x



. (4.2.46)

˜σ =

j6=i

− x

) ⊗ f

(t)



x −



+ λ(x

− x

)



dλ, (4.2.47)

176 CHAPTER 4. THE HIERARCHY OF PHYSICAL MODELS

J =

j6=i

(t) + v

(t)) · f

(t) ⊗ (x

− x

)



x −



+ λ(x

− x

)



dλ (4.2.48)

we then have [30, 41]

∂

m + ∇

· ˜σ = 0, (4.2.49)

∂

E + ∇

J = 0 (4.2.50)

In principle, one way of deriving the continuum models is to perform averaging on the

conservation laws we just discussed. In practice, we face the closure problem: In order

to calculate the averaged ﬂuxes, we need to know many-particle probability density of

the atoms. We get a closed form model if this probability density can be expressed in

terms of the averaged conserved quantities (the mass density, averaged momentum and

energy densities). However, this is usually not the case. Therefore approximations have

to be made in order to solve the closure problem. This point will be made more clearly

in the next section.

4.2.5 Linear response theory and the Green-Kubo formula

The connection between the continuum and atomistic models established above is

quite general but it may not be explicit enough to oﬀer signiﬁcant physical insight. Here

we consider the linear response regime and establish a connection between the transport

coeﬃcients, such as diﬀusion constants and viscosity coeﬃcients, in terms of microscopic

quantities. We can always put a system into the linear resp onse regime by limiting the

size of the external driving force.

First, let us consider th e diﬀusion of a tagged particle in a ﬂuid. We assume that its

probability distribution ρ satisﬁes the diﬀusion equation:

∂

ρ(x, t) = D∆ρ(x, t) (4.2.51)

Assume that the particle is at the origin initially. Since h|x(t)|

i =

|x|

ρ(x, t)dx, we