Ch121a Atomic Level Simulations of Materials and Molecules

1. Ch121a Atomic Level Simulations of Materials and Molecules BI 115 Hours: 2:30-3:30 Monday and Wednesday Lecture or Lab: Friday 2-3pm (+3-4pm) Lecture 11, April 25, 2011 Reactive Force Fields – 1: ReaxFF William A. Goddard III, wag@wag.caltech.edu Charles and Mary Ferkel Professor of Chemistry, Materials Science, and Applied Physics, California Institute of Technology Teaching Assistants Wei-Guang Liu, Fan Lu, Jose Mendoza, Andrea Kirkpatrick

2. Bridging the Gap between QM and MD to describe reactive processes ranging from combustion to catalysis, fuel cells, nanoelectronics, and shock induced chemistry ReaxFF: first principles force field 471. ReaxFF: A Reactive Force Field for Hydrocarbons A.C.T. van Duin, S. Dasgupta, F. Lorant and W. A. Goddard III J. Phys. Chem. A, 105: (41) 9396-9409 (2001) 514. ReaxFF sio Reactive Force Field for Silicon and Silicon Oxide Systems Adri C. T. van Duin, Alehandro Strachan, Shannon Stewman, Qingsong Zhang, Xin Xu, and William A. Goddard, III J. Phys. Chem. A, 107, 3803 (2003) • 533. Shock waves in high-energy materials: The initial chemical events in nitramine RDX Strachan A, van Duin ACT, Chakraborty D, Dasgupta S, Goddard WA Physical Review Letters, 91 (9): art. no. 098301 (2003)

3. Motivation: Design Materials, Catalysts, Pharma from 1st Principles so can do design prior to experiment time ELECTRONS ATOMS GRAINS GRIDS simulations real devices and full cell (systems biology) hours millisec nanosec picosec femtosec Continuum (FEM) Micromechanical modeling Protein clusters MESO MD Deformation and Failure Protein Structure and Function QM distance Å nm micron mm yards To connect 1st Principles (QM) to Macro work use an overlapping hierarchy of methods (paradigms) Need accurate force fields (potentials) with parameters derived fromfirst-principles QM Big breakthrough making FC simulations practical: reactive force fields based on QM Describes: chemistry,charge transfer, etc. For metals, oxides, organics. Accurate calculations for bulk phases and molecules(EOS, bond dissociation) Chemical Reactions (P-450 oxidation)

4. r E = kA(A-A0) A E = kr(r-r0) Ordinary Force Fields Bonds, angles, torsions described as elastic springs Fixed charges, Empirical vdw nonbond terms, Bonds cannot be broken, making the model unsuitable for modeling reactions. Examples: MM3, Dreiding, Amber, Charmm, Gromos, UFF ReaxFF: Allow bonds to break and form and describe barriers for reactions. All parameters from quantum mechanics, no empirical data

5. First Principles Reactive force fields: strategy • Describe Chemistry (i.e., reactions)of molecules • Fit QM Bond dissociation curves for breaking every type of bond • (XnA-BYm), (XnA=BYm), (XnA≡BYm) • Fit angle bending and torsional potentials from QM • Fit QM Surfaces for Chemical reactions (uni- and bi-molecular) • Fit Ab initio charges and polarizabilities of molecules • Pauli Principle: Fit to QM for all coordinations (2,4,6,8,12) • Metals: fcc, hcp, bcc, a15, sc, diamond • Defects (vacancies, dislocations, surfaces) • cover high pressure(to 50% compression or 500GPa) • Generic: use same parameters for all systems(same O in O3, SiO2, H2CO, HbO2, BaTiO3) Require that One FF reproduces all the ab-initio data (ReaxFF) Most theorists (including me) thought that this would not be possible, but we claim to have achieved and validated it for many systems

6. Many Chemical processes: bond breaking for 5000millions atoms. This is far too large for DFT Solution: ReaxFF first principles reactive force field short distance Pauli Repulsion + long range dispersion (pairwise Morse function) Valence energy Electrostatic energy • Based completely on First Principles QM (no empirical parameters) • Valence Terms (EVal) based on Bond Order: dissociates smoothly • Bond distance  Bond order  Bond energy • Forces depend only on geometry (no assigned bond types) • Allows angle, torsion, and inversion terms (where needed) • Describes resonance(benzene, allyl) • Describes forbidden(2s + 2s)and allowed(Diels-Alder)reactions • Atomic Valence Term(sum of Bond Orders gives valency) • Pair-wiseNonbond Terms between all atoms (no “bond” exclusions) • Short range Pauli Repulsion plus Dispersion(EvdW) • QEq Electrostatics allows charges to flow depending on environment and external fields

