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

How to model chemical interactions

Lecture 5

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 1-13

Lectures 14-26

2

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 – p roperty 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 method)

 Lecture 5: How to model chemical interactions 3

Lecture 5: How to model chemical interactions

Outline:

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

2. How to model chemical interactions

2.1 Pair potentials

2.2 How to model metals: Multi-body potentials

Goals of today’s lecture:

 Get to know basic methods to model chemical bonds (starting with simple “pair potentials”)

 Learn how to identify parameters for models of chemical bonds (for pair potentials)

 Limitations of pair potentials – and other, alternative methods

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

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

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 )

7

 A 

How to solve…

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

p r

Probability density distribution

E.g. :

Virtually impossible to carry out analytically Must know all possible configurations

Therefore: Require numerical simulation

Molecular dynamics OR Monte Carlo

8

Monte Carlo scheme

Method to carry out integration over “domain”



Want:

A  



f ( x ) d 

E.g. : Area of circle (=  exact solution)

A C 

 d 2

4





A C  4

 4 A C

d  1

f ( x )   1

inside outside

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

1 / 2

N A : Attempts

made

1 / 2 10

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

How to apply to ensemble average?

 Similar method can be used to apply to integrate the ensemble average

 A 

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

“discrete”

N A A exp   H ( r , p ) /  k T  

p r

1 

H ( p , r ) 

 A 



i  1

 exp  

A A

H ( r , p

B

) /  k

T  

 ( p , r ) 

exp  

Q 

k B T 

i  1

A A B

Computationally inefficient: If states are created “randomly” that have low probability….

 To be computationally more effective, need more complex iteration scheme (replace “ random sampling ” b y “ importance

sampling ”) 11

Challenge: sampling specific types of distributions

 We want to

 Integrate a sharply-peaked

f  x 

 exp    1 0 0 x  12 

function

 Us e Monte Carlo with uniformly-distributed random numbers (e.g. here

A    1

1

f ( x ) dx

from -1 to 1)

0.8

0.6

Random numbers drawn from here

0.4

0.2

- 1 -0. 5 0 .5 1

12

Challenge: sampling specific types of distributions

 We want to

 Integrate a sharply-peaked function

f  x   exp    1 0 0 x  12 

1

 Us e Monte Carlo with uniformly-distributed random numbers (e.g. here from -1 to 1)

A    1

1

f ( x ) dx

 What happens?

 Very few points contribute to the integral (~9%)

 Poor computational efficiency/convergence

Random numbers drawn from here

0.8

0.6

0.4

0.2

 Solution: use a different distribution of random numbers to sample “ importance sampling ”

- 1 -0. 5 0 .5 1

13

Importance sampling

 Core concept: Picking states with a biased probability: Importance sampling (sampling the “correct” way…)

 A 

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

 ( p , r ) 

1 exp  

H ( p , r ) 

  

p r

k B T 

Importance sampling

 Core concept: Picking states with a biased probability: Importance sampling (sampling the “correct” way…)

 A 

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

 ( p , r ) 

1 exp  

H ( p , r ) 

  

p r

k B T 

Notice: Probability (and thus importance) related to energy of state

Importance sampling: Metropolis algorithm

 Leads to an appropriate “chain” of states, visiting each state with correct probability

 Concept:

 Pick random initial state

 Move to trial states

 Accept trial state with certain probability (based on knowledge about behavior of system, i.e ., energy states)

Original reference: J. Chem. Phys. 21 ,1087, 1953

Metropolis-Hastings Algorithm

Concept: Generate set of random microscopic configurations Accept or reject with certain scheme

State A State B

Random move to 17

Metropolis-Hastings Algorithm: NVT

Have: State A (initial state) + en erg y functio n H(A)

Step 1: Gen e rate n e w state B (random move)

Metropolis-Hastings Algorithm: NVT

Have: State A (initial state) + en erg y functio n H(A)

Step 1: Gen e rate n e w state B (random move)

Step 2: if H(B)<H(A) then a = 1 a = true[1]/false[0]

“Downhill” m oves always accepted

else

Draw random number 0 < p < 1

for acceptance

Metropolis-Hastings Algorithm: NVT

Have: State A (initial state) + en erg y functio n H(A)

Step 1: Gen e rate n e w state B (random move)

Step 2: if H(B)<H(A) then a = 1 a = true[1]/false[0]

else

Draw random number 0 < p < 1

for acceptance

“Downhill” m oves

 H ( B )  H ( A ) 

always accepted,

if p  exp  

k T 

uphill moves with finite

else  B 

(“thermal”) probability

endif

a =variable either 0 or 1

Metropolis-Hastings Algorithm: NVT

Have: State A (initial state) + en erg y functio n H(A)

