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 (cont’d)

Lecture 9

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 m odeling paradigm, difference between continuum and at omistic approac h , c a se st ud y: diff us ion)

 Lecture 3 : Basic statistical mechani cs – p roperty calculation I (property calc ulat ion: microsc o pi c stat es vs. macro scopic properties, ensembles, probability density and partition function)

 Lecture 4 : Property calcula tion II (Mont e Carlo, a d v a nc e d prope rty calcul ati on, introduction to chemical interactions)

 Lectu r e 5: Ho w to m odel ch em i c al interactions I ( e xa mp le: mo vie o f coppe r def ormation/ di sloc ati o ns, etc.)

 Lectu r e 6: Ho w to m odel ch em i c al interactions I I (EAM, a bit of ReaxFF—chemical reactions)

 Lectu r e 7: Ap pli c atio n – M D s imulation of materials failure

 Lectu r e 8: Ap pli c atio n – Reacti ve potentials and applications

 Lectu r e 9: Ap pli c atio n – Reacti ve p o tentials and applications (cont’d)

Lecture 9: Reactive potentials and applications (cont’d)

Outline:

1. Notes on fracture application

2. Closure: ReaxFF force field

3. Hybrid multi-paradigm fracture models

Goal of today’s lecture:

 Remarks: Modeling of fracture and relation to diffusion problem

 New potential: ReaxFF, to descri be complex chemistry (bond break ing and formation)

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

1. Notes on fracture application

Consider for pset #2

Brittle fracture mechanisms: fracture is a multi- scale phenomenon, from nano to macro

R e print e d by permis sio n from Macmillan Publish e rs L t d: Nature .

S o urc e : Bueh ler, M., and Z. Xu. "Materia ls S c ien c e: Mind the H e lica l C r ack." Nat u re 464, no. 7285 (2010): 42 -3. © 2010.

Limiting speeds of cracks: linear elastic

continuum theory

Mother-daughter mechanism

Super- Sub-Rayleigh Rayleigh

Mode II

Intersonic

Supersonic

Mode I

C r C s

C l

Limiting speed v

Subsonic

Supersonic

Mode III

C s

C l

Limiting speed v

Linear Nonlinear

c l  ~

c s  ~

c r  0 . 92 c s

Image by MIT OpenCou r seWare.

• C racks can not exceed the limiting speed given by the corresponding wave speeds unless materia l behavior is nonlinear

• C racks that exceed limiting s peed would produce energy (physically impos sible - linear elastic continuum theory )

Subsonic and supersonic fracture

 Under certain conditions, material non linearities (that is, the behavior of materials under large deformation = hyperelasticity) becomes important

 This can lead to different limiting speeds

than described by the model introduced above

Deformation field near a crack

 ( r ) ~ 1

E large (stiff)

c l  ~

small deformation

E small (soft)

Hyperelasticity

Stra in

Image by MIT OpenCou r seWare.

Far away from crack tip

Energy flux concept

L energy

Image by MIT OpenCou r seWare.

Characteristic energy

length

(energy from this distance needs to flow to the crack tip)

L large (stiff)

Crack tip

L energy

 1  Supersonic cracking

Energy flux reduction/enhancement

K domina nce zone

Chara cteristic energy length

Hyperelastic zone

Fracture process zone

New

L energy

Image by MIT OpenCou r seWare.

Energy flux related to wave speed: high local wave speed, high energy flux, crack can move faster (and re verse for low local wave speed)

Physical basis for subsonic/supersonic fracture

 Changes in energy flow at the crac k tip due to changes in local wave speed (energy flux higher in materials with higher wave speed)

 Controlled by a characteristic length scale

L energy

Buehler et al ., Nature , 2003

R e print e d by permis sio n from Ma cmill an Publish e rs L t d: Nature.

S o urc e : Bueh ler, M., F. Abraham, and H. Gao. "H ypere l a s t i c i ty Governs Dynamic Fractur e at a Critica l L e ngth S ca l e." Nature 42 6 (2003): 1 41-6. © 20 03.

11

Integ r ati o n (BCs, ICs)

Characteristic time of diffusion

PD E (continuum equilibrium)

 c 

 t D dx 2

PDE

d 2 c

(Fick’s l a w)

Integ r ati o n (BCs, ICs)

Crack limiting speed

Continuum appr oach (distinct PDE)

Integ r ati o n (BCs, ICs)

© sourc e u n known. All right s res e rved. This con t en t is exclud ed from our Crea t i ve Com m o ns lic e n s e. For more i n format ion, see ht tp:/ /ocw. m it.edu/fai ruse .

Characteristic time of diffusion

Atomistic approach (same PDE)

