1.021 , 3.021, 10.333, 22.00 I ntroduc tion to Modeling and Simulation Part I – C ontinuum and partic le me thods

Review session—preparation quiz I

Lecture 12

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

2

Overview: Material covered so far…

 Lecture 1: B road introduction to IM/S

 Lecture 2 : Introduction t o atomistic and conti nuum modeling (mult i-scale modeling paradigm, difference between continuum and atomis tic approach, case study: diffusion)

 Lecture 3 : Basic statistical mechani cs – p roperty calculation I (property calculati o n: microscopic states vs. macroscopic proper ties, ensembl es, probabi lity density and partiti o n function)

 Lecture 4 : Prope rty calculation II (Monte Carl o, advanced property ca lcu l ati o n, introd uctio n to chemical interacti ons)

 Lecture 5: How to model chemic al intera ctions I (exampl e: movi e of copper deformation/disl ocations, etc.)

 Lecture 6: How to model ch emic al intera ctions II (EAM, a bit of ReaxFF—chemical reacti ons)

 Lecture 7: Appli cation to mo del i ng brittle materials I

 Lecture 8: Appli cation to mo del i ng brittle materials II

 Lecture 9: Appli cation – A pplications to materials failure

 Lecture 10: Appl ications to bi ophysics an d bi onanomechanics

 Lecture 11: Appl ications to bi oph ysics an d bi onanomechanics (cont’ d)

3

 Lecture 12: Review session – p art I

Check: goals of part I

You will be able to …

 Carry out atomistic simulations of various processes (diffusion, deformation/stre tching, materials failure)

Carbon nanotubes, nanowires, bu lk metals, proteins, silicon crystals, …

 Analyze atomistic simulations

 Visualize atomistic/molecular data

 Understand how to link atomistic simulation results with continuum models within a multi-scale scheme

Topics covered – part I

Lecture 12: Review session – p art I

1. Review of material covered in part I

1.1 Atomistic and molecular simulation algorithms

1.2 Property calculation

1.3 Potential/force field models

1.4 Applications

2. Important terminology and concepts

Goal of today’s lecture:

 Review main concepts of atomistic and molecular dynamics, continuum models

 Prepare you for the quiz on Thursday

1.1 Atomistic and molecular simulation algorithms

1.2 Property calculation

1.3 Potential/force field models

1.4 Applications

Goals: Basic MD algorithm (integrat ion scheme), initial/boundary conditions, numeric a l i ssues (supercomputing)

Differential equations solved in MD

Basic concept: atomistic methods

PDE solv ed in MD:

2 

m d r i 

dU ( r )





r { r }

i  1 .. N

( ma  F )

i dt 2

d r i

Results o f M D simulation:

  

r i ( t ), v i ( t ), a i ( t )

i  1 .. N 9

Complete MD updating scheme

(1) Initial conditions: Positions & velocities at t 0 (random velocities so that initial temperature is represented)

(2) Updating method (integration scheme - V erlet)

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

( t 0

)   t  2

 ...

Positions at t 0 -  t

(3) Obtain accelerations from forces

Positions at t 0

Accelerations at t 0

f j , i

 m j a j , i

a j , i 

f j , i

/ m j

 j  1 .. N

(4) Obtain force vectors from potential (s um over contributions from all neighbor of atom j )

d  ( r )

x j , i

  

 12

   6 

F  

d r

f j , i  F r

Potential

 ( r ) 

    







 10

Algorithm of force calculation

6

…

5

j =1

i

4

2

3

r cut

for i=1..N

for j=1..N ( i ≠ j )

[add force contributions]

Loop

Summary: Atomistic simulation – numerical approach “molecular dynamics – M D”

 Atomis tic model; requires atomistic mi crostructure and atomic positions at beginning (initial conditions), initia l velocities from random distribution according to specified temperature

 Step through time by integration scheme (Verlet)

 Repeat force calculation of atom ic forces based on their positions

 Explic it notion of chemical bonds – c aptured in interatomic potential

Scaling behavior of MD code

Force calculation for N particles: pseudocode

 Requires two nested loops, first over all atoms, and then over all other atoms to determine the distance

for i=1..N :

for j=1..N ( i ≠ j ):

 t 1 ~N  t 0

determine distance between i and j

calculate force and energy (if

r ij

< r cut , cutoff radius)

 t 0

add to total force vector / energy

time ~ N 2 : computational disaster 14

Strategies for more efficient computation

