Reactive-Empirical-Force-Fields--Home-Page-Materials-反应力领域-首页材料课件.ppt

1、Reactive Empirical Force FieldsJason Quennevillejasonqlanl.govX-1:Solid Mechanics,EOS and Materials PropertiesApplied Physics DivisionLos Alamos National LaboratoryTimothy C.Germann,Los AlamosAlejandro Strachan,PurdueAdri C.T.van Duin,Caltech William A.Goddard III,CaltechAlexei A.Stuchebrukhov,UC Da

2、vis2019 Summer School on Computational Materials ScienceJuly 31-August 11,2019 University of IllinoisMotivationEmpirical force fields are used in biology,chemistry,physics and materials science to calculate the potential energy surface and atomic forces.Most,like CHARMM and AMBER,assume the same ato

3、mic connectivity(molecular composition)throughout simulation.No Chemistry!Straightforward solution:ab initio or QM/MM (up to 300 atoms for QM system)For materials simulation,we may want 10s of 1000s to millions of atoms and as much as a nanosecond of simulation time.Need a more efficient method!Empi

4、rical Force FieldsEmpirical force fields contain potential energy functions for each atomic interaction in a molecular system.Bond Stretch:Bond Bending:Bond Torsion:Parameters can be taken from experiment(e.g.,vibrational spectroscopy)or from ab initio quantum chemistry calculations.2stretch012RVkRR

5、2bend012Vktorsion1cos;1nVVsns Non-Bonded potentials give the intermolecular interactions:Coulomb:van der Waals:Parameters obtained through ab initio quantum chemistry and liquid simulations.e.g.OPLS(optimized potentials for liquid simulations,W.L.Jorgensen and J.Tirado-Rives,J.Am.Chem.Soc.110,1657(1

6、988).)Empirical Force FieldsLennard-Jones126ijijijijACVRR2Coulombijijq q eVR126;4;4ijijiiiijijiiiAAAACCCC Empirical Valence Bond(EVB)EVB attempts to combine empirical potential energy functions with valence bond ideas to describe chemical reactions efficiently and accurately.EVB Applications Proton

