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

Reactive potentials and applications

Lecture 8

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 st atistical mech anics – p ropert y calculat ion I (property calculation: microscopic stat es vs . macroscopic properties, ensembles, probability density and partition function)

 Lecture 4 : Pr operty calculation II (Monte Carlo, advanced property calculation, introduction to chemical interactions)

 Lecture 5: How to model ch emical interactions I (example: movie of copper deformation/dislocations, etc.)

 Lecture 6: How to model ch emical interactions II (EAM, a bit of ReaxFF— chemical r e actions)

 Lecture 7: Applicat ion – M D si mulation of materials failure

 Lecture 8: Application – R eactiv e potentials and applications

Lecture 8: Reactive potentials and applications

Outline:

1. Bond order force fields - h ow to model chemical reactions

1.1 EAM potential for metals

1.2 ReaxFF force field

2. Hybrid multi-paradigm fracture models

Goal of today’s lecture:

 Learn new potential: ReaxFF, to describe complex chemistry (bond breaking and formation)

 Application in hybrid simulation approaches (combine different force fields)

1. Bond order force fields - h ow to model chemical reactions

Potential energy expressions for more complex materials/chemistry, including bond formation and breaking

Review: 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

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)

Another example: harmonic and harmonic bond snapping potential (see lecture 7)



  1 k ( r  r ) 2

2 0 0



 1

 2

   1

0 ( r  r 0 )

r  r break



 2 0

( r break

 r 0 )

r  r break

Example: calculation of total energy “simply” the sum of all energies of pairs of atoms

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 )

U total  2

12   13

  14

  1 N

...   21

  23

 ...   2 N

 ...  

N  1 , N 

But…are all bonds the same? - valency in hydrocarbons

Ethane C 2 H 6

(stable configuration)

H

All bonds are not the same! Adding another H is not favored

Bonds depend on the environment! 10

Are all bonds the same? – metallic systems

Surface

Bulk

stronger

+ different bond EQ distance

Pair potentials: All bonds are equal!

Reality: Have environment effects; it matter that there is a free sur f ace!

Bonds depend on the environment!

11

Are all bonds the same?

Bonding energy of red atom in is six times bonding energy in

This is in contradiction with bot h experiments and more accurate quantum mechanic al calc ulations on many materials

Bonding energy of atom i

U i 

  ( r ij )

j  1

Are all bonds the same?

Bonding energy of red atom in is six times bonding energy in

This is in contradiction with bot h experiments and more accurate quantum mechanic al calc ulations on many materials

For pair potentials ~ Z

For m e tals ~

Z : Coordination = how many immediate neighbors an atom has

E.g. in metals, bonds get “weaker” as more atoms are added to central atom

Bond strength depends on coordination

energy per bond

potential

~

Nickel

Z

2 4 6 8 10 12 coordination

Transferability of pair potentials

 Pair potentials have limited transferability:

Parameters determined for molecules can not be used for crystals, parameters for specific types of crystals can not be used to describe range of crystal structures

 E.g. difference between FCC and BCC can not be captured using a pair potential

1.1 EAM potential for metals

Note: already used in pset #1, nanowire simulation

Metallic bonding: multi-body effects

 Need to consider more details of chemical bonding to understand environmental effects

Electron (q=-1) Ion core (q=+N)

Delocalized valence electrons moving between nuclei generate a binding force to hold the atoms together: Electron gas model ( positive ions in a sea of electrons )

Mostly non-directional bonding, but the bond strength indeed depends on the environment of an atom, precisely the electron density imposed by other atoms

Concept: include electron density effects

  , j ( r ij )

Each atom features a particular distribution of electron density

Concept: include electron density effects

Electron density at atom i

 i 

   , j ( r ij )

j  1 .. N neigh

  , j ( r ij )

Atomic electron density of atom j

r ij

Contribution to electron density at site i due to electron density of atom j evaluated at correct distance ( r ij )

Concept: include electron density effects

 i   2  ( r ij ) 

F (  i )

j  1 .. N neigh

Electron density at atom i  i 

  , j ( r ij )

Embedding term F (how local electron density contributes to potential energy)

   , j ( r ij )

j  1 .. N neigh

Atomic electron

