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 II

Lecture 4

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

Overview: Material covered so far…

 Lecture 1: Broad introduction to IM/S

 Lecture 2 : Introduction to atomistic and continuum modeling (multi-scale modeling paradigm, difference between continuum and atomistic approach, case study: diffusion)

 Lecture 3 : Basic statistical mechanics – property calculation I (property calculation: microscopic states vs. macroscopic properties, ensembles, probability density and partition function)

 Lecture 4 : Property calculation II (Advanced property calculation, introduction to chemical interactions, Monte Carlo methods)

Lecture 4: Property calculation II

Outline:

1. Advanced analysis methods: R adial distribution function (RDF)

2. Introduction: How to model chemical interactions

2.1 How to identify parameters in a Lennard-Jones potential

3. Monte-Carlo (MC) approach: Metropolis-Hastings algorithm

3.1 Application to integration

3.2 Metropolis-Hastings algorithm

Goal of today’s lecture:

- Learn how to analyze structure of a material based on atomistic simulation result (solid, liquid, gas, different crystal structure, etc.)

- Introduction to potential or force field ( Lennard-Jones )

- P resent details of MC algorithm – b ackground and

implementation 4

1. Advanced analysis methods: Radial distribution function (RDF)

Goals

 Define algorithms that enable us to “make sense” of positions, velocities etc. and time histories to relate with experimentally measurable quantities

 So far: temperature, MSD (mean square displacement function)

 Here: extend towards other properties

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.

Atomistic trajectory – through MSD

Court e sy of Sid Yip. Us ed with per m iss i o n .

Need pos itions over time – what if not available? 8

How to characterize material state (solid, liquid, gas)

 Application: Simulate phase transformation (melting)

© Trivedi Chemistry. All rights reserved. This content is excluded from our Creative Commons license. For more information, see http://ocw.mit.edu/fairuse .

http://ww w .t2i2edu.com/WebMov ie/1Chap1_files /image002.jpg

How to characterize material state (solid, liquid, gas)

Regular spacing

Neighboring particles found at characteristic distances

Irregular spacing

Neighboring particles found at approximate distances (smooth variation)

More irregular spacing

More random distances, less defined

© Trivedi Chemistry. All rights reserved. This content is excluded from our Creative Commons license. For more information, see http://ocw.mit.edu/fairuse .

10

How to characterize material state (solid, liquid, gas)

Regular spacing

Neighboring particles found at characteristic distances

Irregular spacing

Neighboring particles found at approximate distances (smooth variation)

More irregular spacing

More random distances, less defined

Concept:

© Trivedi Chemistry. All rights reserved. This content is excluded from our Creative Commons license. For more information, see http://ocw.mit.edu/fairuse .

 Measure distance of particles to their neighbors

 Average over large number of particles

11

 Average over time (MD) or iterations (MC)

Formal approach: Radial distribution function (RDF)

Ratio of density of atoms at distance r (in control area dr ) by overall density = relative density of atoms as function of radius

Reference atom

g ( r ) 

r

 ( r ) / 

dr

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

Formal approach: Radial distribution function (RDF)

The radial distribution function is defined as

Overall dens ity of atoms (volume)

g ( r )   ( r ) / 

Local dens ity

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 14

of radius r and thickness dr

Radial distribution function

 r considered volume

N ( r   r )

g ( r )  2

 ( r 

 r ) 

Density

  N / V

Note: RDF can be measured experimentally using x-ray or 15

neutron-scattering techniques

Radial distribution function:

Which one is solid / liquid?

© s ourc e u n known. All right s res e rved. This con t en t is exclud ed from our C r eat i ve Co mmo ns lice n se . Fo r mo re i n fo rm a t io n, see http:/ /ocw. m it.edu/fai ruse .

Interpretation: A peak indicates a particularly

favored separation distance for the neighbors to a given particle Thus, RDF reveals details about the atomic structure of the system being simulated

Java applet:

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

Radial distribution function

solid

liquid