7、transport in aqueous acid (CPL,284,71(98);JPCB,102,5547(98)Aqueous acid-base reactions (JPCA,105,2814(01)Enzyme catalysis (Warshel)Nucleophilic substitution reactions Good Introduction:Computer Modeling of Chemical Reactions in Enzymes and Solutions,A.Warshel Wiley-Interscience(02/01/2019)EVB:introd

8、uctionEVB starts with a NN potential energy matrix:N diabatic states(diagonal)N(N-1)couplings(off-diagonal)Each diabatic state looks like a configuration in a standard non-reactive force field.Off-diagonal coupling elements:interaction between each diabatic state and the N-1 remaining states.Diagona

9、lize V adiabatic states.The minimal value is the ground state.coupling termadiabatic ground statediabatic statesCalculation of Forces0022nnmnnmnn mHFRHVaa aRR Diagonalizing the NN EVB matrix yields the ground state as a linear combination of diabatic states.If an is the set of corresponding coeffici

10、ents,the forces can be calculated using the Hellman-Feynman theorem:EVB:diagonal matrix elementsBecause we need to treat bond breaking and formation,the Bond Stretch should be anharmonic:Bending and Torsional potentials can be as before:System-environment interactions treated with standard non-bondi

11、ng potentials:2stretch01 expeVDRR2stretch012RVkRR2bend012Vktorsion1cos;1nVVsns Lennard-Jones126ijijijijACVRR2Coulombijijq q eVR126;4;4ijijiiiijijiiiAAAACCCC Interaction between EVB StatesSystem-system non-bonding interactions more complicated due to the potential for chemical reaction.A functional f

12、orm more flexible than Coulomb+Lennard-Jones is required.The intermolecular interactions(part of the diagonal element)and the coupling terms(off-diagonal)must be parametrized together in order to describe the reaction correctly.In the activated complex,the favorable interaction between the two state

13、s is controlled by the intermolecular interaction.It is normally written in terms of the distance between the two reactant centers.The reaction barrier is controlled by coupling term.This term is generally a function of the reaction coordinate.1.2 1.2 0.97 0.97 0.97 0.97 1161181161181091141.22 1.22

14、0.96 0.96 103103101Optimized geometries of the H2OHOH2+(left)and HOHOH(right)complexes,obtained from first principles(MP2/aug-cc-pVTZ).Application:Proton Transfer in WaterHHOOHHHHHOOHHHHHOOHHHDiabatic States:Adiabatic State:EVB of H3O+H2O Proton Transfer22222222H OH OH OH OO-Hintra0122H OH OH-O-H011

15、exp12eiijjVDRRk3233333332234H O/H OH Ointerdamp4O-O11O-O,H OO-OO-OO-O00O-O,H OO-O,H OO-H,H OO-HO-HO-H00O-H,H OO-H,H O2expexp/expexp/.ijijijqq eeVcRRRDRRRLDRRRL+3333+3332H OH OH OH OO-Hintra0132H OH OH-O-H011exp12eiijjVDRRkH3O+H2O Proton Transfer:Diagonal ElementsOOOO,XYXYI JVP RQ ROO(1)(1)OOOO(1)2-H

16、-H-0-H-(2)(2)OOOO(2)2-H-H-0-H-,expexpXYXYXYXYXYXYXYXYXYXYXYP RVRRVRROOOOOO-H-tanh,011tanh,2XYXYXYXYQ RRRHOO12tXYRRRH3O+H2O Proton Transfer:Coupling ElementsHHOOHHHEVB vs Ab Initio for H3O+/H2O2.4 2.6 2.8 EVBMP2/aug-cc-pVTZEVB SummaryVery good for systems with small number of possible reactionsReacti

17、on barriers are treated explicitlyOffers an empirical description of chemical reactionsGives mixing of diabatic states during reactionCan be difficult to parametrize intermolecular potentials and couplingsLimitation on number of states due to diagonalization(cubic scaling)ReaxFFBond-Order potential,

18、developed at CalTech by Adri van Duin and Bill GoddardPotential parametrized using ab initio calculations(B3LYP/6-31G*)on a“training set”of reactionsWhy bond-order based?non-reactive potentials have atom-types that define connectivity Applications:High Explosives,Propellants,Catalysis,Fuel Cells,Cor

19、rosion,Friction,etc.HN NHN N+H2Background ReferencesBond Order/Bond Length relationship Pauling,J.Am.Chem.Soc.,69,542(1947).Reactive Empirical Bond Order(REBO)Johnston,Adv.Chem.Phys.,3,131(1960).Johnston,Parr,J.Am.Chem.Soc.,85,2544(1963).Other Bond-Order Potentials Tersoff,Phys.Rev.Lett.,56,632(1986

20、);Tersoff,Phys.Rev.Lett.,61,2879(1988).Brenner,Phys.Rev.B,42,9458(1990).Brenner,et al,J.Phys.:Condens.Matter,14,783(2019).ReaxFF van Duin,Dasgupta,Lorant,Goddard,J.Phys.Chem.A,105,9396(2019).Strachan,Kober,van Duin,Oxgaard,Goddard,J.Chem.Phys.,122,054502(2019).User Manual:wag.caltech.edu/home/duin/r

21、eax_um.pdfReaxFF allows for computationally efficient simulation of materials under realistic conditions,i.e.bond breaking and formation with accurate chemical energies.Due to the chemistry,ReaxFF has a complicated potential energy function:bondvalency angletooverpotentialunderpenaltylone pairconjug

22、ationrsionH-bondvdWaalsCoulombEEEEEEEEEEEEReaxFF Potential Energy FunctionCharge equilibration:EEM(Mortier,et al,JACS,108,4315,(86).)HCCHExample:Acetylene Bond Order goes smoothly from 0 1 2 3 as C-C Bond Length shortens from large distance to 1.0 Bond Order,Bond Energy000BOexpexpexpijijijijrrrrrrno

23、t explicitly a function of bond distancebe,1bondebe,1BOexp1BOpijijEDp Bond OrderC-C Distance/Bond orders adjusted to get rid of unphysical bonds.Bond Order CorrectionsOver-and under-coordination of atoms must be avoided.Energy penalty added to the potential energy function for the case where an atom

24、 has more bonds than its valence allows.e.g.,Carbon cant have more than 4 bonds;Hydrogen no more than 1If an atom is under-coordinated,the stabilization of bonding should be used if possible.Bond Angles,Bond TorsionBond Angles and Torsions are intimately tied to the bond types.With a bond order pote

25、ntial,angles and torsions must be written in terms of the bond order.Angle and torsion energy terms must 0 as B.O.0.angle(,)ijkijjlEfBO BOtorsion(,)ijklijjkklEfBO BOBOSee J.Phys.Chem.A,105,9396(2019)for full potential form.CCH 120CCH=180bondangletorsionlone pairconjugationH-bo ndvdWaalsCoulombEEEEEE

26、EEELone Pair Electrons,ConjugationThe creation or reaction of lone-pair electrons should be assigned an energy term.Elone pair corresponds to an energy penalty for having too many lone pairs on an atom(i.e.,overcoordination)Conjugated systems should have added stabilization.Econj has maximum contrib

27、ution when successive bonds have bond-order values of 1.5.bondangletorsionlone pairconjugationH-bo ndvdWaalsCoulombEEEEEEEEE implicit in AMBER/CHARMM-like potentialsthrough atom typeHydrogen bonding extremely important in biological systems but also in many organic solids.X H YHydrogen bonds are cal

28、culated between group X-H and Y,where X and Y are atoms known to form H-bonds(e.g.,N,O)The H-bond energy term is written in terms of the bond-order of X-H,the distance between H and Y,as well as the X-H-Y angle.Can be an expensive part of the calculation because many acceptor(Y)atoms could be availa

29、ble for any given X-H group.All interactions must be calculated out to a cutoff distance(10)in order to remain consistent from timestep to timestep.Hydrogen BondingangleXHY HXH Y(,)Ef BOR-Short-range Pauli Repulsion-Long-range attraction(dispersion)-Coulomb forcesvan der Waals and Coulomb terms are

30、included for all atom pairs(whether bonded or non-bonded)!This avoids changing the potential when chemistry occurs.Such alterations,which are natural in the EVB formalism,would be awkward in ReaxFF.Shielding included for both Coulomb and van der Waals in order to avoid excessive interaction between

31、atoms sharing bond and/or bond angle.Non-Bonded InteractionsEnergy/kcal mol-1Interatomic Distance/Charge EquilibrationThe charge on an atom depends on the molecular species:Atomic charges are adjusted with respect to connectivity and geometry.Many QEq methods available.ReaxFF uses Electronegativity

32、Equalization Method(EEM:Mortier,et al,JACS,108,4315,(86).)The desired charge distribution is that which minimizes,Final Coulomb energy from screened potential all atom-pairs calculatedCoulomb1/333(1/)ijijijq qECrHN NHN N+H2212ijiiiiijiijq qEqJ qrReaxFF/Ab Initio ComparisonReaxFF can decribe a wide v

33、ariety of chemical reactions.Strachan,et al,JCP,122,054502(05).e.g.,unimolecular decomposition of RDX103,680 atoms(4320 molecules,121215 unit cells)256 processorsTATB768 atoms(32 molecules,224 unit cells)Interested in chemical reaction dynamics of high explosives(HE)under shock conditions Want as bi

34、g a system(105 to 106 atoms)as possible in order to study the spread of reactions,temperature distribution,carbon-clustering,etc.Application:HE at High T,PTATBN2H2OCO2Tinitial=1700 K Tfinal=3200 KDecomposition of TATB at High T GRASP(General Reactive Atomistic Simulation Program)developed at Sandia

35、National Lab by Aidan P.Thompson.Objective:Parallel scalable MD code(C+)which enables implementation of a wide range of force field types,particularly reactive force fields,including ReaxFF.CPU Time per timestep:serial code:2.8 seconds parallel code:0.6 seconds (32 CPUs)System size limits:serial cod

36、e:5000 atoms parallel code:500,000 atoms (510 CPUs)ASC Flash:2.0-GHz procs 8 GB memory per nodeReaxFF in ParallelReaxFF SummaryCan simulate chemistry for a wide range of materials significantly faster than ab initio and semi-empirical methodsAccuracy similar to semi-empirical methodsHydrocarbons,CHNO explosives,silicon oxides,etc.Main limitation is governed by the size of reaction training setUsed extensively for explosives under extreme conditions-many possible reactionsSimulation sizes up to a half million atomsTime/Iteration(seconds)Natoms