j i density of atom j

Embedded-atom method (EAM)

Atomic energy

1

new

Total energy

N

 i  

j  1 .. N neigh

2  ( r ij ) 

F (  i )

U total

   i

i  1

Pair potential energy Embedding ener gy

as a function of electron density

 i Electron density at atom i

based on a “pair potential”:

 i 

   , j ( r ij )

j  1 .. N neigh

First proposed by Finnis, Sinclair, Daw, Baskes et al. (1980s)

Physical concept: EAM potential

 Describes bonding energy due to electron delocalization

As electrons get more states to spread out over their kinetic energy decreases

 When an impurity is put into a metal its energy is lowered because the electrons from the impurity can delocalize into the solid.

 The embedding density (electron density at the embedding site) is a measure of the num ber of states available to delocalize onto.

 Inherently MANY BODY effect!

Effective pair interactions

r

+ + + + + +

+ + + + + +

+ + + + + +

r

+ + + + + +

+ + + + +

+ + + + +

+

+

Image by MIT OpenCou r seWare.

Can describe differences between bulk and surface

Summary: EAM method

 State of the art approach to model metals

 Very good potentials available for Ni, Cu, Al since late 1990s, 2000s

 Numerically efficient, can treat billions of particles

 Not much more expensive than pair potential (approximately three times), but describes physics much better

 Strongly recommended for use!

Another challenge: chemical reactions

sp3 sp2

Energy

 stretch

C-C distance

r

Simple pair potentials can not describe chemical reactions

Why can not model chemical reactions with spring- like potentials?

 stretch

 1 k

2

stretch

( r  r 0 )

Set of parameters only valid for particular molec ule type / type of chemical bond

stretch, sp 2 stretch, sp 3

Reactive potentials or reactive force fields overcome these limitations

1.2 ReaxFF force field

For chemical reactions, catalysis, etc.

Key features of reactive potentials

 How can one accurately describe the transition energies during chemic al reactions?

 Us e computationally more efficient descriptions than relying on purely quantum mechanical (QM) methods (see part II, methods limited to 100 atoms )

H  C

 C 

 H 2

H  H 2  C

 C  H 2

inv o lves processes with electrons

A

A- -B

B

Key features of reactive potentials

 Molecular model that is capable o f describing chemical reactions

 Continuous energy landscape during reactions ( k ey to enable integration of equations)

 No typing necessary, that is, atoms can be sp, sp2, sp3… w/o further “tags” – only element types

 C o m p utation ally efficient (that is, should involve finite range interactions), so that large sys tems can be treated (> 10,000 atoms)

 Parameters with physical meaning ( such as for the LJ potential )

Theoretical basis: bond order potential

Concept: Use pair potential that depen ds on atomic environment (similar to EAM, here app lied to covalent bonds)

Modulate strength of attractive part

(e.g. by coordination, or “bond order”)

Abell, Tersoff

Image by MIT OpenCou r seWare.

Changes in spring constant as function of bond order Continuous change possible

= continuous energy landscape during chemical reactions

Theoretical basis: bond order potential

Image by MIT OpenCou r seWare.

D. Brenner, 2000 31

Concept of bond order (BO)

r

BO

sp3 1

sp2 2

sp

3

Bond order based energy landscape

Bond length Bond length

Pauling

Bond or der

Energy Energy

Bond order potential Allows for a more general description of chemistry

All energy terms dependent on bond order

Conventional potential (e.g. LJ, Morse)

Historical perspective of reactive bond order potentials

 1985: Abell: General expression for binding energy as a sum of near nieghbor pair interactions moderated by local atomic environment

 1990s: Tersoff, Brenner: Use Abell formalism applied to silicon (successful for various solid state structures)

 2000: Stuart et al.: Reactive potential for hydrocarbons

 2001: Duin, Godddard et al.: Reactive potential for hydrocarbons “Reax FF”

 2002: Brenner et al.: Second generation “REBO” potential for hydrocarbons

 2003-2005: Extension of ReaxFF to various materials inc l uding metals, ceramics, silicon, polymers and more in Goddard‘s group

Example: ReaxFF reactive force field

William A. Goddard III

California Institute of Technology