Integ r ati o n (BCs, ICs)

Crack limiting speed

12

2. Closure: ReaxFF force field

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

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! 15

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

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 19

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

Adri C.T. v. Duin

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

California Institute of Technology 23

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 26

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 .

Mg-water interaction: How to make fire with water

Video st ill s removed due t o copyri ght restri ct ions; watch the vi deo now: ht tp:/ / w w w . yout ube.com/ watch?v=Q T Ki v M VUcqE .

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

3. Hybrid multi-paradigm fracture models

Focus: model particular fracture properties of silicon (chemically complex material)

Fracture of silicon: problem statement

Pair potential insufficient to describe bond break ing (chemical complex i ty)

Ima g e cour t e sy of NASA.

Multi-paradigm concept for fracture

 ( r ) ~ 1

Need method

“good” for describing rupture of chemical bonds

r

Need method “good” for elastic properties (energy storage)

Image by MIT OpenCou r seWare. 35

What about combining different potentials?

 Empirical models: mathematical functions with parameters (fitted to experiment or quantum mechanics)

 Pair potentials (LJ, Mor s e, Buck., har monic) (le c ture 5)

 Embedded atom models/effective medium theories

 Multi-body potentials (e.g. Tersoff, C H A R MM, etc.)

(lecture 9 and following)

 Reactive potentials (ReaxFF) ( lecture 9)

 Semi-empirical models (explicitly note electronic structure)

 Tight binding

 MINDO (= Modified Inter m ediate Neglect of Differ ential Overlap), NINDO (= Inter m ediate Neglect of Differ ential Overlap)

 Quantum mechanical models: Start from Schroedinger’s equation (and make approximations to be able to solve it)

 Quantum chemistry (Har tree-Fock)

 Density Functional Theory

 Quantum Monte Carlo

“good” for elastic properties (energy storage)

“good” for +

describing rupture of chemical bonds

=

“multi-paradigm

model”

Concept: concurrent multi-paradigm simulations

Crack tips, defects (dislocations)

nonr eac ti v e ato m i s ti c

Reax FF

nonr eac ti v e ato m i s ti c

• Multi-paradigm appr oach : combine different computational methods (different resolution, accuracy..) in a single computational domain

• Decomposition of domain based on suitability of different approaches

• Example : concurrent FE- atomistic-ReaxFF scheme in a crack problem (crack tip treated by ReaxFF) and an interface

Interfaces (oxidation, grain boundaries,..)

problem (interface treated by ReaxFF).

Concurrent multi-paradigm simulations: link nanoscale to macroscale

Concurrent coupling : use of multiple force fields within one simulation domain

Simulation Geometry: Cracking in Silicon

Potential for covalent bonds, not suitable for bond breaking

Fig. 2 in Bueh ler, Markus J., et al. "Multiparad i gm Modelin g of Dynamica l Crack Propagat ion in S ili con Using a React i ve Force Fiel d. " Ph ysical Revie w Let t e rs 96 (2006): 095505. © APS. All rig h ts 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 .

• C ons ider a crack in a silic on crystal under mode I loading.

• P eriodic boundary conditions in the z- direction (corresponding to a plane strain case).

Cracking in Silicon: Hybrid model versus Tersoff based model

Pure T ersof f

Hybrid R eaxFF- T ersof f

Image by MIT OpenCou r seWare.

Conclusion: Pure Tersoff can not describe correct crack dynamic s

How is the handshaking achieved?

Hybrid potential energy model (Hamiltonian)

Weights = describe how much

W i

100%

W 2

W 1

0%

R - R buf

R

R+R tr ans R+R tr ans +R buf

x

R eaxFF

T r ansition la y er

T ersof f

T r ansition region

T ersof f ghost atoms

R eaxFF ghost a toms

a particular FF counts (assigned to each atom)

To obtain forces:

F  

 U tot ( x )

 x

need potential energy

Image by MIT OpenCou r seWare.

Approach: handshaking via mixed Hamiltonians

transition region

42

Assigning weights to atoms

W i

100%

W 2

W 1

0%

R - R buf

R

R+R tr ans R+R tr ans +R buf

x

R eaxFF

T r ansition la y er

T ersof f

T r ansition region

T ersof f ghost atoms

R eaxFF ghost a toms

Image by MIT OpenCou r seWare.

Percentage ReaxFF Percentage Tersoff (relative contribution to total energy)

100% … 100% 70% 30% 0% … 0 %

0% … 0% 30% 70% 100% … 100%

43

Force calculation

W i

100%

W 2

W 1

0%

R - R buf

x

R

R+R tr ans R+R tr ans +R buf

T r ansition region

Potential energy

