1.021 , 3.021, 10.333, 22.00 I ntroduc tion to Modeling and Simulation

Spring 2011

Part I – C ontinuum and partic le me thods

Property calculation I

Lecture 3

Markus J. Buehler

Laboratory for Atomistic and Molecular Mechanics Department of Civil and Environmental Engineering Massachusetts Institute of Technology

Content overview

I. Particle and continuum me thods

1. Atoms, molecul e s, chemistry

2. Continuum modeling approac hes and solution approaches

3. Statistical mechanics

4. Molecular dynamics, Monte Carlo

5. Visualization and data analysis

6. Mechanical proper ties – applic ation: how things fail (and how to prevent it)

7. Multi-scale modeling par adigm

8. Biological systems (simulation in biophysics) – h ow proteins work and how to model them

II. Quantum mechanical methods

1. It’s A Q uantum World: T he Theory of Quantum Mechanics

2. Quantum Mechanics: Practice Makes Perfect

3. The Many-Body Problem: Fr om Many-Body to Single- Particle

4. Quantum modeling of materials

5. From Atoms to Solids

6. Basic pr operties of mater i als

7. Advanced proper ties of materials

8. What else can we do?

Lectures 2-13

Lectures 14-26

Lecture 3: Property calculation I

Outline:

1. Atomis tic model of diffusion

2. Computing power: A perspective

3. How to calculate properties from atomistic simulation

3.1 Thermodynamical ensembles: Micro and macro

3.2 How to calc ulate properties from atomistic simulation

3.3 How to solv e the equations

3.4 Ergodic hypothesis

Goal of today’s lecture:

 Exploit Mean Square Displacement functi on to identify diffusivity, as well as material state & structure

 Provide rigorous basis for property ca lculation from molec u lar dynamic s simulation re sults (sta tistical mechanics)

Additional Reading

Books:

M.J. Buehler (2008): “Atomistic Modeling of Materials Failure” Allen and Tildesley: “Compu ter simulation of liquids” ( classic )

D. C. Rapaport (1996): “The Art of Molecular Dy namics Simulation”

D. Frenkel, B. Smit (2001): “Und erstanding Molecular Simulation”

J.M. Haile (Wiley, 1992), “Mol ecular dynamics simulation”

1. Atomistic model of diffusion

How to build an atomistic bottom-up model to describe the physical phenomena of diffusion?

Introduce: Mean Square Displacement

Recall: Diffusion

 Particles move from a domain with high concentration to an area of

low concentration

 Macroscopically , diffusion measured by change in concentration

 Microscopically , diffusion is process of spontaneous net movement of particles

Result of random motion of particles (“Brownian motion”)

High concentration

Low concentration

c = m/V = c ( x, t )

Ink droplet in water

hot cold

© s o ur c e u n k n o w n. Al l rights res e rv ed. Thi s co ntent i s ex cl u d ed from our Creati ve Co mmo ns lice n se . F o r mo r e in fo r m a t io n , s e e http:/ /o c w .mi t .edu /fai ru se .

Atomistic description

Back to the application of diffusion problem…

 Atomis tic description provides alternative way to predict D

 Simple solve equation of motion

 Follow the trajectory of an atom

 Relate the average distance as functi on of time from initial point to diffusivity

 Goal: Calculate how particle s move “randomly”, away from initial position

Pseudocode

Set particle positions (e.g. cr y s tal lattice) Assign initial velocities

For (all time steps):

Calculate force on each particle (subroutine) Move particle by time step  t

Save p a rticle po sition, velo city , acceleration

Save results Stop simulation

r ( t   t )  2 r ( t )  r ( t   t )  a ( t

)  t 2

 ...

a i 

f i / m

i 0 i 0 i 0 i 0

Positions at t 0

Positions at t 0 -  t

Accelerations at t 0

JAVA applet

Courtesy of t he Center for Polymer S t udi e s at Bosto n Uni v e r si ty. U se d wi th pe rmi s si o n .

URL http://polymer.bu.edu/java/java/LJ/index.html

Link atomistic trajectory with diffusion constant (1D)

Diffusion constant relates to the “abili ty” o f a particle to move a distance  x 2

 x 2

(from left to right) over a time  t

 x 2

D  p  t

Idea – U se MD simulation to measure squar e of dis p lacement from initial

pos ition of particles,

 r 2 ( t ) :

scalar pr oduct

 r 2 ( t )

 1   →

 →  0 )  2

 1    →

 →  0 )    →

 →  0 )  

r i ( t )

i