Step 1: Gen e rate n e w state B (random move)

Step 2: if H(B)<H(A) then a = 1 a = true[1]/false[0]

else

Draw random number 0 < p < 1

for acceptance

if p  exp  

H ( B )  H ( A ) 

else

 

 B 

endif Step 3: if

endif

endif

then accept state B

a =variable either 0 or 1 (used to detect acceptance of state B when a =1)

21

Metropolis-Hastings Algorithm: NVT

Have: State A (initial state) + en erg y functio n H(A)

Step 1: Gen e rate n e w state B (random move)

Step 2: if H(B)<H(A) then a = 1 a = true[1]/false[0]

else

Draw random number 0 < p < 1

for acceptance

repeat

N A times

 H ( B )  H ( A ) 

if p  exp  

k T 

else  B 

endif Step 3: if

endif

then accept state B

a =variable either 0 or 1 (used to detect acceptance of state B when a =1)

endif

 A 

1

N A i  1 .. N

A ( i )

A

Arrhenius law - explanation

E

Consider two states, A and B

H ( B )  H ( A )

H ( B )

H ( A )

H ( B )  H ( A )

A B iteration

State B has higher energy than state A

Otherwis e accepted anyway!

23

Arrhenius law - explanation

H ( B )  H ( A )

H ( B )  H ( A )

H ( B )

H ( A )

H ( B )  H ( A )

E E

iteration radius

Energy difference between states A and B

(“uphill”)

24

Arrhenius law - explanation

Probability of success

exp  

H ( B ) 

H ( A ) 



of overcoming the barrier

Random number 0 < p < 1

k B T 

(equal probability to draw any number between 0 and 1)

Acceptance if:

p  exp  

H ( B ) 

H ( A ) 

 k B T  1

E.g. when exp(..) = 0.8 most choices for p will be below, that is, high er

chance for acceptance

0.8

0.1

…. 0

low barrier

high barrier

H ( B )  H ( A )

Play “1D darts” 25

Summary: Metropolis-Hastings Algorithm

Have: State A (initial state) + en erg y functio n H(A)

Step 1: Gen e rate n e w state B (random move)

Step 2: if H(B)<H(A) then a = 1 a = true[1]/false[0]

else

Draw random number 0 < p < 1

for acceptance

repeat

N A times

 H ( B )  H ( A ) 

if p  exp  

k T 

else  B 

endif Step 3: if

endif

then accept state B

a =variable either 0 or 1 (used to detect acceptance of state B when a =1)

endif

 A 

1

N A i  1 .. N

A ( i )

A

Summary: MC scheme

Have achieved:

 A 

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

p r

Note:

• Do not need forces between atoms (for accelerations)

• O nly valid for equilibrium processes

27

Property calculation with MC: example

A

Iteratio n

“MC time” 28

Complex moves

 Move sets can be adapted for other cases, e.g. not just move of particles but also rotations of side chains (=rotamers), torsions , etc.

E.g. application in protein folding problem when we’d like to determine the 3D folded structure of a protein in thermal equilibrium

29

After: R.J. Sadus

Possible Monte Carlo moves

 Trial moves

 Rigid body trans lation

 Rigid body rotation

 Internal conformational changes (soft vs. stiff modes)

 Titration/electronic states

 …

 Questions:

 How “big” a move should we take?

 Move one particle or many?

Image by MIT OpenCou r seWare.

30

Monte Carlo moves

 How “big” a move should we take?

 Smaller moves : better acceptance rate, slower sampling

 Bigger moves : faster sampling, poorer acceptance rate

 Move one particle or many?

 Possible to achieve more efficient sampling with correct multi- particle moves

 One-particle moves must choose

particles at r andom Image by MIT OpenCou r seWare. 31

2. How to model chemical interactions

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)

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 34

Types of bonding (illustrations)

charge

+/ -

electron dens ity (localized!)

Ionic bonding

Hy drogen bonding

donor/ acceptor

Covalent bonding

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

+

© s ourc e u n known. All righ ts r e s e rved. This con t en t is excluded fr om our Creat i ve Co mmo ns li cense . Fo r m o re inf o rma t io n , see h ttp:/ /ocw.mi t . e du/fai ruse .

Metallic bonding

35

Wax

Cour t e sy of Ruth Ruan e, h t t p ://w ww.wh it ew it ch . i e . Used wi th permission.

Soft, deformable, does not break under deformation 36

Rocks

I m age courte sy o f Wi ki me di a Co mmo ns.

Quite brittle (breaks e.g. during earthquake) 37

Rocks and sand on Mars

What are the properties and composition of extraterrestrial rocks?

Ima g e cour t e sy of NASA.

Gold