Two appr oaches

1. Neighbor lis ts: Store information about atoms in vicinity, calculated in an N 2 effort, and keep information for 10..20 steps

Concept: Store information of neighbors of each atom within vicinity of cutoff radius (e.g. list in a vector); update list only every 10..20 steps

2. Domain decomposition into bins: Decompose system into small bins; force calc ulatio n only between atoms in local neighboring bins

Concept: Even overall system grows, calculation is done only in a loc a l environment (have two nest ed loops but # of atoms does not

increase locally) 15

Force calculation with neighbor lists

 Requires two nested loops, first over all atoms, and then over all other atoms to determine the distance

for

 t 1 ~N  t 0

time ~ N : computationally tractable even for large N 16

Periodic boundary conditions

 Periodic boundary conditions allows studying bulk properties (no free surfaces) with small number of particles (here: N =3 ), all particles are “connected”

 Original cell surrounded by 26 image cells; image particles move in exactly the same way as original particles (8 in 2D)

Image by MIT OpenCou r s e Ware. After Buehl e r. 17

Parallel computing – “ supercomputers”

Imag es removed du e t o copyri g h t r e s t ri ctio ns . Pl e a se see h t t p : //www. san d ia. g ov / ASC / i ma ge s/ library / ASC I- Wh it e. jp g

Supercomputers consis t of a very large number o f individual computing units (e.g. Central Processing

Units, CPUs) 18

Domain decomposition

Each piece worked on by one of the computers in the supercomputer

Fig. 1 c from Buehl e r, M., et a l . "The D y nami cal Complexity of Wor k -Harden i ng: A Large - Scal e Mol e c u l a r Dynami c s Si mul a ti on.”

Acta Me c h S in i ca 2 1 (20 05) : 103-1 1 . © S p r ing e r- Ve rlag . A l l rig h ts rese rved . T h is c o ntent i s e x c l ud ed from our Cre a t i ve Commons lic ense .

For more i n fo rmati on, see h t t p :// oc w . m i t . e d u/ f a ir use . 19

Modeling and simulation

Modeling and simulation

What is a model? What is a simulation?

Modeling vs Simulation

• Mod e lin g : developing a mathem atical representation of a physi cal sit u ati o n

• Si mulation : solving the equations that arose from the development of the model.

© G oo gl e, Inc. Al l ri ghts r e se r v e d . This content is ex cl u d ed from ou r Creative Common s licen se. For more in fo r m a t io n , se e h t t p :// oc w . m i t . e d u/ f a ir use .

© M a s s achuse t t s B a y Tran s p o rtati o n Authori t y. All ri ghts re s e rv ed. Thi s c o nte n t i s ex cl u d ed from our Creative Common s lic ense . F o r mo re in fo r m a t io n , s ee h t t p :// oc w . m i t . e d u/ f a ir use .

Modeling and simulation

M S M S M

Atomistic versus continuum viewpoint

Diffusion: Phenomenological description

Physical problem

Ink dot

© s o ur c e u n k n o w n. Al l rights res e rv ed. Thi s co ntent is ex cl u d ed from ou r Crea t i ve Common s licen se. For mo re in fo r m a t io n , s e e htt p :/ / o cw.mi t .edu /fai ru se

Tissue

BC: c ( r = ∞ ) = 0 24

IC: c ( r=0 , t = 0) = c 0

Continuum model: Empirical parameters

 Continuum model requires paramet er that describes microscopic processes inside the material

 Need experimental measurements to calibrate

Time scale

“Empirical”

or experimental parameter feeding

Length scale 25

How to solve continuum problem: Finite difference scheme

 c 

 t

d 2 c D dx 2

Concentration at i

at old time

i  space

j  time

c  c

  t D  c

 2 c  c

 “explicit” numerical scheme

i , j  1

i , j

 x 2

i  1 , j

i , j

i  1 , j

(new concentration directly from

concentration at earlier time)

Concentration at i

at old time

Concentration at i

at new time

Concentration at i+ 1

at old time

Concentration at i- 1

at old time

26

Other methods : Finite element method

Atomistic model of diffusion

 Diffusion at macroscale (change of concentrations) is result of microscopic processes, random motion of particle

 Atomis tic model provides an alternative way to describe diffusion

 Enables us to directly calculat e the diffusion constant from the trajectory of atoms

 Follow trajectory of atoms and calc ulate how fast atoms leave their initial position

Concept:

follow this quantity over time

27

Atomistic model of diffusion

2 

m d r i 

dU ( r )

 r

{  }

i  1 .. N

i dt 2 d r

i

Images cou rtesy of the Center for Polymer Studi e s at Bo sto n Uni v e r si ty. Use d wi th pe rmi s si o n .

 r 2

Particles Trajectories

Mean Squared Displacement function

Average squares of dis p lacements of all particles

 r 2 ( t ) 

1   

   0 )  2

r i ( t )

i

r i ( t

28

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

(nm)

 r 2

slope = D

Einstein equation

(ps)

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

29

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

30

1.1 Atomistic and molecular simulation algorithms

1.2 Property calculation

1.3 Potential/force field models

1.4 Applications

Goals: How to calc ulate “useful” properties from MD runs (temperature, pressure, RDF, VAF,..); significance of averaging; Monte Carlo schemes

Property calculation: Introduction

Have:

  

r i ( t ), v i ( t ), a i ( t )

i  1 .. N

“microscopic information” Want:

 Thermodynamical properties (temperature, pressure, ..)

 Transport properties (diffusivities, shear viscosity, ..)

 Material’s state (gas, liquid, solid)

Micro-macro relation

T ( t )

Which to pick?

1 1 N  t

T ( t ) 

 m v 2 ( t )

Images cou rtesy of the Center for Polymer Studies at Bosto n Uni v e r si ty. Use d wi th permi s si o n.

3 Nk B

i i

i  1

Ergodic hypothesis: significance of averaging

 A 

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

p r

property must be properly averaged

 Ergodic hypothesis:

Ensemble (statistical) average = time average

 All microstates are sampled with appropriate probability dens ity o ver long time scales

1 

A ( i ) 

A  Ens 

A  Time  N

 A ( i )

Monte Carlo

Average over Monte Carlo steps

NO DYNAMICAL INFORMATION

Monte Carlo scheme

 Concept: Find convenient 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)

 Monte Carlo schemes : Many applications (beyond what is discussed here; general method to solve complex integrations)

Monte Carlo scheme: area calculation

Method to carry out integration (illustrate general concept)

Want:

A  



f ( x ) d 

E.g. : Area of circle (value of  )

 d 2  

A C  4

A C  4

d  1

Monte Carlo scheme

 Step 1: Pick random point 

 Step 2: Accept/reject point based on criterion (e.g. if inside or outside of circle and if in area not yet counted)



 Step 3: If accepted, add

f ( x i )   1

to the total sum otherwise

f ( x i )  0

d  1 / 2

A C  



f ( x ) d 

N A : Attempts

made

Courtes y of Jo hn H. Mathew s. Us ed wi t h pe rmi s si o n . 37

See also: http://math.fullerton.edu/mathews/n2003/MonteCarloPiMod.html

Area of Middlesex County (MSC)

A

SQ

 10 , 000 km 2

x

y

100 km

Image courtesy of Wikimedia Commons.

Fraction of points that lie within MS County

A MSC

1

A SQ N



f ( x i )

Expression provides area of MS County

A i 38

Detailed view - schematic

N A =55 points (attempts) 12

12 points within

A MSC

1

A SQ

A

i

 10000  0 . 22 km 2

 2181 . 8 km 2

39

Results U.S. Census Bureau

823.46 square miles (1 square mile=2.58998811 km 2 )

Taken from: http://quickfacts.census.g ov/qfd/states/25/25017.html

Monte Carlo result:

40

Analysis of satellite images

Image remo v e d du e to co py ri ght restri cti o n s .

G oog le Sa t e ll it e ima g e o f Sp y P o n d , in Ca mb r i d ge , M A .

200 m 41

Metropolis Hastings algorithm

Images r e mo ve d d u e to c o py ri ght restri cti o n s .

Pl ease s e e: Fi g. 2. 7 i n Bue h l e r, Mark u s J. At o m i s t i c Mo de li n g of Mater i al s Fail ur e . S p ri n g er , 20 0 8.

Images r e mo ve d d u e to c o py ri ght restri cti o n s .

Pl ease s e e: Fi g. 2. 8 i n Bue h l e r, Mark u s J. At o m i s t i c M o d e l i n g of Material s Fail ur e . Sp ringe r , 2 008 .

 A 

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

p r 42

Molecular Dynamics vs. Monte Carlo

 MD provides actual dynamical data for nonequilibrium processes (fracture, deformation, instabilities)

