Variable Charge Reactive Potential for Hydrocarbons to Simulate

Article pubs.acs.org/JPCA

Variable Charge Reactive Potential for Hydrocarbons to Simulate Organic-Copper Interactions Tao Liang,† Bryce Devine,‡ Simon R. Phillpot,† and Susan B. Sinnott*,† †

Department of Materials Science and Engineering, University of Florida, Gainesville, Florida 32611, United States U.S. Army Engineer Research and Development Center, Vicksburg, Mississippi 39180, United States

‡

S Supporting Information *

ABSTRACT: A variable charge reactive empirical potential for carbon-based materials, hydrocarbons, organometallics, and their interfaces is developed within the framework of charge optimized many-body (COMB) potentials. The resulting potential contains improved expressions for the bond order and self-energy, which gives a flexible, robust, and integrated treatment of different bond types in multicomponent and multifunctional systems. It furthermore captures the dissociation and formation of the chemical bonds and appropriately and dynamically determines the associated charge transfer, thus providing a powerful method to simulate the complex chemistry of many-atom systems in changing environments. The resulting COMB potential is used in a classical molecular dynamics simulation of the room temperature, low energy deposition of ethyl radicals on the Cu (111) surface (a system with ∼5000 atoms) to demonstrate its capabilities at describing organic-metal interactions in a dynamically changing environment.

I. INTRODUCTION Quantum computational methods have long been used to study systems of chemical and physical interest. They provide the highest materials fidelity relative to other currently available simulation methods. However, they are computationally expensive and cannot reach the time and length scales needed to investigate many problems of interest. Classical molecular simulations, including molecular dynamics (MD) and Monte Carlo (MC) algorithms, have been employed to examine time and length scales beyond the reach of quantum mechanical computations. The key component in classical molecular simulations is the interatomic empirical potential that describes the underlying interactions among the constituent atoms or molecules in gaseous, liquid, and solid states. A variety of empirical potentials have been successfully proposed for different types of bonding environments. For example, the Lennard-Jones (LJ) potentials1,2 are typically used for van der Waals (vdW) interactions, the embedded atom method (EAM) potentials3,4 for metallic bonding, and the Buckingham potentials5,6 for ionically bonded systems. The empirical potentials for covalently bonded systems can be classified into reactive and nonreactive types. Nonreactive potentials, conventionally named force fields in computational chemistry and biology, provide an accurate description of intraatomic or molecular interactions as long as the topology of strong interatomic covalent bonds remains unchanged throughout the simulations. Though they are unable to describe chemical reactions with bond breaking and/or new bond formation, they are widely used in the structural determination of macromolecules, such as proteins, by allowing the configurations to optimize their intramolecular interactions © 2012 American Chemical Society

from associated forces such as dihedral forces. Some widely used representative nonreactive force fields include AMBER,7 CHARMM,8 MM2,9 and GROMACS.10 Reactive potentials, which allow bond formation or dissociation in the simulations, have the capacity to describe the different bonding states of an atom or bond with the same parameter set and thus are sometimes referred to as bond-order potentials. Stillinger and Weber11 (SW) proposed simple reactive potentials for silicon that contain separated pairwise and many-body terms. Derived from chemical pseudopotential theory, the Abell model12 proposed a general charge-free model of bond energetics in which the many-body interactions are coupled with pairwise attractions by a bond-order term. Later, the bond-order term was formulated by Tersoff13−16 to model covalent bonds in silicon. In 1990 Brenner17 extended and parametrized the Abell−Tersoff potential to model carbon materials and hydrocarbon molecules. This potential is known as the Brenner, Tersoff−Brenner, or the reactive empirical bond order (REBO) potential. A second generation of REBO potential (REBO2)18 was proposed in 2002 to provide more accurate descriptions of short-range bonding in solid state carbon materials and hydrocarbons. The REBO2 potential has since been extended to include Si−C−H,19 C−H−O,20,21 and C−H−F22 interactions. It has been shown that the Tersoff or REBO formalism can be extended to metallic systems as they are mathematically equivalent to the EAM potentials.23 Most recently, Liang et al.24 parametrized REBO2 to model the Received: December 14, 2011 Revised: June 25, 2012 Published: June 27, 2012 7976

dx.doi.org/10.1021/jp212083t | J. Phys. Chem. A 2012, 116, 7976−7991

The Journal of Physical Chemistry A

Article

applied to various multicomponent materials, including the Al− Al2O3,39 Si−H−O,40 C−H−O,41 H−C−N−O,42 and Ni−C− H43−46 systems. Very recently, Watanabe47 advanced the concept of the ReaxFF potentials by developing a dynamic bond-order force field (DBO-FF) based on his extended SW potential.34 Compared with the ReaxFF potential, the bond order in DBO-FF is dynamically determined from its “bondorder element” term that reflects the valence electron configuration of each atom. Because the valence electron configuration is an inherent property of an element, the elemental parameters in the DBO-FF are proposed to be more general than those of ReaxFF, as exemplified in the use of the same parameter set for both Si/SiO2 and H2O. The variable charge scheme has also been linked to the Tersoff potentials for short-range interactions. For example, Yasukawa et al.48 developed such a potential for the Si−SiO2 system. In contrast to ReaxFF and DBO-FF, the nonelectrostatic, reactive, many-body effects in the Yasukawa potential are all dynamically coupled with the Tersoff potential’s pairwise attraction by a single bond order term. In addition, unlike the case in the ReaxFF and DBO-FF potentials, the strength of the short-range interactions in the Yasukawa potential is charge dependent, which was constructed on the basis of how the local valence of an ion varies according to its local environment. More importantly, all the parameters in the equations for the charge dependence of the short-range interactions and in the electrostatic energy terms in the Yasukawa potential were solely related to the elements being considered. The Yasukawa potential was refined by Yu et al.49 in the first generation charge optimized many body (COMB), hereafter referred to as COMB1, potential to reproduce the phase order in SiO249 and Cu,50 respectively. COMB1 subsequently evolved into two slightly different parametrizations of the secondgeneration COMB (COMB2) for different areas of study. The first set, labeled here as COMB2A, which included the influence of field effects on the self-energies and adjusted correction terms, was optimized for interfacial and bulk structures. COMB2A has been applied to Si−SiO2 systems, including amorphous SiO2,51 and the Cu−Cu2O52 and Hf− HfO253 systems. However, both COMB1 and COMB2A failed to predict the correct energetics of crystalline systems with low atomic or ionic coordination and some molecular systems. For example, they both predict O22− to be the most stable form of molecular oxygen, whereas first principles calculations correctly predict that O22− is unstable at any bond length. The second set, labeled COMB2B,54 was intended to correctly reproduce small molecules as well as solid-state and interfacial structures. A new electrostatic term that reflects the interactions between nucleus and charge is included in COMB2B to capture the correct dissociation behaviors of molecular anions. In additional to this new electrostatic term, the parameters inside the bond order of the short-range interactions in COMB2B were changed from central atom association to bond type association to increase the flexibility of the formalism. The COMB2B potential has been applied to the Cu−Cu2O54 system and to examine the oxidation of Cu. However, limited by the formalism for the bond order term, COMB2B has difficulty reproducing the energetics of O3 and molecular Cu/O systems. More importantly, the parametrization of COMBs 1 and 2 are focused on specific binary or ternary systems. As a result, their applications are usually limited to simple compounds or binary comprising compounds;

metallic, ionic and covalent bonding present in Mo−S systems. This paper also introduced an analytic function to explicitly capture the effects of atomic coordination on bond order. To model the nonreactive vdW force fields that are present in materials such as polymers, the short-range REBO force fields have typically been coupled with long-range LJ potentials using a cubic spline function, such that the LJ potentials are only active at distances greater than the REBO cutoffs. In 2001, Stuart et al.25 introduced a more sophisticated bond-order dependent switching function to provide a more realistic transition between LJ and REBO potentials. This method was termed the adaptive intermolecular REBO (AIREBO). Subsequently, the REBO + LJ potential family, including AIREBO, has been applied to model many different materials and processes where both short-ranged covalent and longranged van der Waals interactions are important, including fullerenes, carbon nanotubes,26 surface modification,27 polymerization,28 nanoindentation,29 and nanotriobology.30,31 LJ potentials have also been employed to model nonreactive interactions between organic materials and metals. For example, Nguyen et al.32 simulated the ion bombardment of biphenyl coated Cu surfaces using a combination of REBO potentials for hydrocarbons, EAM potentials for Cu, and LJ potentials for Cu−C and Cu−H interactions. However, LJ potentials are insufficient to characterize the strong reactive interactions between hydrocarbons and metals that are present in, for example, organometallics where metal−carbon interactions are a hybrid between covalent, metallic, and ionic bonding. Despite its successes, the REBO formalisms lack explicit electrostatic effects that are present in systems composed of elements with substantial differences in electronegativity. As a result, they are unable to directly examine the effects of charge and are challenged to describe chemical changes in a multicomponent system with dominant ionic features. In contrast, the overall effects of the electronic degrees of freedom have been implicitly included within fitted parameters in some potentials, such as the modified EAM (MEAM33) for metal oxides. Several different types of empirical potentials with fixed charge have been constructed for ionic materials, such as the Buckingham potentials6 for metal oxides, in addition to extended SW34 and Tersoff35 potentials for silica. Though fixed charge schemes work well for bulk-like systems, they are generally less robust in systems that are far from being bulk-like or that require charge distribution in response to changing system conditions. A variable charge scheme has been employed to empirically capture the local environment dependent Coulombic electrostatics, e.g., the electrostatics+ approach (ES+) proposed by Streitz and Mintmire36 for Al2O3. The charge distribution in the ES+ model is determined through charge equilibrium (Qeq) on the basis of the principle of electronegativity equalization (EE). The approach of the ES+ model is customizable to capture the nonelectrostatic effects in a specific system, e.g., ES + EAM for metal/metal oxides proposed by Zhou et al.37 Using an idea that is similar to Qeq for Coulombic interactions, van Duin et al.38 proposed the ReaxFF potentials or reactive force fields, where every nonelectrostatic interaction, such as bond angle, over- and undercoordination, and conjugation, is expressed in a separate and charge-independent energy term, each coupled to a specific bond order value. This feature enables the flexible capture of reactive and nonreactive force fields from different contributions by various potential energy terms. Hence the ReaxFF family of potentials has been 7977