© s ourc e u n known. All right s res e rved. This con t en t is exclud ed from our C r eat i ve Co mmo ns lice n se . Fo r mo re i n fo rm a t io n, see http:/ /ocw. m it.edu/fai ruse .

Interpretation: A peak indicates a particularly

favored separation distance for the neighbors to a given particle Thus, RDF reveals details about the atomic structure of the system being simulated

Java applet:

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

Radial distribution function: JAVA applet

Java a p plet: http://physchem.ox.ac.uk/~rkt/lec tures/liqsolns/liquids.htm l

Ima g e removed for c o pyright rea s on s.

Screen sh ot of th e radi al di stri but ion function Ja va appl et.

18

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 19

Notes: Radial distribution function (RDF)

 Pair correlation function (consider only pairs of atoms)

 Provides structural information

 Can provide information abou t dynamical change of structure, but not about transport properties (how fast atoms move)

Notes: Radial distribution function (RDF)

 Pair correlation function (consider only pairs of atoms)

 Provides structural information

 Can provide information abou t dynamical change of structure, but not about transport properties (how fast atoms move)

Additional comments:

 Describes how - on average - a toms in a system are radially packed around each other

 Particularly effective way of describing the structure of disordered molecular systems (liquids)

 In liquids there is continual movement of the atoms and a single snapshot of the system shows only the instantaneous disorder it is extremely useful to be able to deal with the average structure

Example RDFs for several materials

RDF and crystal structure

1 st nearest neighbor (NN)

3 rd NN 2 nd NN

…

Image by MIT OpenCou r seWare.

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

Face centered cubic (FCC), body centered cubic (BCC)

Image from Wi ki medi a Common s , h t t p : / / c om m o n s . w ik imed ia. o rg

Aluminum , NN: 2.863 Å ( a 0 =4.04 Å)

Copper , NN: 2.556 Å ( a 0 =3.615 Å) )

See also: http://www.webelements.com/

Chromium , N N :

2.498 Å ( a 0 =2.91 Å)

Iron , NN: 2.482 Å ( a 0 =2.86 Å) 24

Hexagonal closed packed (HCP)

a

a

a

c

Image by MIT OpenCou r seWare.

Image by MIT OpenCou r seWare.

Imag e court e sy of the U.S. Navy.

Cobalt

a : 250.71 pm

b : 250.71 pm

c : 406.95 pm

α : 90.000°

β : 90.000°

γ : 120.000°

NN: 2.506 Å

Zinc

a : 266.49 pm

b : 266.49 pm

c : 494.68 pm

α : 90.000°

β : 90.000°

γ : 120.000°

NN: 2.665 Å 25

Graphene/carbon nanotubes

RDF

Imag es removed du e t o copyrigh t restric t ion s. Pl e a se see : h t t p : //w eblog s 3. n r c. n l / t ech n o /w p- c o n t ent/up load s/08042 4_Grafeen/Graphen e_xyz.jpg ht tp:/ /dept s . w ashi ngt o n.edu/pol ylab/i m ages/cn1. j pg

Graphene/carbon nanotubes (rolled up graphene)

NN: 1.42 Å, second NN 2.46 Å … 26

Macroscale view of water

Iceberg

Glacier

Imag e court e sy of dnkemon t oh .

27

Imag e court e sy of blmiers 2 .

RDF of water (H 2 O)

Court e sy of Mar k Tuckerma n. Used wi th permission. 28

http://ww w .nyu.edu/class es/tuckerman/stat.mech/l ectures/lecture_8/node1.html

RDF of water (H 2 O)

O-O distance

Imag es c o u r t e sy of Ma r k Tuckerman. Used wit h permis s i on.

H-O distance

http://ww w .nyu.edu/class es/tuckerman/stat.mech/l ectures/lecture_8/node1.html

RDF of water (H 2 O)

O-O distance

H- bonding ( ≈ 2.8 Å)

Imag es c o u r t e sy of Ma r k Tuckerman. Used wit h permis s i on.

H-O distance