I m age courte sy o f Wi ki me di a Co mmo ns.

Silicon

Ima g e cour t e sy of NASA.

Spider web

Imag e court e sy of U.S. Fi sh and W i l d li fe Servi c e.

41

Very extensible, deformation, yet very strong (similar to steel)

Tree’s leaf

I m age courte sy o f Wi ki me di a Co mmo ns. 42

Very deformable under bending (wind lo ads), but breaks easily under tear

BRITTLE

DUCTILE

Glass Polymers Ice...

Copper , Gold

Shear load

Particularly intriguing…brittle or ductile?

Imag e court e sy of quin n.anya . Li ce nse : CC - B Y .

Image by MIT OpenCou r seWare.

43

Outline

 Goal: model chemical bonds with the objective to enable force calculation (see lecture 2, basic MD algorithm) or energy calculation (see lecture 4/5, MC)

 Two-step approach :

1. Define energy landscape, i.e. defines how distance between particles controls the energy stored in the bond

2. Then take derivatives to obtain forces, to be used in the MD algorithm

“Modeling and simulation” p aradigm:

 First, develop mathematical expressions (modeling)

 Second, use model in numerical solution (simulation, =MD) 44

Models for atomic interactions

 Define interatomic potentials that de scribe the energy of a set of atoms as a function of their coordinates r :

U total

 U total

( r )

Depends on pos ition of all other atoms

r  r j 

j  1 .. N

Models for atomic interactions

 Define interatomic potentials that de scribe the energy of a set of atoms as a function of their coordinates r :

r  r j 

j  1 .. N

U total  U total ( r )

Depends on pos ition of all other atoms

2

i  1 .. N



 

  r 1 , i

, 

 r 2 , i

,  

 r 3 , i 

Change of potential energy due to change of position of particle i (“gradient”)

2.1 Pair potentials

Pair potentials: energy calculation

Simple approximation: Total energy is sum over the energy of all pairs of atoms in the system

r 12

r 25

r ij 

distance between particles i and j

Pair potentials: energy calculation

Simple approximation: Total energy is sum over the energy of all pairs of atoms in the system

r 12

r 25

r ij 

distance between particles i and j

 ( r ij )

Pair wise interaction potential energy for each bond

Pair potentials: energy calculation

Simple approximation: Total energy is sum over the energy of all pairs of atoms in the system

r 12

r 25

r ij 

distance between particles i and j

 ( r ij )

Pair wise interaction potential energy for each bond

En erg y of atom i

U i 

  ( r ij )

j  1 50

Overview - pair potentials: total energy calculation

Simple approximation: Total energy is sum over the energy of all pairs of atoms in the system

Pair wise interaction potential

 ( r ij )

r 12

r 25

Pair wise summation of bond energies

N

r ij 

distance between particles i and j

En erg y of atom i

U i 

  ( r ij )

j  1

51

Example: calculation of total energy

r 12

r 25

two “loops” o ver pairs of all particles

N N

U total 1

   ( r ij )

i  1 , i  j j  1

with

 ij

  ( r ij )

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 (no bond breaking)

How to use a pair potential, e.g. LJ

54

Force calculation – pair potential

Forces on particles can be calculated by taking derivatives from the potential function & by considering all pairs of atoms

Start with force magnitude (STEP 1): Negative derivative of potential energy with respect to atomic distance

F  

d  ( r ) d r

|

r  r ij

  d  ( r ij )

d r ij

   ' ( r ij )

j

f

r ij

f 2

x 2

f 1 x 1

i x 2

Force calculation – pair potential

Forces on particles can be calculated by taking derivatives from the potential function & by considering all pairs of atoms

Start with force magnitude (STEP 1): Negative derivative of potential energy with respect to atomic distance

F  

d  ( r ) d r

|

r  r ij

  d  ( r ij )

d r ij

   ' ( r ij )

Calculate force vector (STEP 2) :

Comp o  nent i of j

f

r ij

f 2

x 2

f 1 x 1

vector

r ij

f  F x i

f  

i f i e i

ij

i x 2

r  

ij r ij

What can we do with this potential?

Bending a copper wire until it breaks

A closer look

Court e sy of Goran Drazi c . Used wi t h permi ssion.

59

http://www2.ijs.si/~goran/sd96/e6sem1y.gif

Case study: plasticity in a micrometer crystal of copper

Simulation details

- 1 ,000,000,000 atoms ( 0.3 micrometer 

side length )

- 12:6 Lennard-Jones ductile material, for copper

[1 1 0]

(1 1 0)

Crack faces

[1 1 0]

Mode 1 tensile loading

- V isualization using energy filtering method ( only show high energy atoms )

Generic features of

atomic

r bonding:

„repulsion vs. attraction“