r i ( t

r i ( t )

i

r i ( t

r i ( t )

r i ( t

time 11

Link atomistic trajectory with diffusion constant (1D)

Diffusion constant relates to the “abili ty” o f a particle to move a distance  x 2

 x 2

(from left to right) over a time  t

 x 2

D  p  t

MD simulation: Measure square of displacement from initial position

of particles,  r 2 ( t ) :

 r 2 ( t )

 1   →

 →  0 )  2

 r 2

r i ( t )

i

r i ( t

t

Link atomistic trajectory with diffusion constant (1D)

Diffusion constant relates to the “abili ty” o f a particle to move a distance  x 2

 x 2

(from left to right) over a time  t

 x 2

D  p  t

MD simulation: Measure square of displacement from initial position

 r 2

of particles,  r 2 ( t ) and not  x 2 ( t ) ….

Replace

D  p

 x 2

 t

D  2  t

Factor 1/2 = no directionality in (equal probability to move forth or back)

Link atomistic trajectory with diffusion constant (1D)

MD simulation: Measure square of displacement from initial position

of particles,  r 2 ( t ) :

 r 2

2

D  R ~

2 t

t

Link atomistic trajectory with diffusion constant (2D/3D)

D  p

 x 2

 t

Higher dimensions

Since:

D  1 1

 r 2

Factor 1/2 = no directio nality in (forth/back)

Factor d = 1, 2, or 3 due to 1D, 2D, 3D (dimensionality)

2 dD  t ~  r 2

C = constant (does not affect D )

Example: MD simulation

Mean Square Displacement function

slope = D

D  1 2 d

lim d

t   d t

  r 2 

 r 2

1D=1, 2D=2, 3D=3

C

Courtesy of S i d Yip. Used with permissi on.

D  1 2 d

lim d

t   d t

 r 2

.. = average over all particles

Example molecular dynamics

 r 2

Particles Trajectories

Mean Square Displacement function

Average square of dis p lacement of all particles

 r 2 ( t )

 1   →

 →  0 )  2

r i ( t )

i

r i ( t

Courtesy of t he Center for Polymer S t udies at Bosto n Uni v e r si ty. U se d wi th permi s si o n . 17

Example calculation of diffusion coefficient

 r 2 ( t )

 1   r ( t )  r ( t

 0 )  2

i i

i

Position of

atom i at time t

Position of

atom i at time t=0

slope = D

 r 2

1D=1, 2D=2, 3D=3

18

Summary

 Molec ular dynamic s provides a po werful approach to relate the diffusion constant that appears in continuum models to atomistic trajectories

 Outlines multi-scale approach: Feed parameters from atomistic simulations to continuum models

Time scale

Co ntinuum model

MD

“Empirical”

or experimental parameter feeding

Summary

 Molec ular dynamic s provides a po werful approach to relate the diffusion constant that appears in continuum models to atomistic trajectories

 Outlines multi-scale approach: Feed parameters from atomistic simulations to continuum models

Time scale

Co ntinuum model

MD

“Empirical”

or experimental parameter feeding

MD modeling of crystals – s olid, liquid, gas phase

 Crystals: Regular, ordered structure

 The corresponding particle motions are small-amplitude vibrations about the lattice site, diffusive mo vements over a local region, and long free flights interrupted by a collis ion every now and then.

 Liquids: Partic les follow Brownian motion (collis i ons)

 Gas: Very long free paths

Image by MIT OpenCou r seWare. A f te r J. A. Barker and D. Hender s o n.

Example: MD simulation results

liquid solid

solid

Courtesy of S i d Yip. Used with permissi on.

Atomistic trajectory

Courtesy of S i d Yip. Used with permissi on.

Multi-scale simulation paradigm

Courtesy of Elsevier, Inc., http://www.scien c edirect. com . Used with p e rmi s sion.

2. Computing power: A perspective

Co u r t e sy El se vier , I n c ., http://www.scien cedirect.com . Used wit h p e rmi s sion.

Historical development of computer simulation

 Began as tool to exploit computin g machines developed during World War II

 MANIA C (1952) at Los Alamos us ed for computer simulations

 Metropolis, R o senbluth, Teller (1953 ): Metropolis Monte Carlo method

 Alder and Wainwright (Livermore National Lab, 1956/1957): dynamics of hard spheres

 Vineyard (Brookhaven 1959-60): dy namics of radiation damage in copper

 Rahman (Argonne 1964): liquid argon

 Application to more complex fluids (e.g. water) in 1970s

 Car and Parrinello (1985 and following): ab-initio MD

 Since 1980s: Many applications, including:

 Karplus, Goddard et al.: Applic ations to polymers/biopolymers, proteins since 1980s

 Applications to fracture since mid 1990s to 2000

 Other engineering applications (nanotechnology, e.g. CNTs, nanowires etc.) since mid 1990s-2000

3. How to calculate properties from atomistic simulation

A brief introduction to statistical mechanics

Molecular dynamics

Follow trajectories of atoms ( classical mechanics, Newton’s laws )

“Verlet central difference method”

r i ( t 0

  t ) 

 r i ( t 0

  t )  2 r i ( t 0

)  t

 a i

( t 0

)   t  2

 ...

a i 

f i / m

Positions at t 0 -  t

Positions at t 0

Accelerations at t 0

Property calculation: Introduction

Have:

→

x ( t ),

→

x ( t ),

→

x ( t )

“microscopic information”

Want:

 Thermodynamical properties (temperature, pressure, stress, strain, thermal conductivity, ..)

 State (gas, liquid, solid)

 …

(properties that can be measured in experiment!)

Goal: To develop a robust framework to calculate a range of “macroscale” properties from MD simulation studies (“microscale information”)

3.1 Thermodynamical ensembles: Micro and macro

Macroscopic vs. microscopic states

C 1

 C 2

T , p , V , N C 3

… C N

Same macroscopic state is represented by many different microscopic configurations C i 33

Definition: Ensemble

 Large number of copies of a system with specific features

 Each copy represents a possible microscopic state a macroscopic system might be in under thermodynamical constraints ( T, p, V, N ..)

Gibbs, 1878

Microscopic states

Microscopic states characterized by

r , p

r   →  ,

  → ˙ 

i  1 .. N

x i p

m i x i

p i

Microscopic states

Microscopic states characterized by

r , p

r   →  ,

  → ˙ 

i  1 .. N

x i p

m i x i

p i

Hamiltonian (sum of potential and kinetic energy = total energy) expressed in terms of these variables

H ( r , p )

 U ( r ) 

K ( p )

U ( r ) 

  i

( r )

K ( p ) 

 1 i

K  1 m v 2

i  1 .. N

i  1 .. N 2 m i

i 2 i i

Ensembles

Result of thermodynamical constraints, e.g. temperature, pressure…

Microcanonical Canonical

Isobaric-isothermal Grand canonical

NV E NV T

Np T TV 

chemical potential (e.g. concentration)

3.2 How to calculate properties from atomistic simulation

Link between statistical mechanics and thermodynamics

?????

Link between statistical mechanics and thermodynamics

Statistical mechanics

Macroscopic conditions (e.g. constant volume, temperature, number of particles…) translate to the microsco pic system as boundary conditions (constraints)

Macroscopic system : defined by extens ive variables, whic h are constant:

E.g. ( N,V,E )= NVE ensemble

Link between statistical mechanics and thermodynamics

Statistical mechanics

Macroscopic conditions (e.g. constant volume, temperature, number of particles…) translate to the microsco pic system as boundary conditions (constraints)

Macroscopic system : defined by extens ive variables, whic h are constant:

E.g. ( N,V,E )= NVE ensemble

The behav ior of the microscopic system is related to the macroscopic conditions. In other words, the distribution of microscopic states is related to the macrosco pic conditions.

To calculate macroscopic properties (via statistical mechanics) from microscopic information we need to kn ow the distribution of microscopic states (e.g. through a simulation)

Example: Physical realization of canonical ensemble ( NVT )

Heat bath (constant temperature)

Coupled to large system, allow energy exchange

NV T

“small” system embedded in “large” heat bath

Constant number of particles = N

Constant volume = V

Constant temperature = T

Macroscopic vs. microscopic states

Canonical ensemble



N , V , T

C 1 r 1 , p 1

C 2 r 2 , p 2

C 3 r 3 , p 3

… C N

r N , p N

Same macroscopic state is represented by many different microscopic configurations

Important issue to remember…

 A few slides ago:

“ To calculate macroscopic properties (via statistical mechanics) from microscopic information we need to know the distribution of microscopic states (e.g. through MD simulation) ”

 Therefore:

We can not (“never”) take a single measurement from a single microscopic state to relate to macroscopic properties

Micro-macro relation

C our tes y o f the Cente r fo r Po lyme r Stud ie s a t Bos t on U n ive r s i ty. U sed w i th pe rmis s i on.

T ( t )

Which to pick?

1 1 N → t

T ( t ) 

 m v 2 ( t )

3 Nk B

i i

i  1

Specific (individual) microscopic states are insufficient to relate to macroscopic properties

Averaging over the ensemble

 Rather than taking single measur ement, need to average over “all” microscopic states that represent t he corresponding macroscopic condition

 This averaging needs to be done in a su itable fashion, that is, we need to consider the specific d i stri bution of microscopic states (e.g. some microscopic states may be more likely than others)

What about trying this….

Property A 1 Property A 2 Property A 3

A

 1  A  A

 A 

C 1

r 1 , p 1

C 2

r 2 , p 2

C 3

r 3 , p 3

macro 3 1 2 3

46

Averaging over the ensemble

 Rather than taking single measur ement, need to average over “all” microscopic states that represent t he corresponding macroscopic condition

 This averaging needs to be done in a su itable fashion, that is, we need to consider the specific d i stri bution of microscopic states (e.g. some microscopic states may be more likely than others)

What about trying this….

Property A 1 Property A 2 Property A 3

A 1  A  A  A 

macro 3 1 2 3

C 1 C 2 C 3

Generally, NO!

47

Averaging over the ensemble

Property A 1 Property A 2 Property A 3

A

 1  A  A

 A 

C 1 C 2

macro 3 1 2 3

C 3

Instead, we must average with proper weights that represent the probability of a system in a particular microscopic state!

(I.e., not all microscopic states are equal)

A macro

  1 A 1

  2 A 2

  3 A 3 

 1 ( r 1 ,

p 1 ) A 1 ( r 1 ,

p 1 ) 

 2 ( r 2 ,

p 2 ) A 2 ( r 2 ,

p 2 ) 

 3 ( r 3 ,

p 3 ) A 3 ( r 3 ,

p 3 )

48

How to relate microscopic states to macroscopic variables?

A ( r , p )

Property due to specific microstate

• Ensemble average, obtained by integral over all microscopic states

• P roper weight

 ( p , r )

- depends on ensemble

A ( r ,

How to relate microscopic states to macroscopic variables?

p ) Property due to specific microstate

To measure an observable quantity from MD simulation we must express this observable as a function of the positions and linear momenta of the particles in the system, that is, r, p

Recall, microscopic states characterized by

r , p

r   → 

  → ˙ 

i  1 .. N

x i , p

m i x i

p i

How to relate microscopic states to macroscopic variables?

 A 

  A ( p , r )  ( p , r ) drdp

p r

Probability density distribution

 ( p , r ) 

1 exp  

H ( p , r ) 

Probability to find system

Q  k B T 

in state ( p,r )

Boltzmann constant

k B

Partition function

exp  

 1.3806503  10  23 m 2

H ( p , r ) 

kg s  2

K  1

  

p r

k B T

 drdp

 51

Illustration/example: phase space

6N-dimensional phase space

Image remo v e d d u e to co py ri g h t re stri cti o n s. S ee the se c o n d i m age at http ://www.ace . g a te ch.e d u /e x p e r ime n ts2/2413/lorenz/fall0 2 / .

r , p

 ( p , r )

Definition of temperature

Classical (mechanics) many-body system :

Average kinetic energy per degree of freedom is related to temperature via Boltzmann constant :

1 → 2

 1

  1

→ 2  1

2 m v i



N f i  1 .. N 

m v i

  k B T



# DOF

N f  3 N

# particles (each 3 DOF for velocities)

Based on equipartition theorem (energy distributed equally over all DOFs)

Definition of temperature

Classical (mechanics) many-body system :

Average kinetic energy per degree of freedom is related to temperature via Boltzmann constant :

→ 2 1  1 → 2  1

m v i

 

f i  1 .. N f 

m v i

  k B T



# DOF  p

Temperature

How to calculate temperature

1 1 N

2 → 2

 T 

   m i v i

 ( p , r ) drdp

p r 3

Nk B

i  1 m i

???

How to solve…

 A 

  A ( p , r )  ( p , r ) drdp

p r

Virtually impossible to carry out analytically

Must know all possible microscopic configurations corresponding to a macroscopic ensemble, then calculate 

Therefore: Require numerical simulation (the only feasible approach…)

Summary: How macro-micro relation works

A  A 

Statistical mechanics

 A 

microscopic A (single point measurement)

 A ( p , r )  ( p , r ) drdp

p r

defines…

3.3 How to solve the equations

Approaches in solving this problem

 Method of choice: Numerical simulation

 Two major approaches:

1. Using molecular dynamics (MD) : Generate microscopic information thr ough dynamical evolution of microscopic system ( i.e., simulate the “real behavior” as we would obtain in lab experiment)

2. Using a numerical scheme/algorithm to randomly generate microscopic states, which, through proper averaging, can be used to compute macroscopic properties. Methods referred to as “ Monte Carlo ”

Monte Carlo (MC) scheme

 Concept: Find simpler way to solve the integral

 A 

  A ( p , r )  ( p , r ) drdp

p r

 Use idea of “random walk” t o step through relevant microscopic states and thereby create proper weighting (visit states with higher probability density more often)

=ensemble (statistical) average

MC algorithm result

Final result of MC algorithm: Algorithm that leads to proper Distribution of microscopic states…

Ensemble (statistical) average

 A 

  A ( p , r )  ( p , r ) drdp

p r

Carry out algorithm for N A steps Average results

..done!

3.4 Ergodic hypothesis

Ergodicity

 MC method is based on directly computing the ensemble average

Define a series of microscopic states that reflect the appropriate ensemble average; weights intrinsically captured since states more likely are visited more frequently and vice versa

 Egodicity : The ensemble average is equa l to the time-average during the dynamic al evolution of a syst em under proper thermodynamical conditions.