H-O covalent bonding ( ≈ 1 Å)

http://ww w .nyu.edu/class es/tuckerman/stat.mech/l ectures/lecture_8/node1.html

2. Introduction: How to model chemical

interactions

31

Molecular dynamics: A “bold” idea

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

Forces between atoms… how to obtain? 32

How are forces calculated?

Force magnitude: Derivative of potential energy with respect to atomic distance

f   d U ( r )

d r

To obtain force vector f i , t ake projections into the three axial directions

Atomic interactions – quantum perspective

How electrons from different atoms interact defines nature of chemical bond

Density distribution of electrons around a H-H molec u le

Ima g e rem o ved du e t o copyrigh t re str i c t ion s . Pleas e s e e t h e an ima t io n of hydrogen bon ding orb i t a ls at h t t p : //w in t e r. grou p. sh e f . ac. uk/ orb itron/MOs/H 2 /1s1s-sig m a/ind e x.html

Much more about it in part II

Concept: Interatomic potential

Ener gy U

r

1/r 12 (or Exponential)

Repulsion

Radius r (Distance between atoms)

e

Attraction

1/r 6

“point particle” representation

Image by MIT OpenCou r seWare.

Interatomic bond - model

Ener gy U

r

1/r 12 (or Exponential)

Repulsion

Radius r (Distance between atoms)

e

Attraction

1/r 6



r 0

Harmonic oscillator

~ k(r - r 0 ) 2

r

Image by MIT OpenCou r seWare. Image by MIT OpenCou r seWare.

Atomic interactions – different types of chemical bonds

 Prima ry bonds (“strong”)

 Ionic (ceramics, quartz, felds par - rocks )

 Covalent ( silicon )

 Metallic (copper, nickel, gold , silver) (high melting point, 1000-5,000K)

 Secondary bonds (“weak”)

 Van der Waals ( wax , low melting point)

 Hy drogen bonds (proteins, spider silk ) (melting point 100-500K)

 Ionic: Non-directional (point charges interacting)

 Covalent: Directional (bon d angles, torsions matter)

 Metallic: Non-directional (electron gas concept)

Difference of material properties originates from different atomic interactions 37

Interatomic pair potentials: examples

 ( r ij ) 

D exp   2  ( r

 r 0

)   2 D exp    ( r

 r 0 ) 

Morse potential

 ( r ij )

 4    

 

 12   6 

   

  ij   

Lennard-Jones 12:6 potential

(excellent model for noble Gases, Ar, Ne, Xe..)

 r 

   6

 ( r )  A exp   i j   C  

ij   

 r 

Buck ingham potential

  

ij 

Harmonic approximation

What is the difference between these models?

Shape of potential (e.g. behavior at short or long distances, around equilibrium)

Number of parameters (to fit) Ability to describe bond breaking

Lennard-Jones potential

Parameters , 

  

 12

   6 

 ( r ) 

     

 

Sir J. E. Lennard-Jones (Cambridge U K )

Lennard-Jones 12:6

Lennard-Jones potential: schematic & parameter meaning

  : well depth (energy stored per bond)

  proportional to point where force vanishes (equilibrium distance between atoms)

r

  

 12

   6 

Lennard-Jones 12:6

 ( r ) 

     

 

2.1 How to identify parameters in a Lennard-Jones potential

(=force field training, force field fitting, parameter coupling, etc.)

Parameter identification for potentials

 Typically done based on more accurate (e.g. quantum mechanical ) results (or experimental measurements, if available)

 Properties used include:

Lattice constant, cohesive bond energy, elastic modulus (bulk, shear, …), equations of state, phonon frequencies (bond vibrations), forces, stability/energy of different crystal structures, surface energy, RDF, etc.

Potential should closely re produce these reference values

 Challenges : mixed systems, different types of bonds, reactions

Multi-scale paradigm

 Show earlier: molecular dynamic s provides a powerful approach to relate the diffusion constant that appears in conti nuum models to atomistic t rajectories

 Force field fitting to identify param eters for potentials (based on quantum mechanic al results) is yet another “step” in this multi-scale paradigm