U tot

 U ReaxFF  U Tersoff

 U ReaxFF  Tersoff

U ReaxFF  Tersoff  w ReaxFF ( x ) U ReaxFF  ( 1  w ReaxFF ) U Tersoff

Recall:

F    U

 x

Image by MIT OpenCou r seWare.

w ReaxFF ( x )  w Tersoff ( x )  1  x

w ReaxFF is the weight of the reactive fo rce field in the handshaking region.

F     w

( x ) F

 ( 1  w

) F    w ReaxFF  U

 U  

ReaxFF  Tersoff

  ReaxFF

ReaxFF

ReaxFF

Tersoff

 x ReaxFF

Tersoff  

D. Sen and M. Buehler, Int. J. M u ltiscale Com p ut. Engrg ., 2007

Hybrid Hamiltonians – force calculation

F     w

( x ) F

 ( 1  w

) F  

ReaxFF  Tersoff

  ReaxFF

ReaxFF

ReaxFF

Tersoff

≈ 0

≈ 0

Slowly varying weights (wide transition region):

 w ReaxFF /  x  0

If U ReaxFF  U Tersoff  0

( i.e., both force fields have similar energy lands cape)

Simplified result: can interpolate forces from one end to the other

D. Sen and M. Buehler, Int. J. M u ltiscale Com p ut. Engrg ., 2007

Energy landscape of two force fields

• S chematic showing the coupling of reactive and nonreactive potentials

• A t small dev iations , energy landsc ape is identic al in nonreac tive

U and reactive models

x

U ReaxFF  U nonreactiv e  0

Summary: hybrid potential energy model

W i

100%

W 2

W 1

0%

R - R buf

R

R+R tr ans R+R tr ans +R buf

x

R eaxFF

T r ansition la y er

T ersof f

T r ansition region

T ersof f ghost atoms

R eaxFF ghost a toms

Image by MIT OpenCou r seWare.

F ReaxFF  Tersoff

  w ReaxFF ( x ) F ReaxFF

 ( 1 

w ReaxFF ) F Tersoff 

w ReaxFF ( x ) 

w Tersoff

( x )  1  x

Fracture of silicon single crystals

Us e multi-paradigm scheme that combines the Tersoff potential and ReaxFF

0 ps

7.0 ps

14.0 ps

Reactive regio n (red) is mo ving with crack tip.

R eaxFF T ersof f

Image by MIT OpenCou r seWare.

Image by MIT OpenCou r seWare.

Quantitative comparison w/ experiment

“empirical” FFs

Fi g. 1c i n Buehle r, M . , et al . "Thre s hol d Crack Sp eed Con t rols Dynamical Frac t u re of Si licon Si ngle Cr ystals." P h ysical Review Let t e rs 99 (2007). © APS. All righ ts r e s e rved. This c o n t en t is excluded from ou r C r ea t ive C o mm on s li c en s e. F o r mor e in f o rmat i o n , s e e http:/ /ocw. m it.edu/fai ruse .

Load: normalized by critical energy release rate to initiate fracture

Crack dynamics

Imag e removed du e t o copyrigh t restric t io n s . Pleas e s e e: Fig. 2 in Bueh ler, M., et al. " T h resh old Crack Speed Con t rols Dy nami cal Fracture of Sil icon Single Cr yst a l s ." Ph ysical Review Lett ers 99 (2 007).

Crack speed: O(km/sec)

=O(nm/ps) (well in reach with MD)

Atomistic fracture mechanism

Fracture initiation and instabilities

Fracture mechanism: tensile vs. shear loading

Shear (mode II) loading:

Crack branching

Tensile (mode I) loading:

Straight cracking

Image by MIT OpenCou r seWare.

M.J. Buehler, A. Cohen, D. Sen, Journal of Algorithms and Computational Technology, 2008

Fracture mechanism: tensile vs. shear loading

Shear (mode II) loading: Crack branching Tensile (mode I) loading: Straight cracking

Image by MIT OpenCou r seWare.

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

Pleas e s e e figur e s in Bu eh ler, M. J . , A. Coh e n, and D. Sen. " Mult i-paradig m Modeling of Fractur e of a Silic on Sing le Cr yst a l Un d e r Mode I I Sh ear L o ad ing . " Journal of Alg o rithms and Computati o nal Techn o l o gy 2 (2008): 203-21.

Summary: main concept of this section

 Can combine different force fields in a single computational domain = multi-paradigm modeling

 Enables one to combine the strengths of different force fields

 Simple approach by interpolating force contributions from indiv idual force fields, use of we ights (sum of weights = 1 at all points)

 ReaxFF based models quite successful , e.g. for describing fracture in silicon, quantitative agreement with experimental results

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 .