7. ReaxFF reactive force field for Reactive Dynamics (RD) Adri van Duin Allow bonds to break and form, describe barriers for reactions. All parameters from quantum mechanics (no empirical data) ReaxFF describes reactive processes (from oxidation to combustion to catalysis to shock induced chemistry)for 1000s to millions atoms We use ReaxFF to prepare the structures of complex heterogeneous systems by processes similar to experimental synthesis (DLC, SiO2/Si)

8. Sigma Bond 1.5Å - 2.3Å First Pi Bond 1.2Å - 1.75Å Second Pi Bond 1.0Å - 1.4Å Bond Order Dependent Potentials

9. Bond Order Dependence on r (CC)   s

10. Bond distance bond order C-C bond   Bond order s Bond distance (Å) Bond order bond energy Bond distance bond order bond energy Use general functional form. Determine parameters from fitting QM bond breaking for many single, double, and triple bonded systems. Bond distance (Å)

11. Van der Waals Include vdW for 1-2 and 1-3 interactions since bonded atoms may dissociate and nonbonded atoms may bond during the dynamics We use a Morse function rather than LJ12-6 (which is too stiff in the inner wall region Or exponential-6 which has problems for small R Here f7 was introduced to modify interactions for small R This was a mistake and should be eliminated

12. vdW Energy

13. Coulomb interacions Electronegativity Hardness Include Coulomb for 1-2 and 1-3 interactions since bonded atoms may dissociate and nonbonded atoms may bond during the dynamics This f7 function provides shielding so that the Coulomb potential goes to a constant, Cqiqj gw at small R, where gw is related to the size of the atoms. ReaxFF uses the Electron Equilibration Method (EEM) of Mortier to determine charges rather than rather than QEq Mortier, W. J.; Ghosh, S. K.; Shankar, S. J. Am. Chem. Soc. 1986,108, 4315; Janssens, G. O. A.; Baekelandt, B. G.; Toufar, H.; Mortier, W. J.;Schoonheydt, R. A. J. Phys. Chem. 1995, 99, 3251.

14. Discussion about charges in ReaxFF I believe that the use of a shielding denominator in the Coulomb energy and the use of EEM in ReaxFF are inconsistent and a fundamental mistake. We should use QEq. In QEq the shielding of charges on different centers is built in both in calculating the Coulomb energy and in calculating the charges. Separating them leads to potential problems for small R, where the charges may be affected by the 1/R form in the potential. In QEq this shielding is based on the covalent radius of the atom and hence need not be optimized.

15. Charge Equilibration (QEq) Charge Equilibration for Molecular Dynamics Simulations A. K. Rappé and W. A. Goddard III J. Phys. Chem. 95, 3358 (1991) I atomic interactions Keeping: Jij Hardness (IP-EA) Electronegativity (IP+EA)/2 1/rij rij ri0 +rj0 • Self-consistent Charge Equilibration (QEq) • Describe charges as distributed (Gaussian) • Thus charges on adjacent atoms shielded • (interactions  constant as R0) and include interactions over ALL atoms, even if bonded (no exclusions) • Allow charge transfer (QEq method) 2 Three universal parameters for each element: 1991: use experimental IP, EA, Ri; ReaxFF get from fitting QM

16. Experiment Mol. Eb re C2H6 90.4 - 1.533 C3H8 85.8 - 1.526 C4H10 86.5 - 1.532 CC Bond Dissociation Energy In ReaxFF the BO-BE term is monotonic and attractive, becoming a constant at small R This is balance by the vdW term which is large and repulsive at small R To finally give the normal bonding function.

17. Deviation of bond orders from saturation Over-coordination Coordination

18. Under-coordination

19. Angle

20. Calculation of 0 SBO2 = 2; SBO>2 SBO2 = 2-(2-SBO)5; 1<SBO<2 SBO2 = SBO 5; 0<SBO<1 SBO2 = 0; SBO  0 2j = j; j < 0 2j = 0; j  0

21. Eangle

22. Torsion • V2 diminishes rapidly when central bond deviates from BO=2 • Sin terms ensure this terms goes to 0 when the angles approach 180 • f4 accounts for over-coordinated central C-C atoms • f3 allows for proper dissociation behavior

23. f3

24. f4

25. V2

26. Etors

27. Conjugation

28. Heats of Formation for 40 compounds/radicals

29. Bond dissociation energies

30. Bond lengths and r0

31. RDX concerted ring breaking

32. RDX N-N dissociation

33. HONO elimination

34. DFT Calculation

35. EoS

36. EoS

37. Over-coordination energy and van der Waals Over-coordination penalty Over-coordination energy Energy (kcal/mol) Total bond order Shielded van der Waals Unshielded van der Waals Energy (kcal/mol) Shielded vdW (w= 0.8) C-C distance (Å)

38. Key features of ReaxFF • To get a smooth transition from nonbonded to single, double and triple bonded systems ReaxFF employs a bond length/bond order relationship1,2. Bond orders are updated every iteration. • Nonbonded interactions (van der Waals, Coulomb) are calculated between every atom pair, irrespective of connectivity. Excessive close-range nonbonded interactions are avoided by shielding. • All connectivity-dependent interactions (i.e. valence and torsion angles) are made bond-order dependent, ensuring that their energy contributions disappear upon bond dissociation. • ReaxFF uses a geometry-dependent charge calculation scheme (EEM) that accounts for polarization effects. 1:Tersoff, PRB 1988;2: Brenner PRB 1990

39. General rules ReaxFF MD-force field; no discontinuities in energy or forces even during reactions. User does not pre-define reactive sites or reaction pathways; ReaxFF automatically handles coordination changes associated with reactions. Each element is represented by only 1 atom type in the force field; force field determine equilibrium bond lengths, valence angles etc. from chemical environment.

40. ReaxFF uses generic rules for all parameters and functional forms ReaxFF automatically handles coordination changes and oxidation states associated with reactions, thus no discontinuities in energy or forces. User does not pre-define reactive sites or reaction pathways (ReaxFF figures it out as the reaction proceeds) Each element leads to only 1 atom type in the force field. (same O in O3, SiO2, H2CO, HbO2, BaTiO3)(we do not pre-designate the CO bond in H2CO as double and the CO bond in H3COH as single or in CO as triple, ReaxFF figures this out) ReaxFF determines equilibrium bond lengths, angles, and charges solely from the chemical environment. Require that One FF reproduces all the ab-initio data (ReaxFF) Most theorists (including me) thought that this would be impossible, hence it would never have been funded by NSF, DOE, or NIH since it was far too risky. (DARPA came through, then ONR, then ARO).

41. Catalysts: Pt-Co Fuel cell cathode, Pt-Ru FC anode VOx catalyzed oxidative dehydrogenation: propane to propene MoVNbTaTeOx ammoxidation catalysts (propaneacrylonitrile) Ni,Co,Mo catalyzed growth of bucky tubes Metal alloy phase transformations (crystal-amorphous) Si-Al-Mg oxides: Zeolites, clays, mica, intercalated polymers Gate oxides (Si-HfO2, Si-ZrO2, Si-SiOxNy interfaces) Ferroelectric oxides (BaTiO3) domain structure, Pz/Ez Hysteresis Loop of BaTiO3 at 300K, 25GHz by MD Decomposition of High Energy (HE) Density Materials (HEDM) MD simulations of shock decomposition and of cook-off MD elucidation of the origins of sensitivity in HE materials Reaction Kinetics from MD simulations (oxidations) ADP-ATP hydrolysis Enzyme Proteolysis Current applications of ReaxFF

42. Applications to energetic materials Sergey Zybin, Peng Xu, Yi Liu, Qing Zhang, Luzheng Zhang, Adri van Duin and William A. Goddard III • Force field development • Training sets for HE-materials • Treating organic crystals • Overview of past and ongoing applications • Predicting chemistry : cookoff of RDX, HMX and TATB and carbon cluster analysis • Prediction of sensitivity for HE materials • Energy release: ISP prediction for RDX/Al and hydrazine materials • New HE-materials: all-Nitrogen

43. Training set for nitramine potential Bond and angle distortion Reaction barriers C-N bond dissociation in H2C=NH Nitromethane C-N-O angle bending First-row elements Charge distributions Condensed phase

44. ReaxFF gives accurate description of complex chemical reactions: Decomposition Mechanism RDX (gas-phase) NO2 dissociation concerted New mechanism QM ReaxFF 8 membered ring HONO elimination Concerted, NO2 and HONO- dissociation pathways (part of the original training set) ReaxFF MD simulation found New unimolecular reaction, confirmed by QM, More important than concerted pathway Strachan, Kober, van Duin, Oxgaard and Goddard, J.Chem.Phys 2005

45. RDX shock simulations: MD with ReaxFF • Simulate High impact shock chemistry using MD simulations • 1st step: thermalize two semi-infinite slabs of RDX (1344 atoms/cell) • Add the shock velocity on top of the thermal velocities • Constant NVE molecular dynamics (adiabatic) ∞ Slab: 32 RDX molecules/cell on ∞ Slab:32 RDX molecules/cell Full-physics, full-chemistry simulations of shock induced chemistry

46. Hydrocarbon combustion and metal-oxide catalyzed hydrocarbon oxidation Adri, Kimberly Chenoweth, Sanja Pudar, Mu-Jeng Cheng, and wag Selective oxidation of propene using multi-metal oxide (MMO) catalysts Develop ReaxFF based on QM-data , use ReaxFF to perform high-temperature simulations on catalyst/hydrocarbon reactions First, establish that ReaxFF can describe non-catalytic hydrocarbon combustion Mixed metal oxide catalyst (BixMoyVzTeaOb)

47. Force field development: hydrocarbon oxidation QM ReaxFF Oxidation reactions QM: Jaguar/DFT/B3LYP/6-311G** Radical rearrangements Rotational barriers Angle strain H3C• + CO • H3C-C=O - total training set contains about 1700 compounds

48. Test ReaxFF CHO-description: oxidation of o-xylene • Oxidation initiates with OOH-formation • Final products dominated by CO, CO2 and H2O Consumed O2 CO2 H2O CO o-Xylene OOH OH • Exothermic reaction • Exothermic events are related to H2O and CO2 formation 2 o-Xylene; 70 O2 in 20x20x20 Angstrom box ReaxFF NVT/MD at T=2500K

49. o-Xylene oxidation: Detailed reaction mechanism O O O O O 2 2 2 2 2 H C O O H C H 2 C H 2 3 H O 2 frame 174 C H frame 128 C H C H 3 3 3 frame 175 O O H C = O H C = O H C O 2 2 2 O O H H O H O H C H 3 C H C H 3 3 frame 176 frame 176 frame 177 H O O O O O O O H O H C H C H frame 180 2 H C = O frame 179 2 O H 2 C H H C = O 2 O H frame 182 H H O H H O H H H O H H H H H H frame 205 C H C O frame 193 2 C O 2 C O 2 C O H O O H 2 2 O H H O Kimberly Chenoweth 2 frame 209 H H C O O H H O H H H O 2 O H O H C O H H H 2 O H H H H frame 232 frame 232 C O H H H O H 2 2 C O frame 234 C O O H 2 2 C O H C O C O 2 H H H H H H H H H H H H H H H H H O O O O O O O O O O O O O O O H H 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 O H 2 O O H C = O • Reaction initiation with HO2-formation • Dehydrogenation occurs at methyl-groups, not at benzyl-hydrogens • Only after H2C=O is formed and dissociated the benzene ring gets oxidized • Ring opens shortly after destruction of aromatic system • Ring-opened structure reacts quickly with oxygen, forming CO2, H2O and CO • ReaxFF gives sensible predictions that can be directly tested against QM

50. Determining the parameters for ReaxFF: MoOx ReaxFF diamond Simple cubic fcc A15 bcc Step 1: get ReaxFF for Mo metal Need to describe the complicated bonding in MoxOy polymorphs QM diamond Energy (kcal/mol) Mo17O47-crystal Simple cubic fcc A15 bcc Energy (kcal/mol)