Can study onset of failure, instabilities

 MC provides information about equilibrium properties

(diffusivities, temperature, pressure) Not suitable for processes like fracture

1

N A i  1 .. N

A ( i ) 

A

A  Ens 

A  Time 

1

N t i  1 .. N

A ( i )

43

t

Property calculation: Temperature and pressure

 Temperature

1 1 N  2

  

T 

3 Nk B

  m i v i 

i  1

i v i

 v i

44

Temperature control in MD ( NVT ensemble) Berendsen thermostat

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

)  t 2

 ...

i 0 i 0 i 0 i 0

1 N  2

Calculate temperature

T 

3 Nk B

  m i v i 

i  1

v ~

Calculate ratio of current vs. desired temperature

Rescale velocities

Interpretation of RDF

46

Formal approach: Radial distribution function (RDF)

The radial distribution function is defined as

Provides information about the density of atoms at a given radius r ;  ( r ) is the local density of atoms

g ( r ) 2  r 2 dr 

Number of particles that lie in a spherical shell 47

of radius r and thickness dr

Radial distribution function: Solid versus liquid versus gas

Image by MIT OpenCou r seWare.

Note: The first peak corresponds to the nearest neighbor shell, the second peak to the second nearest neighbor shell, etc.

In FCC: 12, 6, 24, and 12 in first four shells 48

RDF and crystal structure

1 st nearest neighbor (NN)

5

4

3

2

1

0 0 0 . 2 0 . 4 0 . 6 0 . 8

d i s t a n c e / n m

Image by MIT OpenCou r seWare.

Peaks in RDF characterize NN distance, can infer from RDF about crystal structure

49

http://physchem.ox.ac.uk/~rkt/lectures/liqsolns/liquids.html

Imag es c o u r t e sy of Ma r k Tucker man. Used wi th permi ssion.

Interpretation of RDF

1 st nearest neighbor (NN) 2 nd NN

3 rd NN

copper: NN distance > 2 Å

X

not a liquid (sharp peaks) 51

Graphene/carbon nanotubes

RDF

Images of gr aph e ne/ c ar bo n nan utu bes :

http: / / webl o g s 3 .nr c . nl / tec h n o /w p- conte n t/ up loads / 0804 24_G rafe en/ G raphene _xyz.jpg

http: / /d e pts. was h i n gto n . e d u / p ol yl ab /i mage s / c n 1. j p g

Graphene/carbon nanotubes (rolled up graphene)

NN: 1.42 Å, second NN 2.46 Å … 52

Summary – property calculation

Material properties: Classification

 Structural – c rystal structure, RDF

 Thermodynamic -- equation of stat e, heat capac ities, thermal expansion, free energies, use RDF, temperature, pressure/stress

 Mechanic al -- elastic constants, cohes ive and shear strength, elastic and plastic deformation, fracture toughness, use stress

 Trans port -- diffusion, viscous flow, thermal conduction, use MS D, temperature

Bell model analysis (protein rupture)

Bell model analysis (protein rupture)

a =86.8 6 p N

f ( v ; x b , E b )  a  ln v  b

b= 800 pN

a  k B T

x b

(1)

k B =1.3806505E-23 J/K

Get a and b from fitting to the graphs

k T 

 E  

b   b l n   x

x b 

ex p   b  

 k b T  

(1)

So lve (1) for a , use (2) to solve for E b

  1  1 0 13 1 / sec

a= (800-400)/(0-(-4.6))) pN

=86.86 pN  x b =0.47 Å 56

Determining the energy barrier

b  

k b T



ln   0 x b



exp   b

  

 k b T

 ln   x  

E b 



 bx b

x b  

 ln   x   E b

k b T 

x b 

k b T 

0 b

k b T k T

E b 

ln   x

  bx b

0 b

k b T k T

E   ln   x   bx b  k T

b  0 b



 b

k b T 

b= 800 pN  E b ≈ 9 k cal/mol

1 kcal/mol=6.9477E-21 J

Poss ible mechanis m: Rupture of 3 H-bonds (H-bond energy 3 kcal/mol)

Scaling with E b : shifts curve

f ( v ; x b , E b )  a  ln v  b

E b 

k B  T

k B  T

 E b 

a  b  

x b x b

 ln v 0

v 0   0

 x b

 exp  

58 



k b  T 

Scaling with x b : changes slope

f ( v ; x b , E b )  a  ln v  b