ReaxFF:A Reactive Force Field for Hydrocarbons
Adri C.T.van Duin,†,|Siddharth Dasgupta,‡Francois Lorant,§and William A.Goddard III*,‡
Department of Fossil Fuels and En V ironmental Geochemistry,Drummond Building,Uni V ersity of Newcastle,
Newcastle upon Tyne NE17RU,United Kingdom,Materials and Process Simulation Center,Beckman Institute
(139-74),Di V ision of Chemistry and Chemical Engineering,California Institute of Technology,Pasadena,
California91125,and Institut Franc¸ais du Pe`trole,Geology and Geochemistry Research Di V ision,1-4A V enue
de Bois-Preau,92852Rueil-Malmaison,France
Recei V ed:December4,2000;In Final Form:March30,2001
To make practical the molecular dynamics simulation of large scale reactive chemical systems(1000s of
atoms),we developed ReaxFF,a force field for reactive systems.ReaxFF uses a general relationship between
bond distance and bond order on one hand and between bond order and bond energy on the other hand that
leads to proper dissociation of bonds to separated atoms.Other valence terms present in the force field(angle
and torsion)are defined in terms of the same bond orders so that all these terms go to zero smoothly as bonds
break.In addition,ReaxFF has Coulomb and Morse(van der Waals)potentials to describe nonbond interactions
between all atoms(no exclusions).These nonbond interactions are shielded at short range so that the Coulomb
and van der Waals interactions become constant as R ij f0.We report here the ReaxFF for hydrocarbons.
The parameters were derived from quantum chemical calculations on bond dissociation and reactions of small
molecules plus heat of formation and geometry data for a number of stable hydrocarbon compounds.We find
that the ReaxFF provides a good description of these data.Generally,the results are of an accuracy similar
or better than PM3,while ReaxFF is about100times faster.In turn,the PM3is about100times faster than
the QC calculations.Thus,with ReaxFF we hope to be able to study complex reactions in hydrocarbons.
1.Introduction
The accuracy and speed of modern quantum chemistry(QC) methods allow the geometries,energies,a
nd vibrational energies to be predicted quite accurately for small molecules.However, QC is not yet practical for studying the dynamic properties of larger molecules and solids.Consequently,it is useful to have accurate force fields(FF)to quickly evaluate the forces and other dynamical properties such as the effects of mechanical shock waves or diffusion of small molecules in polymer and mesoporous zeolites.Indeed,for hydrocarbons a number of FF, particularly MM3,1-3provide accurate predictions of geometries, conformational energy differences,and heats of formation. Generic FF such as DREIDING4and UFF5allow predictions for broad classes of compounds,particularly when coupled to charge equilibration6(Q Eq)or other methods for predicting charges.However,in general,these force fields do not describe chemical reactivity.An exception is the Brenner potential,7 which leads to accurate geometries for ground states of hydrocarbons,but is formulated in such a way that it can describe bond breaking.Unfortunately,the Brenner formalism does not include the van der Waals and Coulomb interactions that are very important in predicting the structures and properties of many systems.In addition,the actual potential curves for bond breaking and reactions are often quite poorly described with the Brenner potential.Generalizations of the Brenner FF have included such nonbonded forces,8but without repairing the fundamental problems in the shapes of the dissociation and reactive potential curves.
Two other bond-order-dependent force field methods are noteworthy.The Bond Energy Bond Order(BEBO)method was proposed by Johnston9,10based on Pauling’s relation between bond length and bond order.11The fundamental assumption behind this method is that the path of lowest energy on going from reactant to product is one that conserves total bond order. Originally used for the H+H2reaction surface,it is a very good approximation to more complicated empirical forms such as LEPS surface.12-13While it has recently been extended to more complex reactions,14-15it remains mainly useful for H atom transfer reactions in a linear collision geometry.
The VALBOND method formulated by Landis and colleagues is based on the strength functions of hybrid orbitals.The motivation comes from the need to fit large distortions in the softer angle terms of valence force fields,as well as describing multiple equilibrium angles in transition metal , the90°and180°angles in square planar complexes).Assuming that different ligand atoms,lone pairs,and radical electrons have implicit preference for p-character,VALBOND uses Lewis structure-based allocations to assign hybridizations and the geometries are obtained by minimizing defects in the hybrid orbitals.For a simple force field,it does remarkably well on structures and vibrational frequencies for a wide range of small molecules and transition metal complexes.16-18These methods, however,do not fully address the need to have full chemistry of the breaking and for
ming bonds,in addition to a proper description of the fully bonded equilibrium geometry of complex molecules.
In this paper,we develop a general bond-order-dependent potential in which the van der Waals and Coulomb forces are
*Author to whom correspondence should be addressed.E-mail:wag@
wag.caltech.edu.
†University of Newcastle.
‡California Institute of Technology.
§Institut Franc¸ais du Pe`trole,Geology and Geochemistry Research
Division.
|E-mail:  A.C.T.van-Duin@ncl.ac.uk.
9396J.Phys.Chem.A2001,105,9396-9409
10.1021/jp004368u CCC:$20.00©2001American Chemical Society
Published on Web09/22/2001
included from the beginning and the dissociation and reaction curves are derived from QC calculations.In spirit it is derived from the central force concept used earlier by spectroscopists 19but abandoned because a single harmonic potential between all atom pairs was inadequate for complex molecules.We have kept the central force formalism,where all atom pairs have nonbonded interactions,because it dissociates smoothly,but add local perturbations (bond,angle,torsion,etc.)to describe the complex molecules more accurately.We have attempted to obtain accurate descriptions of quantum phenomena such as resonance,unsaturated valences in radical systems,and chemical reactions.While the current work is restricted to hydrocarbons this approach is easily extended to any molecular system of any class of compounds.In a future paper we will report on our extension to CHNO-systems.Section 2describes the general form of the reactive force field (denoted ReaxFF)and the procedure for optimizing the parameters.Section 3presents the results for a number of systems.Section 4discusses these results,and Section 5presents the conclusions.2.Force Field
Similar to empirical nonreactive force fields,the reactive force field divides the system energy up into various partial energy contributions,as demonstrated by eq 1.
The potential energy functions associated with each of these partial energy contributions are described below.Tables 1-6list the parameters used in these potential functions.
2.1.Bond Order and Bond Energy.A fundamental as-sumption of ReaxFF is that the bond order BO ′i j between a pair of atoms can be obtained directly from the interatomic distance r ij as given in eq 2and plotted in Figure 1.Equation 2consists of three exponential terms:(1)the sigma bond (p bo,1and p bo,2)which is unity below ∼1.5Åbut negligible above ∼2.5Å;(2)the first pi bond (p bo,3and p b0,4)which is unity below ∼1.2Åand negligible above ∼1.75Å,and (3)the second pi bond (p bo.5and p bo,6)which is unity below ∼1.0Åand negligible above ∼1.4Å.
This leads to a carbon -carbon bond with a maximum bond order of 3.For carbon -hydrogen and hydrogen -hydrogen bonds,only the sigma-bond contribution is considered,leading to a maximum bond order of 1
The bond orders BO ′i j are corrected for overcoordination and for residual 1-3bond orders in valence angles using the scheme described in eqs 3a -f.While the 1-3bond order correction,described in eqs 3e and 3f,is applied to all the bonds in the molecule,the overcoordination correction (eqs 3b -d)is only applied to bonds containing two carbon atoms.The final bond orders in the molecule are obtained by m
ultiplying the bond orders from Equation 2by the correction factors from eq 3.Figure 2shows the effects of eqs 3a -f on the bond orders in an ethane molecule in which the C -C bond length is reduced from its equilibrium values of about 1.53Åto 1.0Å.This creates overcoordination on both the carbon and the hydrogen atoms,as the sum of bond orders around all atoms exceeds their valences (4for carbon and 1for hydrogen).As Figure (2)shows,application of Equations (3a -f)removes all of the 1-3bond orders,correcting the overcoordination on the hydrogen atoms,and,in addition,corrects most of the additional overcoordination
TABLE 1:General Parameters
parameter value description
equation λ150.0overcoordination bond order correction 3c λ215.61overcoordination bond order correction 3d λ3  5.021-3bond order correction 3e,f λ418.321-3bond order correction 3e,f λ58.321-3bond order correction 3e,f λ6-8.90overcoordination energy 6λ7  1.94undercoordination energy 7a λ8-3.47undercoordination energy 7a λ9  5.79undercoordination energy 7b λ1012.38undercoordination energy 7b λ11  1.49valence angle energy 8b λ12  1.28valence angle energy 8b λ13  6.30valence angle energy 8c λ14  2.72valence angle energy 8c λ1533.87valence angle energy 8c λ16  6.70valence
angle energy 8d λ17  1.06valence angle energy 8d λ18  2.04valence angle energy 8d λ1936.0penalty energy 9a λ207.98penalty energy 9a λ210.40penalty energy 9b λ22  4.00penalty energy 9b λ23  3.17torsion energy 10b λ2410.00torsion energy 10c λ250.90torsion energy 10c λ26-1.14conjugation energy 11a λ27  2.17conjugation energy 11b λ28
1.69
van der Waals energy
12b
E system )E bond +E over +E under +E val +E pen +E tors +
E conj +E vdWaals +E Coulomb (1)TABLE 2:Atom Parameters As Used in Equations 2,6,7,12,13,and 14a
bond radii
under/over coordination
Coulomb parameters heat increments units r o År o,πÅr o,ππÅp over kcal/mol
p under kcal/mol
ηEV  EV γÅI kcal/mol C    1.399  1.266  1.236
52.229.4
7.41  4.120.69218.6H
0.656
-
117.5
9.14
2.26
0.37
54.3
a
r o (ij ))1/2[r o (i )+r o (j
)].
Figure 1.Interatomic distance dependency of the carbon -carbon bond order.
BO ′ij )exp [p bo,1‚(r ij
r o
)p bo,2
]+exp [p bo,3‚
(r ij
πr o
)p bo,4
]
+
exp [p bo,5‚
(r ij
ππr o
)p bo,6
]
(2)
ReaxFF:A Reactive Force Field for Hydrocarbons
J.Phys.Chem.A,Vol.105,No.41,20019397
on the carbon atoms.
Val i in eqs 3a -3f is the valency of atom i (Val i )4for carbon,
Val i )1for hydrogen).∆′i is the degree of deviation of the sum of the uncorrected bond orders around an atomic center
from its valency Val i ,as described in eq 4.
Equation 5is used to calculate the bond energies from the corrected bond order BO ij .
2.2.Atom Under-/Overcoordination.From the valence
theory of bonding we know that the total bond order of C should not exceed 4and that of H should not exceed 1,except in hypervalent cases.However,as Figure 2shows,even after correction of the original bond orders BO ′i j a degree of overcoordination may remain in the molecule.To handle this we have added an overcoordination penalty term to the force field.
2.2.1.O V er-Coordination.For an overcoordinated atom (∆i >0),eq 6imposes an energy penalty on the system.The form of eq 6,ensures that E over will quickly vanish to zero for under-coordinated systems (∆i <0).∆i is calculated using eq 4,using the corrected bond orders BO ij from eq 3instead of the uncorrected bond orders from eq 2.
2.2.2.Under-Coordination.For an under-coordinated atom
(∆i <0),we want to take into account the energy contribution for the resonance of the π-electron between attached under-coordinated atomic centers.This is done by eqs 7a,b where E under is only im
portant if the bonds between under-coordinated atom i and its under-coordinated neighbors j partly have π-bond character (BO ij,π>0as calculated from the last two terms of eq 2).
2.3.Valence Angle Terms.Just as for bond terms,it is important that the energy contribution from valence angle terms goes to zero as the bond orders in the valence angle goes to zero.Equations 8a -d are used to calculate the valence angle energy contribution.We use the bond-order-dependent form in eq 8a to calculate energy associated with deviations in valence angle Θijk from its equilibrium value Θo .The f 7(BO )term as described in eq 8b ensures that the valence angle energy contribution disappears smoothly during bond dissociation.Equation 8c deals with the effects of over/undercoordination in central atom j on the valence angle energy.The equilibrium angle Θo for Θijk depends on the sum of π-bond orders (SBO)around the central atom j as described in eq 8d.Thus,the
TABLE 3:van der Waals Parameters Used in Equation 12a
atom units r vdW Å kcal/mol R γw ÅC    3.9120.086210.71  1.41H
3.649
0.0194
10.06
5.36
a
Arithmetic combination rules are used for all van der Waals parameters.
TABLE 4:Bond Parameters (D e in kcal/mol)As Used in Equations 2and 3
bond D e p be,1p be,2p bo,1p bo,2p bo,3p bo,4p bo,5p bo,6C -C 145.20.3180.65-0.097  6.38-0.26
9.37
-0.391
16.87
C -H 183.8-0.45412.80-0.0137.65H -H
168.4
-
0.310
10.25
-0.016
5.98
Figure 2.(a)Effect of the bond order correction in eq 2on the C -C and C -H bond orders in an ethane molecule in which the C -C bond is shortened to 1.0Åwith the rest of the geometry fixed.(b)Effects of shortening of the C -C bond length in ethane to 1.0Åon the relaxed C -H bond lengths as calculated by DFT and ReaxFF.Equilibrium C -C and C -H bond lengths are in italics and brackets.
BO ij )BO ′ij ‚f 1(∆′i ,∆′j )‚f 4(∆′i ,BO ′ij )‚f 5(∆′j ,BO ′ij )(3a)
f 1(∆i ,∆j ))12‚(
Val i +f 2(∆′i ,∆′j )Val i +f 2(∆′i ,∆′j )+f 3(∆′i ,∆′j )
+
Val j +f 2(∆′i ,∆′j )
Val j +f 2(∆′i ,∆′j )+f 3(∆′i ,∆′j )
)
(3b)
f 2(∆′i ,∆′j ))exp(-λ1‚∆′i )+exp(-λ1‚∆′j )(3c)
f 3(∆′i ,∆′j ))
1λ2‚ln {1
2
‚[exp(-λ2‚∆′i )+exp(-λ2‚∆′j )]}
(3d)f 4(∆′i ,BO ′ij ))
1
1+exp(-λ3‚(λ4‚BO ′ij ‚BO ′ij -∆′i )+λ5)
(3e)f 5(∆′j ,BO ′ij ))
1
1+exp(-λ3‚(λ4‚BO ′ij ‚BO ′ij -∆′i )+λ5)
(3f)
∆′i )
∑j )1
n bond
BO ′ij -Val i
(4)
E bond )-D e ‚BO ij ‚exp[p be,1(1-BO ij p be,1
)]
(5)
E over )p over ‚∆i ‚
(
1
1+exp(λ6‚∆i )
)
(6)
E under )-p under ‚
1-exp(λ7‚∆i )1+exp(-λ8‚∆i )
‚f 6(BO ij ,π,∆j )
(7a)
f 6(BO ij ,π,∆j ))
1
1+λ9‚exp(λ10‚
∑j )1
neighbors(i )
∆j ‚BO ij ,π)
(7b)
9398J.Phys.Chem.A,Vol.105,No.41,2001
van Duin et al.
equilibrium angle changes from around 109.47°for sp 3hybrid-ization (π-bond )0)to 120°for sp 2(π-bond )1)to 180°for sp (π-bond )2)based on the geometry of the central atom j and its neighbors.In addition to including the effects of π-bonds on the central atom j ,eq 8d also takes into account the effects of over-and under-coordination in central atom j (∆j )on the equilibrium valency angle,including the influence of a lone electron pair.The functional form of eq 8d is designed to avoid singularities when SBO )0and SBO )2.The angles in eqs 8a -d are in radians.
To reproduce the stability of systems with two double bonds
sharing an atom in a valency angle,like allene,an additional energy penalty,as described in eqs 9a and 9b,is imposed for
such systems.Equation 9b deals with the effects of over/undercoordination in central atom j on the penalty energy.
2.5.Torsion Angles.Just as with angle terms we need to ensure that dependence of the energy of torsion angle ωijkl accounts properly for BO f 0and for BO greater than 1.This is done by eqs 10a -c.The V 2-cosine term in eq 10a depends on the bond order of the central bond BO jk .In torsion angles with a central double bond (BO jk )2)the V 2-term is at its maximum (about 30kcal/mol,see Table 6).If BO jk deviates from 2the magnitude of the V 2-term rapidly diminishes.The valence-angle-dependent term sin(Θijk )‚sin(Θjkl )in eq 10a ensures that the torsion energy contribution disappears when either of the two valence angles (Θijk or Θjkl )approaches π.
To avoid excessive torsion contributions in systems containing two attached over-coordinated sp 3-carbon atoms,like an ethane molecule in which the central C -C bond length is reduced from its equilibrium value of about 1.52Åto 1.35Å,we include eq 10c,which reduces the influence of BO jk on th
e V 2-term in eq 10a when atoms j and k are over-coordinated [∆j >0,∆k >0].Equation 10b describes the smooth disappearance of the torsion energy contribution when one of the bonds in the torsion angle dissociates.
2.6.Conjugated Systems.Equations 11a and 11b describe the contribution of conjugation effects to the molecular energy.A maximum contribution of conjugation energy is obtained when successive bonds have bond-order values of 1.5as in benzene and other aromatics.
2.7.Nonbonded van der Waals Interactions.In addition to valence interactions which depend on overlap,there are repulsive interactions at short interatomic distances due to Pauli principle orthogonalization and attraction energies at long distances due to dispersion.These interactions,comprised of
TABLE 5:Valence Angle Parameters As Used in Equations 8a -d
valence angle units Θo,o degree k a
kcal/mol k b
(1/radian)2
p v,1p v,2C -C -C 71.31a 35.4  1.370.01
0.77
C -C -H 71.5629.65  5.29H -C -H 69.9417.37  1.00C -H -C 028.5  6.00H -H -C 00  6.00H -H -H
27.9
6.00
a
This value leads to an equilibrium angle of 180-71.31)108.69°for the single-bond C -C -C valence angle (eq 8d).
TABLE 6:Torsion and Conjugation Parameters (V 2and V 3in kcal/mol)As Used in Equations 10a -c
torsion angle a V 2V 3p t C -C -C -C 21.70.00-2.42C -C -C -H 30.50.58-2.84H -C -C -H
26.5
0.37
-2.33
a
Torsion angles not defined in this table (i.e.,C -H -C -C)are assigned torsion barrier energies of 0kcal/mol.
E val )f 7(BO ij )‚f 7(BO jk )‚f 8(∆j )‚
{k a -k a exp[-k b (Θo -Θijk )2
]}(8a)
f 7(BO ij ))1-exp(-λ11‚BO ij λ12
)
(8b)
f 8(∆j ))
2+exp(-λ13‚∆j )
1+exp(-λ13‚∆j )+exp(p V ,1‚∆j )
‚[
λ14-(λ14-1)‚2+exp(λ15‚∆j )
1+exp(λ15‚∆j )+exp(-p V ,2‚∆j )
]
(8c)
SBO )∆j -2‚{1-exp [-5‚(1
2
∆j
)λ16
]}
+
∑
n )1
neighbors(j )
BO jn ,π
∆j ,2)∆j if ∆j <0∆j ,2)0if ∆j g 0SBO2)0if SBO e 0(8d)
SBO2)SBO λ17if 0<SBO <1SBO2)2-(2-SBO)
λ17
if 1<SBO <2
SBO2)2if SBO >2
Θ0)π-Θ0,0‚{1-exp[-λ18‚(2-SBO2)]}E pen )λ19‚f 9(∆j )‚exp[-λ20‚(BO ij -2)2]‚
exp[-λ20‚(BO jk -2)2](9a)
f 9(∆j ))
2+exp(-λ21‚∆j )
1+exp(-λ21‚∆j )+exp(λ22‚∆j )
(9b)
E tors )f 10(BO ij ,BO jk ,BO kl )‚sin Θijk ‚sin Θjkl
[1
2
V 2
‚exp {p l
(BO jk
-
3+f 11(∆j ,∆k ))2}‚(1-cos 2ωijkl )+
1
2V 3
‚(1+cos 3ωijkl )]
(10a)f 10(BO ij ,BO jk ,BO kl ))[1-exp(-λ23‚BO ij )]‚
[1-exp(-λ23‚BO jk )]‚[1-exp(-λ23‚BO kl )](10b)f 11(∆j ,∆k ))
2+exp[-λ24‚(∆j +∆k )]
1+exp[-λ24‚(∆j +∆k )]+exp[λ25‚(∆j +∆k )]
(10c)
E conj )f 12(BO ij ,BO jk ,BO kl )‚λ26‚
[1+(cos 2ωijkl -1)‚sin Θijk ‚sin Θjkl ](11a)
f 12(BO ij ,BO jk ,BO kl ))exp [
-λ27‚(BO ij -112
)2
]
‚exp [-λ27‚(BO ij -112)2]
‚exp [
-λ27‚(BO kl -112
)2
]
(11b)ReaxFF:A Reactive Force Field for Hydrocarbons
J.Phys.Chem.A,Vol.105,No.41,20019399
van der Waals and Coulomb forces,are included for all atom pairs,thus avoiding awkward alterations in the energy descrip-tion during bond dissociation.In this respect,ReaxFF is similar in spirit to the central valence force fields that were used earlier in vibrational spectoscropy.To account for the van der Waals interactions we use a distance-corrected Morse-potential (eqs 12a,b).By including a shielded interaction (eq 12b),excessively high repulsions between bonded atoms (1-2interactions)and atoms sharing a valence angle (1-3interactions)are avoided.Figure 3shows how the bond energies,derived from eq 5,combine with the van der Waals interactions for diatomic C -C,C -H,and H -H systems to give a dissociation energy curve.
2.8.Coulomb Interactions.As with the van der Waals
interactions,Coulomb interactions are taken into account between all atom pairs.To adjust for orbital overlap between atoms at close distances a shielded Coulomb potential is used.
Atomic charges are calculated using the Electron Equilibration
Method (EEM)approach.20-21The EEM charge derivation method is similar to the Q Eq scheme;6the only differences,apart from parameter definitions,are that EEM does not use an iterative scheme for hydrogen charges (as in Q Eq)and that Q Eq uses a more rigorous Slater orbital approach to account f
or charge overlap.However,the γij in eq 13can be optimized to reproduce the Q Eq orbital overlap correction.The initial values for the EEM-parameters ( ,η,and γ,Table 2)were taken from
Njo et al.,22but these parameters were allowed to change during the FF optimization procedure.Intraatomic contributions of atomic charges,to account for the energy required to polarize the atoms,are taken into account in the energy scheme,23thus allowing a straightforward expansion of the force field for ionic compounds.
With the EEM-parameter values from Table 2,ReaxFF assigns a charge of -0.113to the carbon atoms in cyclohexane and charges of +0.050and +0.063to the equatorial and axial hydrogens,respectively.A Mulliken distribution,from a DFT calculation with the 6-31G**-basis set,results in charges of -0.174,+0.0859,and 0.0876to the cyclohexane carbon,equatorial and axial hydrogens while the Q Eq method 6gets a -0.24,+0.104,+0.137charge distribution.
2.9.Force Field Optimization Procedure.The FF was optimized using a successive one-parameter search technique as described by van Duin et al.24In general,our aim was to reproduce heats of formation to within 4.0kcal/mol,bond lengths to within 0.01Å,and bond angles to within 2°of their literature values.To use the QC data in the FF optimization procedure,structures relating to these data
were added to the FF training set.All molecules used in the heat of formation and geometry data comparisons were continuously energy minimized during the FF optimization while the structures relating to the QC data were kept fixed or were optimized with appropriate bond length or torsion angle restraints.
3.Results
3.1.Non-Conjugated Systems.3.1.1.Energy and Geometry Reproduction.Figure 4a and Table 7show how well the ReaxFF reproduces the heat of formation for nonconjugated closed shell molecules.The heat of formation predictions for ReaxFF are compared with those of the semiempirical
MOPAC-method,
Figure 3.Interatomic distance dependency of the carbon -carbon,carbon -hydrogen,and hydrogen -hydrogen bond-and van der Waals-energy terms of diatomic C -C,C -H,and H -H systems.Energy contributions from Coulomb interactions and energy effects related to under and overcoordination are ignored in the total energy curve.
E vdWaals )D ij ‚{exp [R ij ‚(
1-f 13(r ij )r vdW )]-2‚
exp
[1
2‚R ij ‚(
1-f 13(r ij )r vdW
)]}
(12a)f 13(r ij ))[r ij λ29
reactive web+
(1λw
)λ28]
1/λ28
(12b)
E Coulomb )C ‚
q i ‚q j
[r ij 3+(1/γij )3]
1/3
(13)
Figure 4.Comparison of calculated heats of formation with experi-mental data for nonconjugated (a),conjugated (b),and radical systems (c).
9400J.Phys.Chem.A,Vol.105,No.41,2001van Duin et al.

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。