Article

assumed to be included in the charge-dependent short-range interactions in eq 1, which are given as

the performance of the potential is less good when it is applied to multicomponent systems because the parameters are not very transferrable. In addition, the original Tersoff expression for the short-range bond energy, which is employed in COMBs 1 and 2, is reported18 to be inadequate to capture delocalized bonding in hydrocarbon systems. Moreover, the original Tersoff expression lacks the four-body dihedral term, which is also important for hydrocarbon materials such as polymer chains. In this paper we propose the third generation of COMB (COMB3) with the aim of providing a general, flexible, and robust empirical potential formalism that is capable of treating all different types of bonds within a multicomponent system in an integrated manner. The COMB3 formalism uses terms from COMB2B for electrostatic effects, REBO2 for short-range interactions, and the coordination function that was developed for the REBO potential for Mo−S systems. Several functions in COMB3 have been further refined to provide more flexibility and at the same time realistically capture the interactions within different systems. Though the work presented here is restricted to hydrocarbons and copper, the formalism and parametrization process should be applicable to any molecular or crystalline system. An important feature, as well as a significant challenge, for the COMB3 potentials is that all the element-specific parameters are kept the same for all bond types in metallic, organic, and ionic materials, which enables their simultaneous simulation within a given system. Although developed from different platforms, the ReaxFF potential family and COMB series have been built around the same two fundamental concepts of self-consistent charge equilibration and bond order. The difference in a large number of details regarding these two variable charge empirical potentials represents the different philosophies of the developers rather than the capabilities of the formalisms. The rest of this paper is organized as follows. Section II gives the general form of COMB3. A general guideline for the parametrization of COMB3 potentials is provided in section III. Section IVA is a summary of the fitted and predicted results for this potential which highlights its ability to describe material behavior in hydrocarbons and organometallics. Section IVB gives the results of MD simulations of ethyl radical deposition on the Cu (111) surface to demonstrate the strength of COMB3 in modeling organic-metal interactions.

U es[{q},{r }] = U self [{q},{r }] + U qq[{q},{r }] + U qZ[{q},{r }] + U polar[{q},{r }]

Because the electrostatic energies in the COMB family of potentials use the same concept as in the ES+ potential proposed by Streitz and Mintmire,36 here we only provide a brief summary of these terms in eq 2 to highlight the difference between COMB3 and COMB2B. The detailed origin of the electrostatic energies can be found elsewhere.36 Throughout this paper, the subscripts on the parameters stand for the element types and the number of subscripts on a parameter represents the dimensionality of the parameter, which is defined as the number of elements with which a parameter is associated. For instance, Pi means that the parameter P is only associated with the element type of the ith atom, whereas Pij is two-dimensional and is associated with the interaction between atoms i and j. In general, Pij is not necessarily equal to Pji. U self [{q},{r }] =

∑ [Viioniz(qi) + Vi field(qj ,rij)] i

(3)

Viionize(qi) = Ei0(0) + χi qi + Ji qi 2 + K iqi 3 + Liqi 4 + Lbarr (qi − qilim)qi 4

Vifield

1 = 2πε0

⎛ P χq PijJqj 2 ⎞ ij j ⎜ ∑ ⎜ 3 + 5 ⎟⎟ r rij ⎠ j > i ⎝ ij

(4)

NN

(5)

As specified in eqs 3−5, the energy to form a charge on an atom contains the ionization or affinity energy of an isolated atom (Vioniz(qi)) and a correction function, termed the field effect (Vfield(qi,rij)), which reflects the change of electronegativity and atomic hardness of the atom within its environment; qi and rij represent the charge on the ith atom and the distance between atoms i and j, respectively. Using the notation of Mortier et al.,55 the self-energy of an isolated ion can be approximated as a Taylor expansion series with respected to its charge (eq 4), where Ei0(0) is the energy of the neutral reference state of an element, which is set to zero for all species in the COMB3 potential. The last term in eq 4 acts as a barrier to prevent the charge of an atom from going beyond its charge limits, where Lbarr is zero when qi is within its charge limits and 100 eV/e4 when the charge is outside its charge limits, where e is the charge unit of an electron. The element specified parameter qlim has two values, the lower and upper charge limits, which are normally obtained from the valence electron configuration of an element, such as −4 and +4 for carbon, or assigned to correspond to a reasonable range, such as −2 and +2 for copper. The default values for these charge limits are the parameters QL or QU, as discussed in section IIB. In eq 4, the parameter χ is identified with the electronegativity, whereas J is associated with the chemical hardness or self-Coulomb interaction, which is inherent to each element. As outlined by Toufar et al.,56 the electronegativity and chemical hardness of an atom vary when it is embedded in a molecule or a crystal. These changes are reflected in eq 5, where ε0 is the vacuum permittivity. Unlike in COMB2B, the first term in eq 5 is dependent on qj and the second term is dependent on qj2. The parameters Pχ and PJ in eq 5 are associated with specific

II. DETAILS OF THE THIRD-GENERATION COMB POTENTIAL A. General Form of COMB3. As shown in eq 1, the total potential energy (Utot[{q},{r}]) of a system in COMB3 is composed of the electrostatic energies (Ues[{q},{r}]), shortrange interactions (Ushort[{q},{r}]), van der Waals interactions (UvdW[{r}]), and correction terms (Ucorr[{r}]), where {q} and {r} represent the charge array and coordinate array of the system, respectively. U tot[{q},{r }] = U es[{q},{r }] + U short[{q},{r }] + U vdv[{r }] + U corr[{r }]

(2)

(1)

The electrostatic energies include the energy to form a charge on each atom (Uself[{q},{r}]), the charge−charge interactions (Uqq[{q},{r}]), the charge−nuclear interactions (UqZ[{q},{r}]), and the energies associated with atomic polarizability (Upolar[{q},{r}]). The repulsion between nuclei is not explicitly formulated because their overall effects are 7978



Article

atoms overlap. For consistency, the same damping that is used with the Coulombic interactions is applied here:

bonds in COMB3, rather than specific atoms as in COMB2B; this enables more flexibility in the description of the environmental dependence of field effects. To avoid the Coulombic catastrophe that is introduced by a point charge model, COMB3 makes use of the Streitz and Mintmire charge density distribution functions, in which the charge density of an atom is a function of its charge, the spatial location (r), and atomic position (ri): ρi (r; qi) = Ziδ(|r − ri|) + (qi − Zi)fi (|r − ri|)

(6)

f (|r − ri|) = ξi 3π −1 exp( −2ξi|r − ri|)

(7)

⎛ ⎡ x 2 xy xz ⎤⎞ ⎜ ⎢ ⎥⎟ 1 3 ⎢ ⎜ 2 1 − 2 yx y yz ⎥⎟ Tij = ⎥⎟ 4πε0 | rij⃗ |3 ⎜ | rij⃗ | ⎢ ⎜ ⎢ zx zy z 2 ⎥⎟ ⎣ ⎦⎠ ⎝ × [1 − e−2ξir (1 + 2ξjrij + 2ξj 2rij 2)]

where ri⃗ j is the distance between units i and j and x, y, and z are the Cartesian components of ri⃗ j. Thus, the additional energies that are associated with atomic polarizability are the dipole selfenergy, the dipole−charge interaction, and the dipole−dipole interactions, which correspond to the three terms on the righthand side of eq 16, respectively:

Here, Zi is an effective point core charge treated as a fitted parameter, δ is the Dirac delta function, f(r) is a function that captures the radial decay of the electron density of the S-type orbital,57 and ξ is an orbital exponent that controls the length scale associated with this decay. Thus, the charge−charge interactions can be calculated by the integration over electron densities between pairs of atoms through the Coulomb integral operator Jqq: U qq[{q},{r }] =

∑ ∑ qiJijqqqj i

∫ d3r1 ∫ d3r1

Jijqq = [ρi |ρj ] =

U polar[{q},{r}] =

∑ ∑ (qiJijqZ Zj + qjJ jiqZ Zi) i

j>i

(9)

N

(10)

L=

where JqZ ij is the charge−nuclear coupling operator, which is a shielded nuclear attraction integral as defined by JijqZ = [j|ρi ] − [ρj |ρi ] [j|ρi ] =

∫ d3r

The polarizability P represents the relative tendency of the charge distribution to be distorted in response to an external electric field, which is defined as the ratio of the induced dipole moment, μ⃗, of an atom to the electric field, E⃗ , that produces this dipole moment (eq 13). Using the fluctuating dipole model by Sterne et al.,58 the polarization vector is directly calculated from the electrostatic field generated by the atomic charges, E⃗ qi and the neighboring induced dipoles (eq 13):

∑

q

Ei⃗ =

1 4πε0

N

∑ qj j≠i

μTi =

+

∑ ∑ μi⃗ Tijμj⃗ i

j>i

N

∑ i=1

N

1 MQ qi̇ 2 − U tot[{q},{r }] − υ ∑ qi 2 i=1

∂U tot[{q},{r }] ∂qi

(18)

(19)

By definition, μTi is the local electrostatic chemical potential of the ith atom and μ̅Ti is the average value of μTi of the system; ηD is chosen to be 0.01 in COMB3. The equations of motion for charge evolve iteratively by the Verlet integration algorithm with a time step of 0.005 fs until a convergence criterion is satisfied, which here is that the energy difference between two immediate time steps is less than 10−4 eV. B. Charge-Dependent Short-Range Interactions. The total short-range energy of a system is the summation over the bond energy (Vbond), which is described with the pairwise repulsive term, VR(rij,qi,qj), and pairwise attractive term, VA(rij,qi,qj), that is coupled with a charge independent bond order term bij:

Tij ,μj⃗ ] (13)

∂Jijqq rij⃗ ∂r | rij⃗ |

1 miri̇2 + 2

−mqqï = μTi − μT̅ + ηDqi̇

N j = 1,j ≠ i

i

q

The first term in eq 17 is the kinetic energy associated with the physical movement of the ions. The second term is the kinetic energy of the charges with a fictitious charge mass (MQ) of 0.002 amu, and υ is an undetermined multiplier that enforces the constraint that charge in the system is conserved. The equations of motion for charge are determined from the Euler− Lagrange equations with a damping factor ηD:

(12)

q

∑ μi⃗ ·Ei⃗

(17)

ρi (r)

⃗ ⃗ μi⃗ = PE i ( r ⃗) = Pi[Ei +

∑ i=1

(11)

|r − ri|

2Pi

+

As is the case with previous generations of the COMB potential, the Coulombic integral in eqs 8 and 9 is solved via a direct Wolf summation59 with a universal Coulombic cutoff radius, which is set to 11 Å in COMB3. The partial charge on each atom is treated as a dynamic variable in the COMB family of potentials. To effectively reach the charge equilibration, an extended Lagrangrian approach is used to express the time dependence of the total energy of the system, as indicated in eq 17, in which mi, ri, and qi correspond to the mass, position, and charge of the ith atom, respectively.

where r1 and r2 indicate the centers of ρi(r) and ρj(r). Here, r12 is the distance between density distributions. The charge− nuclear interactions are calculated by U qZ[{q},{r }] =

μi⃗ 2

(16)

ρi (r1) ρj (r2) r12

∑ i

(8)

j>i

(15)

(14)

where Pi is the polarizability tensor and Tij is the damped dipole field tensor. Because the atomic polarizability is assumed to be isotropic, Pi thus reduces to a scalar value that is only associated with the elemental nature of the atom. The induced dipole field tensor, Tij, is employed as a damped function that diminishes as 7979


The Journal of Physical Chemistry A U short[{q},{r }] =

⎛ QO i a Bi = ⎜⎜1 − ΔQ i ⎝

∑ ∑ Vijbond i

=

Article

j>i

∑ ∑ {Fc(rij)[V R(rij ,qi ,qj) i

j>i

− bijV A(rij ,qi ,qj)]}

ΔQ i =

(20)

Here, F c (r ij ) is the cutoff function. Both V R and V A exponentially decay with interatomic distance rij, as is the case with the Tersoff potential. The charge dependence on the pairwise repulsion and attraction in eq 20 makes use of the Yasukawa potential, where a charge dependent correction function, Di(qi), is added to the exponential coefficient of the short-range energies to reflect the change in atomic radius with charge:

{

V R (qi ,qj ,rij) = Aij × exp −λijrij +

QO = i

{

×

Di(qi) = DUi + |b Di(Q U − qi)|n Di b Di = (DLi − DUi )1/ n Di /(Q U − Q L ) i

i

(23)

(25)

ln DUi − ln(DUi − DLi ) ln Q U − ln(Q U − Q L ) i

(26)

DU and DL are parameters that reflect the change in atomic radius with charge, and QU and QL are the atomic charges that correspond to the limits of the valence shell. In addition, the change in the bond order with charge is reflected in the charge dependent function, Bij*(qi,qj), which equals unity under charge neutral conditions, and zero at the lower (QL) or upper (QU) limits of charge: Bij*(qi ,qj) = (Bi*B*j )1/2

bij =

1 σ−π (bij + bjiσ − π ) + bijπ 2

bijπ = bijRC + bijDH

(36)

(37)

Here, bπij is used to account for nonlocal conjugation effects in organic materials. In particular, bRC ij is used to determine the radical character and bDH depends on the dihedral angle for ij carbon−carbon bonds. Because bσ−π and bσ−π have different ij ji values, the overall bσ−π is taken as their mean value and contains contributions from the local bond angles (gij(cos(θijk)) and the atomic coordination (Pij(Ni1,Ni2,Ni3,...)) for central atoms i and j, respectively:

(28)

|a Bi|1/ n Bi ΔQ i

(34)

and

(27)

Bi* = [a Bi − |b Bi(qi − Q O )|n Bi ]

}

In eq 34, the single attraction term in the Tersoff type potential has been expanded to a summation over three attraction terms in a manner that is similar to what is done in the REBO2 potential. However, this modification makes it difficult to determine the charge dependence of the attraction in eq 22 and apply the mixing rules in eq 23, both of which are constructed on the basis of a single attraction term. Taking the mixing rule (eq 23) as an example, BCH is the square root value of the product between BCC and BHH. However, there are now three BCC values. To generalize eqs 21−23, an additional set of pairwise parameters for the carbon−carbon bond is developed, specifically A*CC, B*CC, λ*CC, and α*CC; these are used for the mixing rule in eq 23 and are further used in eqs 33 and 34. This set of carbon parameters describes molecular hydrocarbons well but does not reproduce the elastic constants in graphite and diamond carbon structures. In a manner that is analogous to REBO2,18 the bond order term, bij, in eq 20 modifies the short-range attraction on the basis of the local environment of the bond between atoms i and j, which is written as a sum of terms:

(24)

i

n n exp( −αCC rij)] ∑ [BCC n=1

Equation 23 makes use of traditional mixing rules to provide the default values for AB bond parameters from AA and BB bond parameters. Although the mixing rules fail where elements A and B have substantially different electronegativities, they provide a good starting point for the parametrization of A−B short-range interactions. Di(qi) is specific to each element type. Specifically,

b Bi =

}

3

Bii Bjj αij = 0.5(αii + αjj)

i

1 * λCC[Di(qi) + Dj(qj)] 2

{

}

i

(32)

1 * A VCC = Bij*(qi ,qj) exp αCC[Di(qi) + Dj(qj)] 2

In eqs 21 and 22, parameters Aij, Bij, λij, and αij are determined for each bond type with the default values of

i

1 (Q − Q L ) i 2 Ui

(33)

(22)

n Di =

(31)

(21)

1 × exp −αijrij + [αiiDi(qi) + αjjDj(qj)] 2

λij = 0.5(λii + λjj)

1 (Q − Q L ) i 2 Ui

{

V A(qi ,qj ,rij) = Bij Bij *(qi ,qj)

Bij =

(30)

R VCC = A CC exp −λCCrij +

}

Aii A jj

⎟ ⎟ ⎠

where nBi is equal to 10 for all elements. As described in ref 18, a single attraction term in eq 22 is insufficient to describe changes in bond energies in carbon systems. Consequently, eqs 21 and 22 have been modified as follows for carbon−carbon bonds:

1 [λiiDi(qi) 2

+ λjjDj(qj)]

Aij =

n Bi ⎞−1

(29) 7980



Article

N

The function Fc(rij) in eqs 20, 38, and 43 is the Tersoff’s cutoff function, which smoothly terminates the short-range interactions with a region defined by two cutoff radii:

bijσ − π = {1 + [ ∑ [Fc(rik)Nikcrossζ(rij ,rik)gij(cos(θijk))] k ≠ ij

+ Pij(Ni1 ,Ni2 ,Ni3 ,...)]ηi }−1/2ηi

ζ(rij ,rik) = exp[βijmi (rij − rik)mi ]

⎧1 rij ≤ R ijmin ⎪ ⎪ ⎡ ⎛ r − R min ⎞⎤ ⎪1 ij ij ⎟⎥ R ijmin < rij ≤ R ijmax Fc(rij) = ⎨ ⎢1 + cos⎜⎜π max min ⎟⎥ R R − ⎪ 2 ⎢⎣ ⎝ ij ij ⎠⎦ ⎪ ⎪0 rij > R ijmax ⎩

(38) (39)

6

gij(x) =

∑ bijang_nx n n=0

min min where the default Rmin ij is the arithmetic mean of Rii and Rjj , max with an analogous expression for Rij . From a mathematical view, the Tersoff cutoff function provides a cosine function transition between 1 and 0 within two cutoff limits, which is also used in eqs 45 and 47. For convenience, only the two cutoff limits of the Tersoff cutoff functions in eqs 45 and 47 are specified. Given the inherit complexity of the different types of bonding in hydrocarbon systems, the bond-order term in COMB3 uses the idea proposed in ref 18 for REBO2 which has proven to be successful for describing hydrocarbons. In particular, the angular function for carbon−carbon interactions is given as

The parameter Ncross is the coordination number dependence ik by element type, which is defined below. The function ζ(rij,rik) is an asymmetric term because it returns a value different from unity if rij differs from rik. More specifically, this term weakens the longer bond in a way that is comparable to a screening function. The parameter βij is normally set equal to αij in eq 22 as a default but is modified as needed to control the strength of ζ(rij,rik) for cases such as carbon−carbon bonding when the asymmetry is complex. The angular function g(x) is a sixthorder polynomial as shown in eq 40, where the parameters bang_0 to bang_6 are two-dimensional and specific to the bond type between atoms i and j. The coordination function Pij(Ni1,Ni2,Ni3) is a tricubic spline function for the bond types in which the bond properties have a complex dependence on the coordination numbers of the central atom, such as is the case in carbon-based bonds. The atom type used to calculate the three variables (Ni1, Ni2, and Ni3) for this spline function are categorized according to the elemental electronegativity with respected to the carbon atom. In particular, the first dimension (Ni1) in this tricubic spline function is the number of neighbors of the carbon atom, NiC, the second dimension (Ni2) is the number of neighbors of the atoms with smaller electronegativity than carbon, such as NiH + NiCu, and the third dimension (Ni3) is the number of neighbors of the atoms with larger electronegativity than carbon, such as NiO + NiF. For the rest of the bond types, the coordination term employs the same analytic function that was introduced in ref 24. Pij = c0 Ωi + c1e c2Ωi + c3

* (cos(θ )) + Fc(Ωi)[γ (cos(θ )) gCC(cos(θ )) = gCC CC * (cos(θ ))] − gCC

(41)

(42)

The superscript represents the element identification number and Nil represents the total number of lth element neighbors of atom i, l

Nil =

∑ f ikC (rik)Nilcross k ≠ i ,j

(45)

where Fc(Ωi) is the Tersoff cutoff function on Ωi, in which two cutoff limits Ωmin and Ωmax are 2.2 and 2.7, respectively. In eq 45, gCC * (cos(θ)) is a combination of three sixth-order polynomials of cos(θ), which cover bond angles 0−109.45°, 109.45−120°, and 120−180°, respectively. An additional sixthorder polynomial, γCC, is used for Ωi < 2.2 or a coordination number i

⎡⎛ vdW ⎞12 σ vdW ⎢⎜ ij ⎟ 4εij ⎢⎜ ⎟ r ij ⎝ ⎠ j>i ⎣⎢

NN

=

∑∑ i

⎛ σ vdW ⎞6 ⎤ ij ⎟⎥ − ⎜⎜ ⎟⎥ ⎝ rij ⎠ ⎥⎦

(49)

Here, εvdW and σvdW reflect the strength and equilibrium ij ij distance of the vdW interactions, respectively. The vdW interactions are truncated and shifted to zero at the Coulombic cutoff radii. To avoid excessively high repulsion between shortrange bonded atoms, a cubic spline function is added from 0.95 σij to smoothly terminate the vdW interactions at Rmax ij , the upper cutoff radii for short-range bonds. For any binary vdW bond, the εvdW and σvdW are taken as the geometric and ij ij arithmetic means of their values in ii and jj vdW bonds, respectively. D. Correction Terms. The energy terms described in this section are primarily used to modify the energy contribution for some specified angles, and are therefore designated as correction terms. They consist of Legendre polynomials (LPs) up to sixth order and a bond bending term: N

U corr[{r }] =

N

6

∑ ∑ ∑ {∑ K ijkLP (cos(θijk)) n

i

j>i k≠i

n=1

θ 2 BB + K ijk [cos(θijk) − cos(K ijk )] }

(50)

III. PARAMETERIZATION OF COMB3 This section discusses the details of parametrization of the COMB3 potential. Despite the changes from prior versions of COMB, the COMB3 formalism is compatible with the previous

The LPs, whose detailed forms can be found in ref 60 are a set of symmetric energy penalties on the bond angles. In contrast, the bond bending term provides an asymmetric energy penalty that depends on a specific bond angle Kθijk. All the parameters in 7982



Article

versions of the COMB family of potentials. For copper−copper bond parameters, the COMB3 potential directly takes the parameters from COMB2B.54 The resulting properties for metallic copper from these parameters are the same as those in ref 54. However, the parameters in the previous versions of the COMB family of potentials, especially the one-dimensional parameters such as the copper one-dimensional parameters in COMB2B, are not readily transferred to multicomponent systems. Thus the one-dimensional parameters were redone in this work. As indicated in section II, the COMB3 potential has a complex formalism for the potential energies that has different physical origins. To ensure transferability, a wide range of pure and binary systems with differing local environments are considered in the parametrization process. The properties of these systems are obtained from published experimental data and/or the electronic structure calculations. For two-dimensional lattices or three-dimensional crystalline structures, the electronic structure calculations are performed with plane-wave DFT calculations using the software VASP (Vienna Ab initio Simulation Package)61 with the projector augmented wave (PAW)62,63 method within the generalized gradient approximation (GGA). For molecular systems, the electronic structure calculations are performed using the Gaussian0964 computational chemistry software package. Different theoretical levels of ab initio methods in Gaussian09 are employed for different purposes. Unless otherwise specified, the geometric optimization for molecular systems is carried out using B3LYP/6311G*,65,66 the nonequilibrium properties are calculated with CASSCF(M,N)/6-31+G(d),67−69 and the properties of excited states such as are used to determine ionization energies are calculated with the CCSD(t)/cc-pVTZ70,71 level of theory. To unify the fitting database, rather than using total energies of these systems, the parametrization focuses on reaction energies, i.e., the energy differences between products and reactants, where the reactants are normally chosen to be the most stable phases of the given species. Furthermore, these reaction energies are scaled with respect to available experimental values. The resulting databases from the theoretical calculations are thus all consistent with available experimental values, which are labeled as ‘target’ values in subsequent sections of this paper. Based on the origin of the potential energies, the overall parametrization process in COMB3 is performed in three steps: (1) fitting the parameters for the short-range interactions in pure systems; (2) fitting parameters for electrostatic energies and the charge dependence of short-range interactions; (3) fitting the parameters for the short-range interactions in binary systems. A. Parameterization of Pure Systems. Because the electrostatic energy terms for an uncharged system are zero, the COMB3 potential for such a system is reduced to a Tersoff type of potential with forms that are similar to forms presented in refs 18 and 24. Therefore, the parametrization of the COMB3 potential starts from fitting the pure systems and the charge on each atom is assumed to be zero. As is the case for other Tersoff-type potentials, parametrization in this step (Table 1) consists of two substeps: (1) fitting the pairwise parameters, Aii, Bii, λii, and αii, and (2) fitting the many-body parameters, bang_0 to bang_6 and c0 to c3. The carbon systems used for this stepwise parametrization include dimer and trimers, two-dimensional square and triangle lattices, and three-dimensional graphite, diamond, simple cubic

Table 1. Parameter Details for the Pure Systems parameters

descriptions

Aii, Bii, λii, αii ηi, mi βii bang_0 to bang_6 for ii bond

pairwise term on bond order asymmetric term bond angle term

c0 to c3 for ii bond vdW εvdW ii , σii BB KLP to KLP iii 1 iii 6, Kiii , and Kθiii

coordination term vdW interaction correction functions

default values

locations eqs 33 and 34 eqs 38 and 39 eq 39 eq 40

none 1 αii 0 except

eq 41 eq 49 eq 50

bang_0 =0.5 0 0, 4 0

(SC) and face centered cubic (FCC) crystals. The fitting database for carbon systems is listed in Table 2, where the data points for the structures with the superscript a are taken from refs 72 and 73. Table 2. Database for Carbon Systems (SC = Simple Cubic and FCC = Face Centered Cubic) bond length (Å)

energy per atom (eV/atom)

energy per bond (eV)

notes

dimer

1.31

−3.06

−6.12

used for the first substep

isolateral C3 square triangle graphitea diamonda SCa FCCa linear C3

1.37

−4.31

−4.31

1.62 1.71 1.42 1.54 1.93 2.18 1.32

−5.23 −3.84 −7.40 −7.36 −4.70 −2.86 −4.573

−2.64 −1.28 −4.93 −3.68 −1.57 −0.46 −6.860

120° 130° 140° 150° 160°

1.32 1.32 1.32 1.32 1.32

−4.547 −4.558 −4.564 −4.568 −4.571

−6.820 −6.836 −6.846 −6.853 −6.857

structures

a

C3 C3 C3 C3 C3

used for the second substep

References 72 and 73.

Table 2 indicates that the calculated bond length for graphite from ref 73 is exactly the same as the experimental value74 (1.42 Å), whereas the cohesive energy is 0.03 eV/atom smaller than the experimental value (−7.37 eV). In the fitting, the energy per atom for graphite is scaled down to the experimental value. The same scaling factor, 7.37/7.40, is applied to the energies for diamond, SC and FCC structures. As specified in eqs 33 and 34, there are two sets of pairwise parameters for carbon systems. Consequently, the first substep of this parametrization is performed twice for carbon. The first uses one attraction term, as shown eq 22. The second uses eqs 33 and 34, where the short-range attraction consists of three terms. The many-body parameters for carbon systems are obtained using the same fitting strategy as in the REBO2 vdW potential.18 The εvdW CC and σCC are fit to the experimental 75 interlayer binding energy (−0.05 eV/atom) and the lattice parameter c (6.71 Å) of graphite, respectively. After adding vdW interactions to the carbon systems, the energetics and the elastic properties of graphite are slightly altered, and the final angular functions gCC * (cos(θ)) and γCC(cos(θ)) in eq 45 are 7983



Article

charge-associated parameters of an atom. The involved parameters for this step are listed in Table 5.

adjusted to reproduce the elastic properties of carbon systems with vdW interactions. No three-dimensional parameters are used for pure carbon systems. The properties of graphite, diamond and graphene from experiments, REBO2, ReaxFF, and COMB3 are reported in Table 3. In general, the COMB3 properties of selected carbon

Table 5. Parameters for the Charge Dependent Energy Terms parameters

Table 3. COMB3 Properties of Carbon Systems exp/ DFT

properties Ec (eV/atom) a (Ǻ ) c (Ǻ ) bulk modulus (GPa) shear modulus (GPa) C11 (GPa) C12 (GPa) C13 (GPa) C33 (GPa) C44 (GPa) C66 (GPa) interlayer binding energy (eV/ atom) vac formation (eV) Ec (eV/atom) a (Ǻ ) bulk modulus (GPa) shear modulus (GPa) C11 (GPa) C12 (GPa) C44 (GPa) vac formation (eV) Ec (eV/atom) γ11 (N/m) γ12 (N/m) γ66 (N/m)

a

REBO2

ReaxFF

−7.41 2.46

−7.34 2.50 6.60

Graphite −7.37c 2.46c 6.71c 277c 214c 1060c,d 180c,d 15c,d 36.5c,d 4c,d 440c,d −0.05e 7.6a Diamond −7.34c 3.57c 451c,f 567c,f 1075c,f 139c,f 567c,f 7.2a Graphene −7.32e 358c,g 55c,g 152c,g

b

COMB3

−0.04

−7.41 2.46 6.74 279 218 1050 185 1 40 1 440 −0.05

7.6

7.6

−7.36 3.57

−7.32 3.61

1080 130 740 7.2 −7.37

Xi, Ji, Ki, Li Pχii, PJii ξi, Zi Pi DUi, DLi, QUi, QLi

−7.27 3.57 470 581 1080 158 661 7.2 −7.36 360 55 152

Reference 18. bReference 38. cReference 76. dReference 77 and 78. Reference 75. fReference 79. gReference 80.

e

polymorphs are in good agreement with experimental values. The predictions for C13, C33, and C44 in graphite have larger quantitative errors, arising from the influence of vdW interactions in the region between the upper cutoff radii Rmax for short-range interactions and σvdW CC for the vdW interactions. For hydrogen, the pairwise parameters are fit to the stretched potential energy surface (PES) of H2 and isolateral H3. The many-body parameters are fit to the bending PES of H3 from 120° to 180°. The gHH(cos(60°)) is fit to the energy of the isolateral H3. The resulting properties are given in Table 4. B. Parameterization of Charge-Dependent Terms. The second step of the parametrization is used to determine the

H 2(r ) + e− → H 2−(r )

bond length (Å)

properties

exp/target

COMB3

exp/target

COMB3

−2.358 −0.61

−2.259 −0.68

0.74 1.09

0.74 1.19

eq 4 eq 5 eqs 6−7 eq 13 eqs 24−32

defaulted value none 0 2, 0 0 −0.5, 0.5, 4, −4

(51)

E(H2−(r))

the reaction energy, E (r) = − E(H2(r)), consists of the self-ionization energy, the molecular polarization, charge− charge interactions, nuclear−charge interactions, and charge dependence on the short-range interactions, in which the selfionization energy, the molecular polarization and the shortrange interactions are determined with the previously fit parameters. The only unknown parameters that influence the energy for such a reaction are the above parameters for the rea

Table 4. Properties of Hydrogen

H2 isolateral H3

ionization energy field effects Coulomb interaction atomic polarizability charge dependence on short range

locations

The fitting process for this step is performed in three substeps: (1) fitting Xi, Ji, Ki, and Li to the electron affinity, first, second, and third order of ionization energies of an isolated atom; (2) fitting the atomic polarizability, Pi, to the molecular polarization that is perpendicular to the bond in a dimer system; and (3) fitting the rest of the parameters to the reaction energies of inserting or detaching an electron from the dimer systems. In the first substep, data points are obtained either from experiments81 or from CCSD(t)/cc-pVTZ calculations. Because the ionization energy is expressed as a Taylor expansion about its charge, the relative electronegativity of elements is primarily related to Xi, which is the energy derivative with respect to charge at a charge of 0. Given that charge in this model can be a fractional number, the electron affinity is adjusted to properly reproduce the relative magnitude of Xi values of the elements, i.e., the electronegativity trend of the elements. The original ES+ model is capable of reproducing the polarization along with the bond direction of a dimer and predicts zero polarization at directions that are perpendicular to the bond. Here, an isotropic atomic polarizability is added to account for the molecular polarization that is perpendicular to the bond direction. Thus the parameter Pi can be directly determined from the atomic polarization in a bonded dimer system. However, the isotropic atomic polarizability is unable to reproduce the asymmetric polarization that is observed in the C2 dimer. So for carbon, the parameter Pi is fit to the mean value of the C2 polarizations of its two perpendicular directions. The parameters that are involved in the last substep of this parametrization step are Pχii, PJii, ξi, Zi, DUi, DLi, QLi, and QUi, in which QLi and QUi represent the lower and upper charge limits of an element. Except for QL for hydrogen, the QLi and QUi are directly determined according to the valence electron configuration of the element, for example −4 and +4 for carbon, or assigned to a reasonable range, such as −2 and +2 for copper. In the case of the reduction reaction associated with adding an electron to the H2 dimer with interatomic distance r:

a

energy per atom (eV/atom)

descriptions

7984



Article

third substep parametrization. By varying the interatomic distance of H2, the resulting energy curve of H2− is the dissociation energy curve of H2−. The dissociation energy curves for H2− and Cu2− as well as H2+ and Cu2+ are calculated with B3LYP/6-311G* and for C2− and C2+ are calculated with CASSCF/6-31+G(d). The reaction energies as a function of distance are then calculated from these dissociation energy curves and used for the fitting. Due to the high energy of the antibonding orbital in H, the reaction in eq 51 is unlikely to occur and hydrogen has not been observed as H− in any chemical reactions. However, the charge in the COMB3 potential is treated as a variable, which allows it be presented as a fractional value. As for the reaction in eq 51, H−0.5−H−0.5, which corresponds to H2−, has a lower energy than H2 because a single H atom has an electron affinity of 0.75 eV. To account for the antibonding orbital in hydrogen, QL for hydrogen is preset to −0.5 and subjected to the parametrization steps. Figure 1 gives the reaction energies as a function of interatomic distance of detaching (H2+) or inserting (H2−) an electron from or to H2. Figure 2. Reaction energies of detaching or inserting an electron from or to (a) C2 and (b) Cu2.

states among the C2 orbitals, which cannot be reproduced with an empirical potential that lacks the explicit treatment of electrons. Thus the potential is fit to the global energy minimum of C2− and the overall trend of the reaction energies for adding an electron to C2. In general, the COMB3 potential predicts the correct trends for both positively and negatively ionized carbon and copper dimers with a wide range of interatomic distances. After these charge-associated parameters are determined, the bond angle parameters bang_0 to bang_6 from the first step fitting might be slightly adjusted to account for the charge distribution in the asymmetric pure systems. C. Parameterization of the Binary Systems. The parameters for the last step of the parametrization are listed in Table 6. Compared with COMB2B, four more parameters Pχij, PJij, Pχji, and PJji (eq 5) are introduced to adjust the charge of the atom in the binary systems. Similar to the parametrization process in the first step, the potential starts the third step with fitting the pairwise terms Aij, Bij, λij, αij, Pχij, PJij, Pχji, and PJji to the phase orders and charges of a variety of binary systems, which is

Figure 1. Reaction energies of detaching or inserting an electron from or to H2.

As interatomic distance increases, the reaction energies for H2+ decrease. In theory, the energy of detaching an electron from H2 at infinite distance is equal to the first ionization energy of a single H atom because H2 is no different than two single isolated H atoms at that distance. However, because the charge is allowed to be a fraction in COMB3, the reaction energies for H2+ in the COMB3 approach correspond to the smaller value of 2 × H0.5 and H+. From the mathematical view of the self-ionization expression used (eq 4), two half ionized atoms normally have a smaller total ionization energy than the first ionization energy of an atom. As a result, the predicted charge in COMB3 from the electrostatic parameters obtained from this methodology normally range from 50% to 80% of the Mulliken charge from Gaussian calculations or the Bader charge from VASP calculations. Instead of reproducing the charge from higher level calculation methods, the COMB3 potential is thus focused on the electronegativity trend for the elements and the charge trend in multiple systems. Within the bonding distance, the predicted reaction energies for both H2 anions and cations are in good agreement with the high fidelity electronicstructure calculation methods considered. To obtain both the positive and negative ionization curves of H2 correctly, it is important to accurately model chemical reactions that involve hydrogen atoms. Parts a and b of Figure 2 illustrate the positive and negative ionization energies of C2 and Cu2, respectively. In the case of carbon anions, the multiple minima (Figure 2a) from the CASSCF calculations correspond to different electron partition

Table 6. Parameters for the Binary Systems parameters

descriptions

Pχij, PJij, Pχji, PJji Aij, Bij, λij, αij

field effects pairwise short range

βij , βji bang_0 to bang_6 for ij and ji bonds c0 to c3 for ij and ji bonds PCC and PCH

on symmetric term bond angle term

eq 5 eqs 33−34 eq 39 eq 40

coordination term

eq 41

coordination term for hydrocarbons radical spline spline on torsion angle correction functions

eq 46 eq 48 eq 50

πCC φCC LP BB KLP ijk 1 to Kijk 6, Kijk , and Kθijk 7985

defaulted values

locations

Pχii, PJii, Pχjj, PJjj eq 23 αij, αij that of ii and jj bonds that of ii and jj bonds 0 0 0 0



Article

Table 7. Heats of Formation ΔH (kcal/mol) from Experiments, REBO2, ReaxFF, and COMB3a

a

species

expb

Δ(REBO2)c

CH2 CH3 methane C2H acetylene ethylene H3C2H2 ethane cyclopropene CH2CCH2 propyne cyclopropane propene n-C3H7 i-C3H7 propane cyclobutene 1,3-butadiene CH3CHCCH2 1-butyne 2-butyne cyclobutane 1-butene cis-butene i-C4H9 t-C4H9 n-butane isobutane 1,3-pentadiene 1,4-pentadiene cyclopentene 1,2-pentadiene 2,3-pentadiene cyclopentane 2-pentene 1-butene, 2-methyl 2-butene, 2-methyl n-pentane isopentane neopentane benzene cyclohexane naphthalene root-mean-square error

93.2 35.8 −15.99 135 54.33 14.52 28.3 −16.52 68.3 47.4 45.8 16.8 8.3 25 22.2 −19.7 41.5 29.8 42 42.8 38 12.7 4.9 2.3 19.4 15.2 −23.5 −25.4 22.8 29.1 13.9 38.4 36.5 −10.6 −1.1 −2 −3.2 −27.5 −28.81 −31.3 24 −20 31.04

−3.34 −0.32 −0.71 −1.58 −0.37 −0.86 −5.53 −0.79 30.41 5.65 5.67 0.4 −0.69 −4.76 −5.55 −0.12 27.54 20.25 4.98 6.15 10.97 −0.28 0.18 −0.75 −1.62 −5.09 1.16 3.1 21.15 0.89 2.16 −1.03 −2.7 −6.49 0.12 0.56 −1.59 2.64 4 6.78 −2.25 −1.29 6.47 8.50

Δ(ReaxFF)d 1.11 3.55 −19.95 8.56 −3.80 −3.13 1.69 22.12 −6.17 2.23 6.12 2.64 −4.02 −4.55 5.04 −3.64

−3.91

−2.19 −4.38 0.31 1.45 −2.78 −2.21

0.86

−0.73 2.07 −1.16 −1.53 5.14

Δ(COMB3) 0.96 0.78 0.77 1.05 0.81 3.72 −1.04 −1.21 5.14 2.57 −7.53 −2.91 0.82 −0.12 −1.92 0.47 −1.76 11.79 −2.30 −8.13 −22.34 1.11 0.62 −2.87 −4.80 4.34 2.09 0.01 8.42 2.12 −6.29 −0.83 −7.11 −2.55 −8.61 −0.09 −7.53 3.92 10.31 7.43 0.45 4.22 4.48 5.72

The values from the empirical potentials are reported as the difference from the experimental value. bReference 83. cReference 18. dReference 38.

followed by fitting the many-body parameters to the properties of the binary systems. For hydrocarbons, the stretching PES of CH and CH2 is used to determine the pairwise short-range terms; the parameters Pχij, PJij, Pχji, and PJji use their default values. The parameters on the bond angle term are fit to the bond bending PES of CH2. As specified in section IIB, the coordination functions for hydrocarbons, PCC and PCH, are expressed as two-dimensional spline functions. In addition, two three-dimensional spline functions (πCC in eq 47 and φCC in eq 48) are used to determine the energy dependence on the radical effects and torsion angle of the bonded C atoms. The fitting species and procedure for these spline functions are the same as for

REBO2. The resulting heats of formation for hydrocarbons are reported in Table 7 in section IV. Lastly, the potential is fit to the energy barriers for the reactions given in eqs 52 and 53, where the reaction in eq 52 is for H to bridge two methyl groups and the reaction in eq 53 is that of a H atom bonding to a H atom and a methyl radical. The energetics and geometries of the transition states for these two reactions made use of the values reported in Figures 11 and 12 in ref 38. In particular, KθH−C−C and KBB H−C−C are used to get the energy of the transition state of the reaction in eq 52 and θ BB KθH−C−H, KBB H−C−H, KH−H−C and KH−H−C are fit to the energy of the transition state of the reaction in eq 53. •

7986

•

H3C−H−•CH3

H3C + H−CH3 ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ H3C−H + •CH3

(52)


The Journal of Physical Chemistry A •

H−H−•CH3

H + H−CH3 ⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯⎯→ H 2 + •CH3

Article

body interactions into the charge-dependent, pairwise shortrange terms via the bond order term, where the parameters for many-body interactions have a reduced definition on the parameter dimension. For example, the parameters in the threebody bond angle term in COMB3 are specified to be dependent on the bond type between atoms i and j, which significantly reduces the number of parameters. As a result, the COMB3 potential only includes representative structures in its training database, which is substantially less than the number of structures used in the training database for the ReaxFF potentials. The excluded structures, as well as the major transition states in most cases, can be properly predicted by COMB3 with one-dimensional and two-dimensional parameters and proper bond order values, as indicated in sections IIID and IV. In a few cases, the COMB3 potential is less able to capture the energetics and geometries of selected transition states, for example the transition state in a bridging methylene between two methyl groups (CH3−CH2−CH3). At the same time, the charge dependence on short-range interactions enables the COMB3 potential to properly capture the energetics of charged molecules as indicated in Figures 1 and 2. Despite these differences, both potentials contain variable charge equilibration schemes and bond order terms that allow the connectivity for the chemical bonds in the real materials to be changed and dynamically determined according to their local environment. With regard to modeling organic-metal interactions, the ReaxFF potentials demonstrated their capabilities by accurately predicting several surface reactions associated with the deposition of hydrocarbons on Ni surfaces.45,46 The analogous capability of the COMB3 potential is demonstrated in section IV, where the results of deposition of selected hydrocarbon molecules on Cu surfaces are discussed. Though the systems are slightly different, some similar chemical phenomena, such as hydrocarbon dissociation and the formation of new carbon− metal bonds, are predicted by both potentials. D. Summary. This section provides a step by step parametrization scheme for the COMB3 potentials. The associated parameter sets for C−H−Cu systems are provided as the Supporting Information for this paper and will be available within the LAMMPS82 (Large-scale Atomic/Molecular Massively Parallel Simulator) MD simulation software. This systematic parametrization scheme has related each energy term in COMB3 to certain properties of specified systems and provided a guideline for the further parametrization of COMB3 potentials. The predicted results have been compared with extensive sets of data from experiments, literature and/or quantum mechanically based calculations. This comparison indicates that the COMB3 formalism is capable of producing the correct energetics, geometries, and elastic properties as well as charge trends for carbon-based materials, hydrocarbons, copper, and multicomponent compounds involving these elements.

(53)

With these parameters, COMB3 is capable of reproducing the reaction barriers for the reactions in eqs 52 and 53 from DFT calculations, which are 0.6 and 2.3 eV, respectively. In the case of a bridging methylene between two methyl groups as illustrated in Figure 10 in ref 38, the transition state occurs at an interatomic carbon distance from methylene and one side of methyl is 1.52 Å and the other side of methyl is 2.04 Å, which is outside the cutoff radius (2.0 Å) for both REBO2 and COMB3. As a result, the energy barrier for methylene to bridge two methyl groups is 3.2 eV for REBO2 and 2.1 eV for COMB3, which is a significantly larger than that predicted by either DFT or ReaxFF (0.4 eV for both). In the case of Cu−C bonds, the potential is fit to several selected molecular structures such as the CuC dimer, CCu2, C2Cu, and copper metallocene. In the case of CuH, the potential is fit to the stretching PES of CuH only. Though the training database includes information from a wide range of pure and binary structures, many force fields that occur at interfaces, such as C−Cu−Cu or C−C−Cu bond angles at Cu/graphene interfaces, are excluded from the training database. However, one of the primarily goals of developing this potential is to apply COMB to multifunctional and multicomponent systems. So far, the interactions that have been considered only involve one or two different types of elements. This is addressed as follows with a proper description of the interactions within the pure and binary systems, such that the interactions that are missed in the training database in most cases can be manifested appropriately. For example, the model can properly describe the force fields involved in Cu metallocene, a ternary system, without any three-dimensional parameters between Cu−C−H. A more accurate description of these excluded many-body interactions can be achieved by parametrizing the correction terms with the opening threedimensional parameters for specific applications. The only special case where a different approach is needed is the interactions in multicomponent organic materials such as, for example, CHONF systems. As described in section IIB, COMB3 employs spline functions to account for coordination effects within organic systems. The spline functions are limited to be three-dimensional, in which the first dimension is the coordination number of carbon atoms, the second dimension is the coordination number of the atoms with lower electronegativity than carbon, and the third is the neighboring atoms that are more electronegative than carbon. To effectively reduce the dimensions for such multicomponent organic systems, twodimensional parameters Ncross are employed to account for the il element dependence on the three dimensions in the spline functions. With such a design, the COMB3 potential should be capable of capturing the fundamental chemistry and physics associated with a multicomponent hydrocarbon system. It is worthwhile to compare the fitting and performance of ReaxFF and COMB3 potentials that were both designed to model similar heterogeneous multicomponent systems. We therefore compare the COMB3 potential presented here with ReaxFF potentials parametrized for C−H38 and Ni−C−H43,44 systems. With regard to potential fitting, the many-body short-range interactions in the ReaxFF potentials are expressed in a separated and charge-independent energy term in which the parameter dimension for N-body interactions are specified to be N. In contrast, the COMB3 potential integrated the many-

IV. RESULTS AND DISCUSSION A. Predicted Structures and Energetics. Table 7 compares REBO2,18 ReaxFF,38 and COMB3 values for the heat of formation with the experimental values for an extensive set of hydrocarbons that include radicals, saturated, and large unsaturated species. In general, the heats of formation predicted by the COMB3 potential deviate from experimental values by no more than 15 kcal/mol. The root-mean-square errors for COMB3, ReaxFF, 7987



Article

Thus, compared with quantum mechanically based approaches, the COMB3 potential provides a reasonably accurate and much faster approach for calculating the properties of heterogeneous systems. This opens up new possibilities for simulating larger systems than may be considered with electronic structure methods, as exemplified in the next section. B. Ethyl Radical Deposition on the Cu(111) Surface. The deposition of polyatomic ions or clusters at the hyperthermal level (collision energies ranging from 20 to 2000 kcal) is widely used for thin-film growth and surface modification. Here, atomistic simulations are employed to examine the chemistry at the surface following cluster−surface impacts. To feature the organic−inorganic interactions, ethyl radicals or •CH3CH2 are selected as the incident molecules and the Cu(111) surface is chosen to be the substrate used for the deposition. To validate the potentials, we have calculated the binding energies of ethyl radicals on Cu surfaces and compared them with the results of DFT calculations from ref 86; the results are given in Table 9.

and REBO2 are 5.72, 5.14, and 8.50 kcal/mol, respectively, which indicates that these three potentials are comparable at predicting the heat of formations of these molecules. Table 8 gives the fitting results of the heat of formation at 0 K, carbon−copper bond length and the charge on Cu of CuC Table 8. COMB3 Properties of CuC, CuC2, and Cu Metallocene heat of formation (kcal/mol) molecules

B3LYP

CuC CuC2 Cu metallocene

184 323 55

C−Cu length (Å)

COMB3 B3LYP 182 295 55

1.77 1.85 2.28

charge on Cu (e)

COMB3 B3LYP 1.89 1.93 2.25

0.37 0.71 1.18

COMB3 0.17 0.28 0.65

dimer, CuC2, and Cu metallocene. To directly compare the carbon−metal interaction in organometallics, the energy of the dissociation reaction of the Cu metallocene into two cyclopentadienyls and a Cu atom is calculated, as shown in Figure 3. The dissociation energy that is illustrated in Figure 3 is 122 kcal/mol for COMB3 and 122 kcal/mol for B3LYP.

Table 9. Binding Energies of •CH3CH2 on Cu Surfaces

a

surface

DFTa (kcal/mol)

COMB3 (kcal/mol)

Cu(100) Cu(110) Cu(111)

−36.2 −39.9 −30.9

−35.5 −39.3 −32.5

Reference 86.

The predicted binding energies from DFT calculations are obtained for a half monolayer (ML) of CH3CH2. Table 9 indicates that COMB3 predicts binding energies of −35.5, −39.3, and −32.5 kcal/mol for Cu(100), Cu(110), and Cu(111) surfaces, respectively, which are in good agreement with the predictions of DFT calculations. The good agreement indicates the ability of the former to model organic−inorganic surface chemistry. To more realistically analyze low energy deposition, a Cu(111) surface with dimensions of 3.1 × 3.1 × 5.0 nm (Figure 4) is considered within classical MD simulations to model ethyl radical deposition. Such deposition simulations are possible but extremely time-consuming for smaller systems with DFT-MD. As indicated in Figure 4, the bottom three Cu layers are fixed, the next six layers are thermostatted with a Nose− Hoover algorithm,87,88 and the top 15 layers are active, or have no applied constraints. Within these 15 active layers, the atoms

Figure 3. Dissociation of Cu metallocene in COMB3 (colored by charge).

As indicated in Figure 3, the carbon−metal interaction within organometallics such as copper metallocene, is a hybrid of ionic, covalent, and metallic bonding. From the molecular orbital point of view, the bonding between organic and inorganic fragments can be explained by the isolobal analogy proposed by Hoffmann,84 which is that these fragments must have similar numbers of electrons, orbital shapes, symmetry properties, and approximate energy levels of the frontier orbitals to form organic−inorganic bonds. This analogy is also implied in the 18-electron rule, which is a rule of thumb used to predict the stability of the metal carbonyls in organometallic chemistry. Both the isolobal analogy and the 18-electron rule85 have structural implications, especially for the coordination number of the structures. This structural implication is empirically approximated in COMB3 by the coordination number dependence in the field effects (eq 5) to adjust the charge and in the bond order term (eq 41). With these dependences on coordination number, COMB3 is capable of reproducing the charge trend, bond lengths, and energetics of interactions between Cu and carbon atoms from the B3LYP calculations, as shown in Table 8 and Figure .3. For example CuC and CuC2 have heats of formation of 182 and 295 kcal/mol, which indicate the Cu−C bonds inside these two low-coordinated molecules are unstable. In contrast, the Cu−C bonds inside the Cu metallocene, where the coordination number of Cu is 10, have binding energies of −122 kcal/mol or −5.3 eV/mol, which suggests that relatively strong Cu−C bonds are formed in Cu metallocene.

Figure 4. System configuration of ethyl radical deposition on the Cu(111) surface. 7988



Article

Figure 5. Snapshots for CH3CH2 dissociation into CH3 and CH2 on the Cu(111) surface. All the images are colored by charge according to the scale.

Figure 6. Snapshots for CH3CH2 impacting another CH3CH2 on the Cu(111) surface (colored by charge). (a, b) Longer carbon chain formation. (c−f) C interstitial (highlighted with an arrow) formation in the Cu tetrahedral interstitial site. (g, h) Cu atom knocked out from the surface to accommodate a C interstitial.

reaction CH3CH2 ⇒ CH3 + CH2 in the absence of Cu is predicted to be about 110 kcal/mol by DFT and about 100 kcal/mol by COMB. In the presence of a single Cu atom, this energy is predicted to be substantially decreased to about 90 kcal/mol by DFT and to about 85 kcal/mol by COMB. In Figure 6a, the center two molecules are the CH3CH2 radicals that are involved in chemical reactions, whereas the three molecules to the right are previously deposited molecules. When the CH2 carbon atom from one CH3CH2 hits the CH2 carbon atom from another CH3CH2 radical (Figure 6a), CH3CH2CH2CH3, n-butane, is formed as indicated in Figure 6b. As indicated in Table 7, the heat of formation for reaction CH3CH2 + CH2CH3 ⇒ CH3CH2CH2CH3 is about −80 kcal/ mol from experimental data and −76 kcal/mol as predicted by the COMB3 potential. Though the formation of n-butane is energetically favorable, only one n-butane is formed out of 12 incident ethyl molecules, which suggests the formation reaction above is dependent on the relative positions of the ethyl radicals on the surface. The energy released from the formation reaction together with the incident kinetic energy is sufficient to inject a C atom into the Cu surface, break the carbon−carbon bond and thereby form a C interstitial, as indicated in Figure 6c−f. To accommodate the strain introduced by a C interstitial, one Cu atom is slightly lifted out from the surface as shown in Figure 6g−h.

in the topmost six layers are allowed to undergo charge transfer whereas the atoms in the next nine layers are held at fixed neutral charge. Two beams of 12 CH3CH2 molecules are deposited in the direction normal to the surface at two incident kinetic energies of 115 and 230 kcal, i.e., 5 and 10 eV, per molecule. The relaxation time between deposition events in both beams is set to 0.5 ps and the temperature of the surface is maintained at 300 K throughout. In the simulations, 10 out of 12 incident ethyl radicals deposited with kinetic energies of 115 kcal, and 9 out of 12 radicals deposited with energies of 230 kcal, were predicted to remain on the surface. When a CH2 on an ethyl radical directly impacts a Cu atom, or a C atom already on the surface, chemical reactions are predicted to occur, as illustrated in Figures 5 and 6. In particular, when the carbon atom in CH2 impacts a Cu atom on the surface (Figure 5a), the instantaneous charges are +0.5 e on Cu, −0.3 e on the CH2, and −0.4 e on the CH3 of the CH3CH2, and this accounts for the charge-associated colors in this figure. The charge transfer and the resulting increase in the electrostatic energies are adequate to break the carbon−carbon bond in CH3CH2, as indicated in Figure 5b−d. This charge transfer and the decrease in the dissociation energy are in qualitative agreement with the predictions of DFT (B3LYP/6311G*) of a CH3CH2 interacting with a single Cu atom, where the charges change by +0.29 e on the Cu, −0.05 e on the C on the CH3, and −0.30 e for the C of the CH2. The dissociation 7989



Article

and Simulation of Nuclear Reactors under U.S. Department of Energy Contract No. DE-AC05-00OR22725.

This case study demonstrates the capabilities of COMB potentials for simulating dynamic reactions associated with organic−inorganic surface chemistry.

■

IV. CONCLUSIONS The third generation of the charge optimized many body (COMB3) potential, with inspiration from the MoS2 REBO potential and the second generation REBO2 potential for hydrocarbons, has improved on the previous versions of the COMB formalisms by adding a coordination function, a radical contribution, and a dihedral term to the bond-order term and by changing the expressions of the self-energy term. These enhancements enable COMB3 to appropriately and autonomously determine the charge state of an atom and largely capture the basic chemistry and physics of the different bond types including unconjugated, conjugated, and radical-containing organic compounds, the complex bonding in organometallics, as well as the ionic, metallic, and van der Waals bonds that can coexist in the multicomponent and multifunctional systems. A step by step parametrization scheme is provided, which is helpful to understand the origins of and the physics behind each potential energy term and provides a general guideline for the parameterizing the COMB potentials. The fitted and predicted results from COMB3 for solid carbon materials, hydrocarbons, and Cu−hydrocarbon systems are compared with available REBO2, ReaxFF, and quantum chemical calculations and experimental results. This comparison shows that the COMB3 potential is indeed able to reproduce the energetics, geometries, and elastic properties associated with the nonreactive and reactive behavior of these compounds. This new potential enables the modeling of computational chemistry in heterogeneous, many-atom systems at system sizes that are larger than those that are currently tractable with firstprinciples approaches. A case study of low energy deposition of ethyl radicals on the Cu(111) surface is presented. Several surface reactions and the associated charge transfer instances are properly predicted by COMB3. This case study example further demonstrates the strength of COMB3 in describing bond formation and breaking between organic and inorganic species in a constantly changing local environment.

■

(1) Jones, J. E. Proc. R. Soc. London, Ser. A 1924, 106, 441. (2) Jones, J. E. Proc. R. Soc. London, Ser. A 1924, 106, 463. (3) Daw, M. S.; Baskes, M. I. Phys. Rev. Lett. 1983, 50, 1285. (4) Daw, M. S.; Baskes, M. I. Phys. Rev. B 1984, 29, 6443. (5) Buckingham, R. A. Proc. R. Soc. London, Ser. A 1938, 168, 264. (6) Lewis, G. V.; Catlow, C. R. A. J. Phys. C: Solid State Phys. 1985, 18, 1149. (7) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; K.M., M. J.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (8) Brooks, B. R.; Brooks, C. L., III; Mackerell, A. D.; Nilsson, L.; Petrella, R. J.; Roux, B.; Won, Y.; Archontis, G.; Bartels, C.; Boresch, A.; Caflisch, S.; Caves, L.; Cui, Q.; Dinner, A. R.; Feig, M.; Fischer, S.; Gao, J.; Hodoscek, M.; Im, W.; Kuczera, K.; Lazaridis, T.; Ma, J.; Ovchinnikov, V.; Paci, E.; Pastor, R. W.; Post, C. B.; Pu, J. z.; Schaefer, M.; Tidor, B.; Venable, R. M.; Woodcock, H. L.; Wu, X.; Yang, W.; York, D. M.; Karplus, M. J. Comput. Chem. 2009, 30, 1545. (9) Allinger, N. L. J. Am. Chem. Soc. 1977, 99, 8127. (10) Gunsteren, W. F.; Huenberger, P. H.; Mark, A. E.; Smith, P. E.; Tironi, I. G. Comput. Phys. Commun. 1995, 91, 305. (11) Stillinger, F. H.; Weber, T. A. Phys. Rev. B 1985, 31, 5262. (12) Abell, G. C. Phys. Rev. B 1985, 31, 6184. (13) Tersoff, J. Phys. Rev. Lett. 1986, 56, 632. (14) Tersoff, J. Phys. Rev. Lett. 1988, 61, 2879. (15) Tersoff, J. Phys. Rev. B 1988, 37, 6991. (16) Tersoff, J. Phys. Rev. B 1988, 38, 9902. (17) Brenner, D. W. Phys. Rev. B 1990, 42, 9458. (18) Brenner, D. W.; Shenderova, O. A.; Harrison, J. A.; Stuart, S. J.; Ni, B.; Sinnott, S. B. J. Phys.: Condens. Matter 2002, 14, 783. (19) Schall, J. D.; Gao, G.; Harrison, J. A. B. Phys. Rev. B 2008, 77, 115209. (20) Ni, B.; Lee, K.-h.; Sinnott, S. B. J. Phys.: Condens. Matter 2004, 16, 7261. (21) Fonseca, A. F.; Lee, G.; Borders, T. L.; Zhang, H.; Kemper, T. W.; Shan, T.-R.; Sinnott, S. B.; Cho, K. Phys. Rev. B 2011, 84, 075460. (22) Jang, I.; Sinnott, S. B. J. Phys. Chem. B 2004, 108, 18993. (23) Brenner, D. W. Phys. Rev. Lett. 1989, 63, 1022. (24) Liang, T.; Phillpot, S. R.; Sinnott, S. B. Phys. Rev. B 2009, 79, 245110. (25) Stuart, S. J.; Tutein, A. B.; Harrison, J. A. J. Chem. Phys. 2000, 112, 6472. (26) Sinnott, S. B. J. Nanosci. Nanotech. 2002, 2, 113. (27) Hanley, L.; Sinnott, S. B. Surf. Sci. 2002, 500, 500. (28) Hsu, W.-D.; Tepavcevic, S.; Hanley, L.; Sinnott, S. B. J. Phys. Chem. C 2007, 111, 4199. (29) Pearson, J. D.; Gao, G.; Mohammed, A. Z.; Harrison, J. A. Comput. Mater. Sci. 2009, 47, 1. (30) Heo, S.-J.; Sinnott, S. B. J. Appl. Phys. 2007, 102, 064307. (31) Sinnott, S. B.; Heo, S.-J.; Brenner, D. W.; Harrison, J. A. Computer simulations of nanometer-scale indentation and friction. In Springer handbook on nanotechnology; Bhushan, B., Ed.; SpringerVerlag: Heidelberg, Germany, 2007; pp 1051. (32) Nguyen, T. C.; Ward, D. W.; Townes, J. A.; White, A. K.; Krantzman, K. D.; Garrison, B. J. J. Phys. Chem. B 2000, 104, 8221. (33) Baskes, M. I. Phys. Rev. B 1992, 46, 2727. (34) Watanabe, T.; Fujiwara, H.; Noguchi, H.; Hoshino, T.; Ohdomari, I. Jpn. J. Appl. Phys. 1999, 38, L366. (35) Billeter, S. R.; Curioni, A.; Fischer, D.; Andreoni, W. Phys. Rev. B 2006, 73, 155329. (36) Streitz, F. H.; Mintmire, J. W. Phys. Rev. B 1994, 50, 11996. (37) Zhou, X. W.; Wadley, H. N. G.; Filhol, J.-S.; Neurock, M. N. Phys. Rev. B 2004, 69, 035402. (38) van Duin, A. C. T.; Dasgupta, S.; Lorant, F.; Goddard, W. A. J. Phys. Chem. A 2001, 105, 9396.

ASSOCIATED CONTENT

S Supporting Information *

Tables of one- and two-dimensional parameters in COMB3, C−C bond and angular parameters, knot values, and threedimensional parameters in eq 50. This material is available free of charge via the Internet at http://pubs.acs.org.

■

REFERENCES

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]fl.edu. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS T.L., B.D., and S.B.S. gratefully acknowledge the support of the National Science Foundation (CHE-0809376). The work of S.R.P. was supported by the Consortium for Advanced Simulation of Light Water Reactors (www.casl.gov), an Energy Innovation Hub (http://www.energy.gov/hubs) for Modeling 7990



Article

(39) Zhang, Q.; Cagin, T.; Van Duin, A.; Goddard, W. A., III. Phys. Rev. B 2004, 69, 045423. (40) van Duin, A. C. T.; Strachan, A.; Stewman, S.; Zhang, Q.; Xu, X.; Goddard, W. A., III. J. Phys. Chem. A 2003, 107, 3803. (41) Chenoweth, K.; van Duin, A. C. T.; Goddard, W. A., III. J. Phys. Chem. A 2008, 112, 1040. (42) Chakraborty, D.; Muller, R. P.; Dasgupta, S.; Goddard, W. A., III. J. Phys. Chem. A 2000, 104, 2261. (43) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A., III. J. Phys. Chem. C 2010, 114, 4939. (44) Mueller, J. E.; van Duin, A. C. T.; Goddard, W. A., III. J. Phys. Chem. C 2010, 114, 5675. (45) Liu, B.; Lusk, M. T.; Ely, J. F. Surf. Sci. 2012, 606, 615. (46) Meng, L.; Sun, Q.; Wang, J.; Ding, F. J. Phys. Chem. C 2012, 116, 6097. (47) Watanabe, T. J. Comput. Electron. 2011, 10, 2. (48) Yasukawa, A. JSME Int. J. Ser. A 1996, 39, 313. (49) Yu, J. G.; Sinnott, S. B.; Phillpot, S. R. Phys. Rev. B 2007, 75, 085311. (50) Yu, J.; Sinnott, S. B.; Phillpot, S. R. Philos. Mag. Lett. 2009, 89, 136. (51) Shan, T.-R.; Devine, B. D.; Hawkins, J. W.; Asthagiri, A.; Phillpot, S. R.; Sinnott, S. B. Phys. Rev. B 2010, 82, 235302. (52) Shan, T.-R.; Devine, B. D.; Phillpot, S. R.; Sinnott, S. B. Phys. Rev. B 2011, 83, 115327. (53) Shan, T.-R.; Devine, B. D.; Kemper, T. W.; Sinnott, S. B.; Phillpot, S. R. Phys. Rev. B 2010, 81, 125328. (54) Devine, B.; Shan, T. R.; Cheng, Y. T.; McGaughey, A. J. H.; Lee, M.; Phillpot, S. R.; Sinnott, S. B. Phys. Rev. B 2011, 84, 17. (55) Mortier, W. J.; van Genechten, K.; Gasteiger, J. J. Am. Chem. Soc. 1985, 107, 829. (56) Toufar, H.; Nulens, K.; Janssens, G. O. A.; Mortier, W. J.; Schoonheydt, R. A.; DeProft, F.; Geerlings, P. J. Phys. Chem. 1996, 100, 15383. (57) Slater, J. C. Phys. Rev. 1930, 36, 57. (58) Stern, H. A.; Kaminski, G. A.; Banks, J. L.; Zhou, R. H.; Berne, B. J.; Friesner, R. A. J. Phys. Chem. B 1999, 103, 4730. (59) Wolf, D.; Keblinski, P.; Phillpot, S. R.; Eggebrecht, J. J. Chem. Phys. 1999, 110, 8254. (60) Thijsse, B. J. Nucl. Instrum. Methods Phys. Res. B 2005, 228, 198. (61) http://cms.mpi.univie.ac.at/vasp/. (62) Blochl, P. E. Phys. Rev. B 1994, 50, 17953. (63) Kresse, G.; Joubert, J. Phys. Rev. B 1999, 59, 1758. (64) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian 09, Revision A.02; Gaussian, Inc.: Wallingford, CT, 2009. (65) Becke, A. D. J. Chem. Phys. 1993, 98, 5648. (66) Lee, C.; Yang, W.; Parr, R. G. Phys. Rev. B 1988, 37, 785. (67) Bernardi, F.; Bottini, A.; McDougall, J. J. W.; Robb, M. A.; Schlegel, H. B. Faraday Symp. Chem. Soc. 1984, 19, 137. (68) Frisch, M. J.; Ragazos, I. N.; Robb, M. A.; Schlegel, H. B. Chem. Phys. Lett. 1992, 189, 524. (69) Wiberg, K. B.; Hadad, C. M.; LePage, T. J.; Breneman, C. M.; Frisch, M. J. J. Phys. Chem. 1992, 96, 671. (70) Purvis, G. D.; Bartlett, R. J. J. Chem. Phys. 1982, 76, 1910. (71) Dunning, T. H. J. Chem. Phys. 1989, 90, 1007.

(72) Yin, M. T.; Cohen, M. L. Phys. Rev. Lett. 1983, 50, 2006. (73) Yin, M. T.; Cohen, M. l. Phys. Rev. B 1984, 29, 6996. (74) Donohue, J. The Structures of Elements; Wiley: New York, 1974. (75) Zacharia, R.; Ulbricht, H.; Hertel, T. Phys. Rev. B 2004, 69, 155406. (76) Mounet, N.; Marzari, N. Phys. Rev. B 2005, 71, 205214. (77) Blakslee, O. L.; Proctor, D. G.; Seldin, E. J.; Spence, G. B.; Weng, T. J. Appl. Phys. 1970, 41, 3373. (78) Seldin, E. J.; Nezbeda, C. W. J. Appl. Phys. 1970, 41, 3389. (79) Grimsditch, M. H.; Ramdas, A. K. Phys. Rev. B 1975, 11, 3139. (80) Klintenberg, M.; Lebegue, S.; Ortiz, C.; Sanyal, B.; Fransson, J.; Eriksson, O. J. Phys.: Condens. Matter 2009, 21, 335502. (81) http://webelements.com. (82) Plimpton, S. J. J. Comput. Phys. 1995, 117, 1. (83) http://cccbdb.nist.gov. (84) Hoffmann, R. Angew. Chem., Int. Ed. 1982, 21, 711. (85) Jensen, W. B. J. Chem. Educ. 2005, 82, 28. (86) Li, X.; Gellman, A. J.; Sholl, D. S. J. Chem. Phys. 2007, 127, 144710. (87) Hoover, W. G. Phys. Rev. A 1985, 31, 1695. (88) Nose, S. J. Chem. Phys. 1984, 81, 511.

7991


Variable Charge Reactive Potential for Hydrocarbons to Simulate

Recommend Documents