Time scale

force fie l d fitting

diffusion calculation s

MD

Continuum model

“Empirical”

or experimental parameter feeding

Quantu

44

Derivative of LJ potential ~ force

0.2

0.1

F = _ d  ( r ) = -  '

d r

0

-0.1

 ( r )

F max, LJ

-0.2

2

r 0

3

4

5

r (A)

o

Image by MIT OpenCou r seWare.

EQ r 0

relates to equilibrium spacing crystal 45

Properties of LJ potential as function of

parameters , 

Equilibrium distance between atoms r 0 and maximum force



 r 0

F max, LJ

 2 . 394 



first derivative zero (force)

second derivative

zero (=loss of convexity, spring constant=0)

FCC

r 0 

distance of nearest ne ighbors in a lattice r 0

Image from Wi ki medi a Common s , ht tp:/ /commons.wi ki m edi a . o rg

Determination of parameters for atomistic interactions

 Example (based on elastic properties) of FCC lattice

 Approach: Express bulk modulus as function of potential parameters

 Second derivative of potential is related to spring constant

(=stiffness) of chemical bonds

  1 / 4

Young’s modulus

Shear modulus

K  E /( 3 ( 1  2  ))

E  8 / 3 

  r 2 k

/ 2 / V

V  3 / 4

47

Determination of parameters for atomistic interactions

 Example (based on elastic properties) of FCC lattice

 Approach: Express bulk modulus as function of potential parameters

 Second derivative of potential is related to spring constant

(=stiffness) of chemical bonds

  1 / 4

Young’s modulus

Shear modulus

K  E /( 3 ( 1  2  ))

E  8 / 3 

  r 2 k

/ 2 / V

V  3 / 4

 2  ( r )

k   r 2

  ' '

Determination of parameters for atomistic interactions

 Example (based on elastic properties) of FCC lattice

 Approach: Express bulk modulus as function of potential parameters

 Second derivative of potential is related to spring constant

(=stiffness) of chemical bonds

  1 / 4

Young’s modulus

Shear modulus

K  E /( 3 ( 1  2  ))

E  8 / 3 

  r 2 k

/ 2 / V

V  3 / 4

k 

 r 2

  ' '

Bulk modulus copper E = 140 GPa

Lennard-Jones potential – example for copper

0.2

0.1

0

-0.1

-0.2

2

r 0 3 4 5

r (A)

Image by MIT OpenCou r seWare.

LJ potential – parameters for copper 50

3. Monte Carlo approaches

How to solve…

 A 

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

p r

Probability density distribution

Virtually impossible to carry out analytically Must know all possible configurations

Therefore: Require numerical simulation

Molecular dynamics OR Monte Carlo

3.1 Application to integration

“Random sampling”

Monte Carlo scheme

Method to carry out integration over “domain”



Want:

A  



f ( x ) d 

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

 d 2 

A C  4 A C  4

  4 A C

d  1

f ( x )   1

inside outside

Conventional way…

 Evaluate integrand at predetermi ned values in the domain (e.g. quadratic grid)

 Evaluate integral at discrete points and sum up

…..



What about playing darts..

Public domain imag e.

Alternative way: integration through MC

…..

Playing darts: Randomly select point in domain Evaluate integral a these points

Sum up results to solve integral

Monte Carlo scheme for integration

 Step 1 : Pick random point 

in 

 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 to the total sum

f ( x i )  1

A   

A  

C  f



( x ) d

C 16

N A : Attempts

made

Court e sy of John H. Ma th ews. Used wi th permi s sion.

58

http://math.fullerton.edu/mat hews/n2003/MonteCarloPiMod.html

Java applet: how to calculate pi

 http://polymer.bu.edu/java/jav a/montepi/montepiapplet.html

Court e sy of the Cent er for Pol ymer Studies a t Bo sto n Uni v e r sity. U s e d wi th per m is sio n .