Court e sy of Bill Goddard. Used wit h permis s i on.

Adri C.T. v. Duin

California Institute of Technology 35

ReaxFF: A reactive force field

E system

 E bond

 E vdWaals

 E Coulomb

 E val , angle

 E to rs

 E over

2- body

 E under

3- body 4-body

multi-body

Total energy is expressed as the sum of various terms describing indiv i dual chemical bonds

All expressions in terms of bond order

All interactions calculated between ALL atoms in system…

No more atom typing: Atom type = chemical element

Example: Calculation of bond energy

E system



 E vdWaals

 E Coulomb

 E val , angle

 E tors

 E over

 E under

E   D  BO

 e x p  p

 1  BO p be , 1  

bo n d e ij

 b e , 1

ij 

Bond energy between atoms i and j does not depend on bond distance

Instead, it depends on bond order

Bond order functions

BO goes smoothly from 3-2-

1-0

(1)

(1)

(3)

(2)

Fig. 2.21c in Bueh ler, Ma r k us J. At omisti c Modeling

of Materials Failure . Sp ring er, 2008. © Spring e r. All righ t s res e rved. This c o nt ent is exclud ed from our Creat i ve Commons lic e n s e. For more infor m at ion, s ee http:/ /ocw. m it.edu/fai ruse .

(2) ( 3)



  

    

  

  



B O ij

 ex p   







ex p  

 

  r ij 

r

  ex p  

   

  r ij  

r

 

  0   

 0  

Characteristic bond distance

All energy terms are expressed as a function of bond orders 38

Illustration: Bond energy

Imag e removed du e t o co pyrigh t restric t ion s.

Pleas e s e e slid e 10 in van Duin, Ad ri. "Dish i n g Out th e D i rt on R e axFF. " ht tp:/ / w w w . wag . cal t ech.edu/home/dui n/FFgroup/Di r t.ppt .

vdW interactions

E sy st e m

 E bond

 E v dW aal s 

E C oul om b 

E va l , angl e 

E tors 

E ov e r

 E unde r

 Accounts for short distance repulsion (Pauli principle orthogonalization) and attraction energies at large distances (dispersion)

 Included for all atoms with shielding at small distances

Imag e removed du e t o co pyrigh t restric t ion s.

Pleas e s e e slid e 11 in van Duin, Ad ri. "Dish i n g Out th e D i rt on R e axFF. " ht tp:/ / w w w . wag . cal t ech.edu/home/dui n/FFgroup/Di r t.ppt .

Resulting energy landscape

Contribution of E bond and vdW energy

Sour c e : van Dui n , C. T . Adri , et al. "Re axFF: A Reac ti ve Forc e Fi el d for Hydro ca r b o n s. " Jo ur nal of Phys i c al Chemi s t ry A 10 5 (2 001) . © A me ri c a n Chemic al Soc i e t y. All rig h ts res e rved . T h is content is ex cl uded from our Cr eative Common s licen se. For more in fo r m a t io n , se e h t t p ://ocw.m i t .e d u / f a i r u se .

Current development status of ReaxFF

: not currently described by ReaxFF

Period ic tab l e c o ur t e sy of Wikime dia Comm on s .

42

Mg-water interaction: How to make fire with water

Imag es removed du e t o copyrigh t restric t ion s.

Imag es from video sh owing explosive reac tio n of magne sium, silver nitra t e, and water, whi c h can be accessed here: http:/ / w w w . youtu be.com/ w atch?v=Q T K i v MVUcqE .

Mg

http://video.google.com/vi deoplay ?docid= 4697 996 292 9490 459 21& q= magnes iu m+ w a ter&total=46&s t ar t=0&nu m=50&so=0&ty pe=search&plindex = 0

Mg – w ater interaction – R eaxFF MD simulation

MIT OpenCou rse Wa re ht t p :// ocw.mit.edu

3.021 J / 1.021J / 10. 3 33J / 1 8 .36 1 J / 22.0 0 J Introd uctio n to Mo de li ng a n d Sim u lati o n

Sp r i ng 2012

F o r in fo r m a t ion about c i ting these ma te ria l s o r our Ter m s o f use , vis i t htt p :// oc w. mit.edu/ t e rms .