2. Fundamentals of Molecular Dynamics Method

(1)

2.1 Equation of Motion and Potential 2.2 Examples of Potential

2.2.1 Lennard-Jones Potential 2.2.2 Potential Models for Water

2.2.3 Potential Models for Carbon and Silicon 2.3 Integration of Equation of Motion

2.4 Boundary Condition

2.5 Initial Condition, Temperature and Pressure Control 2. Fundamentals of Molecular Dynamics Method 2. Fundamentals of Molecular Dynamics Method

Φ

−∇

=

i i i

dt m d r F

2 2

Newton's equation of motion

) (

_ij

¦¦

_>

=

i ji

φ r Φ

Pair Potential Approximation

Equation of Motion and Potential Equation of Motion and Potential

Approximation of Schrödinger Equation

∂ ψ

∂ψ H i t =

Pair Potential Approximation Pair Potential Approximation

2 2

dt m

_i

d

ⁱ

j ij i

F r F = ¦ =

Newton’s Equation

Small Droplets Small Droplets

Only 256 molecules

864 molecules

Lennard-Jones (12-6) Potential (1)

σ 2σ

-ε

Distance between Molecules r Repulsion

Attraction

( )

φ ε σ σ

r = § r r

©¨ ·

¹¸ − §

©¨ ·

¹¸

ª

¬«

«

º

¼» 4 »

1 2 6

Potential

r

σ

62

Lennard-Jones Potential (2) Lennard-Jones Potential (2)

Parameters for inert molecules.

σ [nm] ε [J] ε/k_B [K]

Ne 0.274 0.50×10^-21 36.2

Ar 0.340 1.67×10^-21 121

Kr 0.365 2.25×10^-21 163

Xe 0.398 3.20×10^-21 232

( ) ( ) ( )

_»^»

¼ º

««

¬

ª − −

¸¸¹·

¨¨©§

−

°¿+

°¾

½

°¯

°®

¸

¹

¨ ·

©

−§

¸¹

¨ ·

©

= § ⁻ ⁻ ⁻ C^*⁻⁶

*12 C 2

C

*6 C

*12 C 6 12

4 7 3

6

4 r r

r r r r r

r ε σr σ

φ

Cut-Off of potential: rc= 2.5～5.5σ

(2)

Property Reduced Form

Length r* = r/σ

Time t* = t/τ = t(ε/mσ²)^1/2 Temperature T* = k_BT/ε

Force f* = fσ/ε

Energy φ* = φ/ε

Pressure P* = Pσ³/ε Number density N* = Nσ³ Density ρ* = σ³ρ/m Surface Tension γ* = γσ²/ε

Lennard-Jones Potential (3) Lennard-Jones Potential (3)

Non-dimensional Form for L-J System

ε σ

= τ m/

= ε N E E

Lennard-Jones Potential (4) Lennard-Jones Potential (4)

Phase Diagram T-ρ

0 0.2 0.4 0.6 0.8

0.6 0.8 1 1.2 1.4 1.6

Number density ρ* = ρσ³

Temperature T* = kT/ε

0 –0.1 –0.2 0.1

0.2 0.3 0.4 0.5

–0.3 –0.4

–0.5

Lennard-Jones Potential (5) Lennard-Jones Potential (5)

Phase Diagram p-v

0 10 20

0 0.1 0.2

Volume v* = 1/ρ*

Pressure p* 1.35

1.4 1.5

1.6

1.3

1.2 1.1

1.0

Small Droplets Small Droplets

Only 256 molecules

864 molecules

Potential for Water (BNS, ST2) Potential for Water (BNS, ST2)

r_OH rOM

O

+qH

-q_M

∠HOH

¦¦

»+

»

¼ º

««

¬ ª

¸¸¹·

¨¨©§

¸¸¹ −

¨¨© ·

= §

i j ij

j i

r q R q

R S

R ¹² ₀

6

12 OO 12

12 OO OO 2 1

12( , ) 4 ( ) 4

πε σ

ε σ

φ R R

BNS potential

by Ben-Naim and Stillinger (1972)

ST2 potential

by Stillinger and Rahman (1974)

( )

⁼ ^°

=

∠HOH 2cos⁻¹1/ 3 109.47

Potential for Water (SPC, SPC/E) Potential for Water (SPC, SPC/E)

¦¦

»+

»

¼ º

««

¬ ª

¸¸¹·

¨¨©§

¸¸¹ −

¨¨© ·

= §

i j ij

j i

r q q R

R ₀

6

12 OO 12

12 OO OO 2 1

12( , ) 4 4

πε σ

ε σ

φ R R

SPC potential

by Berendsen et al. (1981)

SPC/E potential

by Berendsen et al. (1987)

-q_O

+q_H

∠HOH rOH

( )

⁼ ^°

=

∠HOH 2cos⁻¹1/ 3 109.47

(3)

Potential for Water (TIP4P) Potential for Water (TIP4P)

¦¦

»+

»

¼ º

««

¬ ª

¸¸¹·

¨¨©§

¸¸¹ −

¨¨© ·

= §

i j ij

j i

r q q R

R ₀

6

12 OO 12

12 OO OO 2 1

12( , ) 4 4

πε σ

ε σ

φ R R

TIP4P potential

by Jorgensen et al. (1983)

°

=

∠HOH 104.52 r_OM

O

+q_H

+qH

∠HOH -q_M r_OH

Potential for Water (MCY, CC) Potential for Water (MCY, CC)

MCY potential

by Matsuoka et al. (1976)

CC potential

by Carravetta & Clementi (1984)

°

=

∠HOH 104.52 Fitting to ab initio

r₂₆

O6

H₁

H₂ H₃

H₄ -q_M

-qM

+q_H

+q_H O₅

r₁₄

r₅₆ r16

[ ]

[

êxp(êxp( ⁾⁾ êxp(êxp( ⁾⁾ êxp(êxp( ⁾⁾ êxp(êxp( ⁾⁾

]

) exp(

) 4 exp(

) , (

45 4 35 4 26 4 16 4 4

45 3 35 3 26 3 16 3 3

24 2 23 2 14 2 13 2 2

56 1 1 0 2

1 12

r b r

b r

b a

r b r

b r

b a

r b r

b r

b a

r b r a

q q

i j ij

j i

− +

−

− +

=

¦¦

_πε

φ R R

Potential for Water (Comparison) Potential for Water (Comparison)

ST2 SPC/E TIP4P CC

rOH [nm] 0.100 0.100 0.095 72 0.095 72

∠HOH [°] 109.47 109.47 104.52 104.52

σOO [nm] 0.310 0.316 6 0.315 4 N/A

εOO ×10^-21 [J] 0.526 05 1.079 7 1.077 2 N/A rOM

qH

a [nm]

[C]

0.08 0.235 7 e

0 0.423 8 e

0.015 0.52 e

0.024 994 0.185 59 e qM [C] -0.235 7 e -0.847 6 e -1.04 e -0.371 18 e

aCharge of electron e = 1.60219×10^-19 C

SPC/E, TIP4P, CC are currently used

Potential for Water (Polarizable & Flexible) Potential for Water (Polarizable & Flexible)

c/f Dipole moment:

Isolated water: 1.85 D SPC/E: 2.351 D Rigid potentials are EffectivePair Potential

Polarizable Potentials

-q_O

+qH

+q_H

Harmonic or Morse Springs Dang & Pettitt (1987)

Anderson et al. (1987) TIP4P: Dang & Chang (1997) MCY: Niesar et al. (1990) SPC: Bernardo et al. (1994) Flexible Potentials

Droplet of Water

Droplet of Water Potential for Covalent System (C, Si) Potential for Covalent System (C, Si)

{ }

¦ ¦

_< ⁻

=

i jij

ij ij ij

ij V r b V r

r f

) (

A

* R

C( ) ( ) ( )

Φ

( )

{

^e

}

e C

R exp 2

) 1 ( )

( Sr R

S r D f r

V − −

= − β

( )

{

^e

}

e C

A exp 2/

) 1 ( )

( S r R

S S r D f r

V − −

= − β

( )

°°

¯

°°®

+

>

+

<

»¼ −

«¬ º

ª −

−

<

=

D R r

D R r D R D R r

D R r r

f

0 / ) 2( 2sin 1 2 1

1 )

C(

π

2

* ij ji

ij

b

b b +

=

δ

θ

−

≠ ¸¸

¹

·

¨¨

©

§

¿¾

½

¯® +

=

¦

n

ijk j i k

ik C n

ij a f r g

b 1 ( ) ( )

) , (

2 2

2 2 2

) cos 1 (

)

(θ = + − + − θ

h d

c d g c

Tersoff (1988, 1989), Brenner (1990)

(4)

Potential for Covalent System (C, Si) Potential for Covalent System (C, Si)

Tersoff’s Silicon

0.15 0.2 0.25 0.3

–4 –2 0 2 4

Distance r_ij [nm]

Potential Energyφij [eV]

(R_e, D_e) θijk=45°

θijk=90°

θijk=126.7°

and 2–body θijk=180°

R

rij θ_ijk

i

j k

{

⁽ ⁾ ⁽ ⁾

}

) (

'= f_C r_ij V_R r_ij −b_ijV_A r_ij

φ

Potential for Covalent System (C, Si) Potential for Covalent System (C, Si)

Tersoff (Si) Tersoff (C) Brenner (C) De [eV] 2.6660 5.1644 6.325

Re [nm] 0.2295 0.1447 0.1315

S 1.4316 1.5769 1.29

β [nm^-1] 14.656 19.640 1.5 A 1.1000×10^-6 1.5724×10^-7 1.1304×10^-2

N 7.8734×10^-1 7.2751×10^-1 1

δ 1/(2n) 1/(2n) 0.80469

C 1.0039×10⁵ 3.8049×10⁴ 19

D 1.6217×10¹ 4.384 2.5

H -5.9825×10^-1 -5.7058×10^-1 -1

R [nm] 0.285 0.195 0.185

D [nm] 0.015 0.015 0.015

Potential for Covalent System (C, Si) Potential for Covalent System (C, Si)

Example: Brenner Carbon (modified)

Integration of Newton’s Equation Integration of Newton’s Equation

( ) i( ) (i ) ( ) i i

it t 2rt rt t t²F(t) m

r +∆ = − −∆ +∆

( )^t {i(^t ^t) (i^t ^t)} ^t

i =r +∆ −r −∆ 2∆

v

( )

i i i

i m

t t t t

t t F

v

v ∆ ∆¸+∆

¹

¨ ·

©

§ −

¸=

¹

¨ ·

©

§ +

2 2

( ) ( ) ¸

¹

¨ ·

©

§ + +

=

+ 2

t t t t t

t i i

i

∆ ∆

∆ r v

r Order of ∆t

0.005 τor 10 fs with argon 0.5 fs for covalent Carbon Verlet’s Method

Leap Flog Method (Modified Verlet)

Gear’s predictor-corrector method is also sometimes used

Boundary Condition (Periodic) Boundary Condition (Periodic)

Basic Cell

Replica

Replica Replica

Potential must be Cut-Off at L/2

Ewald sum method for Coulomb Term

Boundary Condition (Gas) Boundary Condition (Gas)

Mirror Boundary

=Simple Reflection Example

Periodic Boundary

(5)

Boundary Condition (Solid Wall) Boundary Condition (Solid Wall)

°¿

°¾

½

°¯

°®

¸

¹

¨ ·

©

− §

¸¹

¨ ·

©

= §

3 INT 9 3 INT INT INT S

S 2 15

45 ) 2

(z π mρ ε σ σz σz

Φ

One-Dimensional Function by Bulk Integration

Φ(z)

Boundary Condition (Solid Wall) Boundary Condition (Solid Wall)

°¿

°¾

½

°¯

°®

¸

¹

¨ ·

©

− §

¸¹

¨ ·

©

= §

4 INT 10 INT 2

0 2 INT

INT 2 5

15 3 ) 4

(z πε Rσ σz σz

Φ

0 2 4 6 8

-4 -2 0 2 [×10-21]

Integrated Potential Lennard-Jones

εINT=0.575x10-21 J σINT=3.085 Distance from surface [ ]

Potential energy[J]

Å Å

One-Dimensional Function by Surface Integration

Boundary Condition (Phantom Method)

Phantom molecules

Fixed molecules

move

t T k

s kg

b

F= ∆

×

= ⁻

σ α α

2

/ 10 184 .

5 ¹²

σF vertical：2k horizontal：0.5k vertical：2k horizontal：3.5k

α

F

Boundary Condition (Phantom Method)

Phantom (Not Shown)

Initial Condition Initial Condition

m T k v = 3

_B _C

/

Constant velocity for random direction or Maxwell-Boltzmann distribution

Make Sure Potential is not Too Large

Temperature Control Temperature Control

/ C

' v T T v_i= _i

dt md dt

m_id ⁱ _i _i rⁱ r=F−ζ

2 2

( )

Q E E dt

d k k

2 − 0

ζ =

Anderson method [Anderson (1980)]

Velocity Scaling

Nosé-Hoover Thermostat [Nosé(1984), Hoover(1985)]

Replace Velocity of Randomly Selected Molecule to Maxwell-Boltzmann Distribution

(6)

Pressure & Stress Control Pressure & Stress Control

Parrinello and Rahman (1980, 1981) Andersen (1980)

Berendsen et al. (1984)

Extension of Anderson: Change Shape of Box

P

C P t

P dt dP/ =( − )/

r r'=χ¹^/³

(^P ^P)

t t

C p

T −

−

= β ∆ χ 1

Change Box Size as if Piston is Connected