eReaxFF: A Pseudoclassical Treatment of Explicit Electrons within Reactive Force Field Simulations

: We present a computational tool, eReaxFF, for simulating explicit electrons within the framework of the standard ReaxFF reactive force ﬁ eld method. We treat electrons explicitly in a pseudoclassical manner that enables simulation several orders of magnitude faster than quantum chemistry (QC) methods, while retaining the ReaxFF transferability. We delineate here the fundamental concepts of the eReaxFF method and the integration of the Atom-condensed Kohn − Sham DFT approximated to second order (ACKS2) charge calculation scheme into the eReaxFF. We trained our force ﬁ eld to capture electron a ﬃ nities (EA) of various species. As a proof-of-principle, we performed a set of molecular dynamics (MD) simulations with an explicit electron model for representative hydrocarbon radicals. We establish a good qualitative agreement of EAs of various species with experimental data, and MD simulations with eReaxFF agree well with the corresponding Ehrenfest dynamics simulations. The standard ReaxFF parameters available in the literature are transferrable to the eReaxFF method. The computationally economic eReaxFF method will be a useful tool for studying large-scale chemical and physical systems with explicit electrons as an alternative to computationally demanding QC methods.


I. INTRODUCTION
Electrons, which are intrinsically of a quantum nature, govern all forms of molecular structures and chemical reactions.In order to investigate the mechanism of chemical reactions and physical properties of substances at an atomistic level, computational chemistry has been used extensively as a complement to experimental techniques over at least the last four decades.Historically, the development of computational methodologies for describing chemical reactions at an atomistic scale has primarily concentrated on quantum chemistry (QC) methods.These methods, which possess explicit electronic degrees of freedom, provide the most accurate and detailed description of chemical reactions, but their inherent complexity in formalism and high computational cost limit applications to relatively small length-and time-scales.Empirical formalisms have been introduced in QC methods to improve computational efficiency, at the expense of accuracy.Approximate QC based methodsfor example, density functional theory (DFT) or Tight-Binding (TB) schemeshave been developed as computationally advantageous alternatives to truly ab initio methods.Even with such approximations, nonadiabatic dynamics using time-dependent DFT are quite expensive and are limited to a relatively small number of atoms (∼1000) and short time scales (∼1 ps). 1,2orce field (FF) based methods have also been developed as an alternative to QC methods.FF methods use simplified functional forms and provide much better computational performance, enabling simulations to access much larger scales, of order nanometers in length and nanoseconds in time.However, most of the available force field methods provide a substitute for the electronic structure and are typically nonreactive in nature, implying the absence of explicit electrons and their impact on system geometry and energy, as well as the inability to simulate bond breaking and formation.To bridge the gap between minimal computational cost and the ability to simulate reactions, a genre of force field methods known as "reactive force fields" was developed.The reactive force field methods 3−6 possess the capability of simulating reactions onthe-fly; however, the implicit treatment of electrons is inadequate for the description of many physical systems, such as redox reactions of rechargeable battery interfaces, polarization behavior in ferro-/piezo-electric materials, and solar cells.An accurate description of these phenomena requires an explicit treatment of electron degrees of freedom, which must be within the classical framework of the classical potentials.
Recently, a number of force field methods have been developed that include aspects of explicit valence electrons, like the electron force field (eFF) 7,8 and LEWIS 9−11 force field.The eFF has been parametrized only for a limited number of elements in the periodic table, and the primary application areas are materials that are subjected to extreme pressure and temperature apart from the complexity of the reactive systems. 7,8The ability of the LEWIS method to describe atomic electron affinities (EA) and ionization potentials (IP) of the first three row elements is encouraging, still such description for molecules is not available. 12,13The lowtemperature dynamics of liquid water 9,11 and prediction of ground state configurations of some chemical species 14 indicate the promising capabilities of the LEWIS method, but, the challenge of describing complex reactions or high-temperature dynamics with an acceptable level of accuracy is yet to be demonstrated.There also exists alternate methods to couple covalent and electrostatic interactions through split-charge equilibration (SQE) 15,16 for describing charge transfer in reactive dynamics, 17 without an explicit electron treatment.To the best of our knowledge, no classical force field method with access to explicit electrons has been employed to simulate complex reactive events.
−24 ReaxFF has relatively sophisticated schemes to treat electrons in its energy functionalnot explicitly rather implicitly it accounts for the effect of electron distribution in chemical bonding.Such as an over-or undercoordination energy term performs by evaluating the difference between the number of electrons around an atom involved in a bonded interaction and the number of valence electrons of the atom, a polarizable charge calculation scheme obtains charges on every atom based on atomic electronegativity and hardness, and a lone-pair term includes the energy effects of breaking up a lone electron pair on an atom. 18ReaxFF is computationally significantly less intensive than DFT or TB methods and is capable of nanoseconds of MD simulations with millions of atoms. 25−30 The inability of ReaxFF to give accurate EAs or IPs for various molecules restricts its applicability to the description of redox reactions. 31A detailed investigation of electron flow associated reactions or electron dynamics requires further extensions to the ReaxFF concept.
To extend the ReaxFF method, we introduce an explicit electron-like or hole-like particle description into this method.The treatment of an explicit electron within the classical framework can only capture the "particle" nature of the "waveparticle dual" nature of the electron.The wave properties are treated implicitly in the functionals of the ReaxFF method.This electron version of the ReaxFF method will be addressed as "eReaxFF".In this study, we perform the force field training, and as a proof of concept we consider two model hydrocarbon radicals and study electron transfer (ET) dynamics to compare our eReaxFF results with those obtained from Ehrenfest dynamics simulations.
This paper is organized as follows.In section II, we briefly introduce the fundamental concepts of the ReaxFF method, followed by the general theory of the eReaxFF method and implementation of the Atom-condensed Kohn−Sham DFT approximated to second order (ACKS2) 32 charge calculation scheme.Section III describes the relevant force field training results.In section IV, we present an application of the eReaxFF method to the model hydrocarbon radicals that includes highlights from the eReaxFF MD simulations and a comparison with the Ehrenfest dynamics results.

II. COMPUTATIONAL METHODS
II.A. ReaxFF Description.ReaxFF is a general bond order 3,33 19 This combination of covalent and Coulomb interactions has made ReaxFF applicable to a wide range of materials, including covalent, 27,30,34 metallic, 23 and liquid−solid interfaces. 28The electrostatic interactions are decoupled from the covalent contributionsi.e.charge is not coupled with the valency of an atom type.ReaxFF uses geometry dependent charge calculation schemes.Most of the currently published ReaxFF force field parametrizations are based on the Electronegativity Equalization Method (EEM) 35 method.We recently developed an extension of the EEM-method − the Atom-condensed Kohn− Sham DFT approximated to second order (ACKS2) 32 − for calculations of atomic charges. 22II.B.General Theory of eReaxFF.We developed our eReaxFF method within the basic framework of the ReaxFF.All 2-, 3-, and 4-body and nonbonded interaction terms of ReaxFF are retained.Next, we added new energy functionals to compute pairwise electrostatic interactions for explicit electrons and modified existing over-and undercoordination energy terms.
We incorporated a limited pseudoclassical explicit electron/ hole degrees-of-freedom scheme that is complementary to the implicit treatment of electrons in the bonded interactions of ReaxFF.The electron or hole is represented as an additional particle that carries a −1 (electron) or +1 (hole) charge, respectively.The nuclei are treated as point charges and electrons as Gaussian wave functions (ψ∞ exp(−α(r−r′) 2 ).The energy is a function of the position of the electron and atomic-centers and depends on the spread of the Gaussian function, α.The pairwise electrostatic interaction between the electron and core-charge is described as 7 where Z i is the nuclear charge, R ij is the distance between the electron and nucleus, and α (Gaussian exponent) and β are constants that depends on the atom type.The core-charge is the charge corresponding to the atomic number of an atom.Electron−electron interactions are treated through Coulomb point charges and short-range Gaussian repulsion functions.
In ReaxFF, valency and number of lone-pair electrons are treated as constants for an atom type.It is necessary to allow changes in the atom valency and number of lone pair of electrons when explicit electron/hole degrees of freedom are introduced into the system.The diffuse nature of an added explicit electron/hole and its relation to a particular atom is described as a function of the distance between electron/hole and the nuclear position of the host atom.We chose an exponential function to determine the number of electrons in the host atom; at the same time, this function ensures that an electron can virtually split itself among its neighboring atoms, that is, it resembles a partial delocalization in a molecule.This function has the following form where R ij is the distance between the atom-center and the electron/hole, and p val is a general parameter in the force field.
The effect of an explicit electron in the vicinity of an atom can be explained as follows: an electron uptake by an oxygen atom, for instance, decreases its valency and increases its number of valence electrons by one.This changed oxygen atom behaves like a pseudofluorine atom, that is, it acquires an electron configuration similar to that of fluorine.The acceptance of a hole acts in the opposite manner, the number of outer shell electrons is reduced, and the valency of the oxygen atom increases from two to three, thus the electron deficient oxygen atom resembles the electron configuration of nitrogen.Accordingly, we have established a set of rules in the eReaxFF code that dictates the effect of explicit electron/hole on the respective atom types.The atom valency and the number of lone-pair electrons become a dynamic variable, and they are evaluated at every step of the dynamics.
In ReaxFF, the bonded and nonbonded interactions are calculated independently, that is, charge is not coupled with the number of valence electrons of an atom type which is used in bond order correction and over-and undercoordination energy terms.For example, a Li + -dimer is, in principle, able to make a Li−Li bonddespite the positive chargeas the number of valence electron remains one.With the charge-valency coupling, the Li + −Li + bond is additionally counteracted by the overcoordination terms, as the positively charged Li-atom loses its valence electron.In our eReaxFF implementation, the explicit charges of an atom are coupled with the atom valency, and the over-and undercoordination energies are calculated based on the corrected valence electrons.
Modification in the number of valence electrons of an atom type modifies the degree of over-or undercoordination (Δ i ), thus corresponding energy penalties.We explain how the degree of over-or undercoordination is calculated below.A decrease in the valency of an atom at an existing bonding environment increases the overcoordination, resulting in a larger overcoordination energy penalty, and therefore reduces the BO associated with the atom and weakens relevant bonds.
For example, in a methane molecule, an electron uptake by the carbon will decrease its valency from four to three, and at this bonding situation the carbon atom is further overcoordinated by ∼1 which leads to a significant overcoordination energy penalty and total bond order of carbon reduction to ∼3.0 resulting in the elongation of C−H bonds.
In our eReaxFF implementation, the amount of modified degree of over-or undercoordination used in the corresponding energy penalty calculation is designed as a function of both atom and bond types, that is, it is contingent on the atom type as well as on the local bonding environment.This atom and bond-type dependent treatment of over-and undercoordination enables the eReaxFF method to capture EAs or IPs of various species using the same set of atomic-Gaussian parameters.The modified over-and undercoordination energy functionals are given by We used new functional forms for an explicit electron correction on the degree of over/undercoordination.It is essential to choose a functional form that is efficient in computation and smoothly differentiable to ensure energy conservation.The former condition is satisfied by using simple exponentials, and the latter is ensured since the functional form is analytically differentiable.The explicit electron correction depends on the BOs and the electron occupancy of the host atom.The bond type parameter (p ij xel1 ) and the corresponding BOs are used to adjust the number of electrons available to the host atom (eq 4b).The atom type dependency is introduced through eq 4c, where p i xel2 is an atom type parameter.Finally, the explicit electron correction on the Δ i is expressed in eq 4d as Δ i xel .The Δ i xel is used in eq 4e, 4f, and 4g to calculate the explicit electron modified over-or undercoordination energy penalty.In our current eReaxFF implementation, a variable valency is considered in the over-and undercoordination and lone-pair energy functionals.In future developments, variable valency will be extended to the three-and four-body terms, that is, the valence angle and torsion energy expressions.The flow scheme of the eReaxFF method is shown in Figure 1.The fictitious mass of the electron and hole are assigned as equivalent to that of a hydrogen atom (1 a.m.u).Two other computational methodseFF 7 and LEWIS 9 where electrons are treated semiclassically, also uses 1 a.m.u as a mass of the electron.The eReaxFF code is currently available only in a serial version, and a parallel implementation in the commercial Amsterdam Density Function (ADF) 36 package is planned.II.C.ACKS2 Charge Calculation Scheme.The EEM methodwidely used within the ReaxFF approachis known for its spurious long-range charge transfer, for instance, between two or more molecules, even when they are wellseparated. 37It allows charge to be redistributed over all atoms as if the total system is an electric conductor. 38,39The longrange charge transfer of the EEM scheme is problematic for the eReaxFF method, as it impedes the accurate charge description of the reduced/oxidized molecules.The metallic polarizability of the EEM charges results in an almost complete charge compensation of the explicit electrons or holes.The ACKS2 charge calculation scheme 32 largely eliminates the issue of unrealistic long-range charge smearing.Therefore, we implemented this method in our eReaxFF code.Just like the EEM energy, the ACKS2 energy is quadratic in the atomic charges.Additional variables and quadratic energy terms are introduced to control the range over which charge can delocalize The variables are q i (atomic charges) and u i (fluctuations in the atomic Kohn−Sham potential) and are recomputed at each time step for the current geometry.The first two terms are identical to the EEM energy in the conventional ReaxFF, whereas the last two terms are new and account for nonlocal contributions to the electronic kinetic energy. 32In ReaxFF, χ i and η ii are atomic parameters, which are fixed for every chemical element in the calibration process.The parameter q k 0 is a reference charge: it is zero for all atoms and it is set to −1 or +1 for the explicit electron or hole, respectively.The remaining matrix elements depend on the interatomic distance R ij as follows: The off-diagonal element η ij represents the Coulomb interaction with short-range damping controlled by the atomic parameters γ i and γ j . 40The off-diagonal element X ij can be interpreted as a bond softness (or atom pair softness): it determines to what extent atom i and j can exchange charge directly.In the bonding region, X ij increases with increasing R ij until it reaches its maximum, X soft , and then goes smoothly to zero at the cutoff (C i + C j )/2.The atomic softness cutoff parameters C i and the universal softness parameter X soft are fixed in the calibration procedure, with the cutoff parameters remaining short-ranged.
The energy penalty for a direct charge transfer between atoms i and j is quadratic in charge transfer and is proportional to 1/X ij . 32When X ij is zero, charge transfer must pass through intermediate atoms, with similar energy penalties for all intermediate atom pairs.When no intermediate pairs are present, for example, between two separated molecules, no charge transfer is allowed, and the total charge of each isolated fragment is equal to the sum of the reference charges of the constituting particles.When two distant atoms reside in the same macromolecule, the number of required intermediate pairs, and thus also the energy penalty for charge transfer, scales linearly with the interatomic distance, which effectively inhibits long-range charge transfer.
The constrained minimization and maximization in the ACKS2 energy can be carried out simultaneously by solving a set of linear equations, which is similar to the EEM scheme.In addition, constraints are added to fix the charge of the explicit electron or hole, to −1 or +1, respectively.When some (or even all) X ij elements go to zero, the equations remain well behaved, which is attractive compared to the BOP/SQE model, in which charge transfer is disabled by letting a bond-hardness parameter diverge to + ∞. 38 The polarization catastrophe can be avoided by imposing the following inequalities on the parameters: η ii > γ i /4πε 0 and X soft > 0. Under these conditions, the matrices η ij and X ij are guaranteed to be positive and negative semidefinite, respectively, for any possible geometry.As such, the minimum and maximum in the ACKS2 energy are well-defined.
Originally, the ACKS2 model was derived from Kohn−Sham DFT and only included fluctuating atomic charges (and their conjugate potential variables).Recently, it was shown that a similar formalism can be derived from any variational electronic structure theory and that the ACKS2 model is easily extended with atomic multipole moments. 32lthough EEM and ACKS2 are conceptually similar, EEMbased ReaxFF force fields still need to be refitted to work with the ACKS2 method.Fortunately, given the similarities between EEM and ACKS2 similarities, this refitting is relatively straightforward−as the EEM-based ReaxFF parameters are typically a very good initial guess for the ACKS2/ReaxFF parameters.Apart from this EEM/ACKS2 transition, the current eReaxFF implementation is fully transferable with the ReaxFF parameter sets available in the literature.

III. FORCE FIELD TRAINING
The accuracy of the description of reduction reactions depends primarily on the ability of the force field method to reproduce EAs of various chemical species.We trained our eReaxFF force field to reproduce the EAs of various saturated, unsaturated, and radical species.The objective of choosing all three types of species is to demonstrate that eReaxFF can describe EAs in a wide range of bonding environments.The optimization of the parameters was performed via a successive one-parameter search technique 41 to minimize the following expression for the error where x lit is the target value for the EA, x eReaxFF is the eReaxFF calculated value, and σ i is the weight assigned to data point, i.
We performed training to optimize only those parameters that are related to the explicit electron description, that is, the atomic-Gaussian, as well as atom and bond parameters as they appear in eq 4d and retained all other original ReaxFF parameters.In eReaxFF, the electron affinity is defined as where E X and E X − are the energy of a species in a neutral state, and in a state with an additional electron, respectively, and E el is the energy of an electron.In ReaxFF simulations, the EA is calculated by subtracting the energy of a negatively charged species from the neutral one with the constraint of setting the molecular charge equal to −1e.In eReaxFF simulations, the EA is calculated by adding an electron to a molecule.In our description, the electron is a negatively charged particle; therefore, it is mostly localized at a particular site in the molecule during geometry optimization.As described in the previous section, a limited degree of delocalization is still available among the neighboring atoms through the modification of their valence.This delocalization influences the local Electron affinity data of various species as calculated using the eReaxFF method and comparison with the experimental 42,43 and DFT data.
chemistry of the molecule.Quantum calculations indicate a higher-degree of delocalization of the added electron over the entire molecule, compared to the eReaxFF method.This particular limitation is a result of the treatment of the electron as a pseudoclassical particle.
The eReaxFF force field training results of the EAs of various species are summarized in Figure 2 and compared to both ReaxFF and literature data. 42,43eReaxFF qualitatively reproduces the literature data of the EAs, while the ReaxFF completely fails to capture EAs of most of the species considered in the training set.eReaxFF significantly underestimates the EAs of the three unsaturated species, ethylene, propene, and isobutene, compared to experimental data. 43To investigate this discrepancy, we performed DFT calculations with the M06-2X/aug-cc-PVTZ basis set functional implemented in the Jaguar 7.5 program, 44 which exhibits that DFT also underestimates the EAs of these three species.Interestingly, eReaxFF results are in reasonable agreement with the DFT predicted results.Overall, eReaxFF provides an improved description of the EAs for a wide variety of species, which indicates that it can be effectively to study electron dynamics in hydrocarbon-related species.We next employ this newly developed force field to investigate electron transfer (ET) in model hydrocarbon radicals.The physical effects that dictate the ET can be explained using Marcus theory of the electron transfer. 46According to this theory, the rate of ET is determined by the electronic coupling between initial and final states, ⟨i|H|f⟩, the difference of their Gibbs free energies, the reorganization energy, and the temperature.In molecular dynamics method, the entropic and thermal contributions are incorporated directly, while force-field approaches parametrize the potential energies.In addition to improving on the latter, in eReaxFF, we attempted to parametrize and fit the effect of electronic coupling via the introduction of explicit electron and associated electron−nuclear interaction parameters.

IV. RESULTS AND DISCUSSION
IV.A. eReaxFF MD Simulations.We chose two representative hydrocarbon radicals, C 12 H 19 , to study electron dynamics, both of which consist of a conjugated part (polyacetylene), an aliphatic chain, and a radical site.These molecules are shown in Figure 3.To model an excited state, we inject an electron at the conjugated part of the radicals.Structural relaxation simulations were carried out using NVT (constant volume, temperature)-MD simulation at T = 1 K temperature.In the MD simulations, electrons are treated as  classical particles, and Newton's equations of motion are solved using the velocity Verlet algorithm.
We investigated the effect of temperature on the time-scale required for the electron to transfer from the polyacetylene to the radical site.The latter corresponds to the ground-state configuration for the electron in these species.We conducted NVT MD simulations at three different temperatures, 400, 500, and 600 K, with a Berendsen thermostat 45 and a damping constant of 100 fs.An MD time step of 0.1 fs was used.The potential energy profiles during the simulation of the C 12 H 19 • radical at different temperatures are shown in Figure 3(c).The electron diffuses to the polyacetylene or aliphatic chain at these temperatures with similar energy jumps.The localized electron at the radical site has a significantly lower energy than at the other two sites.From our MD simulations, we find a relatively faster ET from the polyacetylene part to the beginning of the aliphatic chain.Subsequently, the electron experiences slightly higher stability at the intersection of these two chains which reduces the ET rate from the intersection to the radical site.The increase in temperature enhances the probability of the electron to overcome the potential well, hence to transfer to the radical site.The potential energy profile at 600 K shows that, due to the higher thermal energy at this high temperature after the initial transfer of the electron from the polyacetylene to the radical site, the electron continues to hop between the three parts, polyacetylene, aliphatic, and radical sites.In Figure 4(a) we present the time-averaged electron localization around the contact point of the polyacetylene and aliphatic chains at different temperatures.With increasing temperature, electron localization around the contact point decreases, which illustrates the shorter time-scale requirement for ET at an elevated temperature.In Figure 5 (a) we present the time-scale required for the full ET process from the aliphatic to the radical site at various temperatures.Temperature plays a significant role in reducing the time-scale required for ET process.Snapshots from the ET process on the C 12 H 19 • radical are shown in Figure 5(b).
We also studied the effect of aliphatic chain length on the ET dynamics by considering two hydrocarbons with different length aliphatic chains, C 14 H 23 ̇• and C 12 H 19 • .Figure 5(a) illustrates these effect: the increase in the aliphatic chain length slows down the ET.In a similar vein, Figure 4 shows the timeaveraged electron localization around the contact point between the conjugated and aliphatic chains is higher in the radical with longer aliphatic chain.This higher electron localization can be attributed to the slower ET rate.Moreover, the longer aliphatic chain increases the diffusion length for the electron to arrive at the radical site, which also contributes to the larger time-scale required for the ET to the radical site.Likewise, in the C 12 H 19 • ̇case, increasing the temperature decreases the electron localization around the contact point, which results in an accelerated ET at higher temperatures.
IV.B. Ehrenfest Dynamics Simulations.In order to provide a qualitative validation of our eReaxFF MD simulation results, we carried out Ehrenfest dynamics simulations on the C 12 H 19 • radical.−51 In this approach, the nuclei are treated as classical particles, 52 while the electronic system is treated quantum-mechanically with real-time time-dependent density functional theory (rt-TDDFT). 45or the ED-rt-TDDFT simulations, we prepared the radical C 12 H 19 • by running ground-state MD simulations for 1.0 ps at 600 K. Then we introduced one extra electron into the system and found the electronic ground state without relaxing the geometry.Next, we simulated the electron injection (schematically presented in Figure 6(a)) into the polyacetylene part of the molecule by exciting the electron from the ground state of C 12 H 19 − into the first excited state.In the ground state about 0.6 of the extra electron charge is located on the aliphatic part of the radical molecule and 0.4 on the conjugated part.By exciting the molecule into the first excited state about −0.5 |e| is transferred to the conjugated part.We show in Figure 7, the obtained excited state as charge-density difference of the excited and ground states.
We then run ED-rt-TDDFT for ∼4.0 ps with 25 as (attosecond) time steps at 600 K and monitored the Hirshfeld charge on each atom.During the simulation, numerical errors were monitored, and the accumulated error was found to be small (total energy drift of 0.05 eV).We present snapshots of charge dynamics in Figure 6.Initially, the excited electron is localized in the polyacetylene part; after about 0.3 ps −0.1 |e| of charge is transferred to the adjacent CH 2 groups in the aliphatic part, with a small amount of charge appearing on the radical site (terminal CH 2 group); subsequently, the transfer is slow, steady, and almost linear in time, and at 3.0 ps about −0.2 |e| is transferred to the aliphatic side.In Figure 8 we show the timedependent excess Hirshfeld charge, apportioned into the polyacetylene and aliphatic (which includes radical site) parts.
The most interesting feature of the 600 K trajectory is the ultrafast electron transfer which is observed at about 3.5 ps time in trajectory.At this time about −0.2 |e| is transferred in a time interval of <7 fs.The analysis of potential energy surfaces around this point of the trajectory reveals that the system passes through the conical intersection of the first-excited and ground-state potential energy surfaces (see Figure 9), which facilitates the ultrafast electron transfer.We also observe that before passing through the conical intersection the electronic wavepacket already has substantial ground state character as can be seen in Figure 7(b).After the crossing the ground state retains the same character as before, while the electronic densities of the time-propagated and ground-state systems are virtually indistinguishable, see Figure 7(c), although the wavepacket retains some of the excited-state character which is manifested in coherent oscillations at the end of trajectory, see Figure 8(a).
To study the effect of temperature on the ET dynamics, we performed ED simulations at 300 K. Figure 8  Hirshfeld charges on the polyacetylene and aliphatic parts of the radical for this temperature.The slope of the charge transfer profile is higher in the case of 600 K compared to the 300 K, which indicates a faster rate of ET at an elevated temperature.
In order to compare our ED simulations to time scales estimated in eReaxFF simulations we fitted an exponential to the smooth part (leaving out ultrafast ET) of the 600 K ED time-dependent Hirshfeld charge shown in Figure 8    trajectories from ED simulations would have to be sampled to give a proper statistical average for a more meaningful comparison.We conclude that the qualitative trends of the ET time-scale as a function of temperature are comparable, and the estimate from the 600 K ED trajectory suggests even reasonable qualitative agreement.In both simulation techniques, an accelerated ET rate is observed with increasing temperature.We speculate that at higher temperature, the potential energy surface is sampled at a higher rate thus increasing the chance of encountering the region of high coupling as we observed in the ED trajectory.

V. CONCLUSIONS
In this work, we developed a pseudoclassical description of an explicit electron within the framework of the ReaxFF reactive force field.In the eReaxFF development, simple, yet robust functional forms are chosen to describe electron interactions so that this method can produce significant computational gains.The capabilities of the method are demonstrated by calculating electron affinities of various types of saturated, unsaturated, and radical species which qualitatively reproduce experimental values for all the species considered in this study.We studied electron transfer dynamics as a proof of concept for this newly developed method.The eReaxFF results are qualitatively in good agreement with results obtained from the Ehrenfest dynamics (rt-TDDFT) simulations.This agreement is very encouraging since computationally eReaxFF is several orders of magnitude faster than the Ehrenfest dynamics or other ab initio MD simulations.Future development of eReaxFF will be focused on the complex interfacial reactions of Li-ion batteries.Investigations are currently underway to describe reduction and oxidation reactions of the anode-and cathode-electrolyte interfaces, respectively, of Li-ion batteries with ethylene carbonate electrolytes.In these studies, we will further illustrate the dynamics of the both electron and hole.The current version of the eReaxFF method cannot describe the interactions between electron−hole pairs and as such is unable to describe the mechanism of the photon-driven processes such as photochemistry.We hope to incorporate this capability in future versions of the method.However, we anticipate that eReaxFF would be an effective method to develop a physicsbased description for piezo/ferro electric materials and their response to the application of an external field.eReaxFF should provide a powerful tool for understanding the dynamics of explicit electrons in physical and chemical systems where QC methods are still daunting computationally.

Notes
The authors declare no competing financial interest.

Figure 1 .
Figure 1.Flow diagram of the eReaxFF method: the covalent and nonbonded interactions are coupled through the explicit electron/hole.

Figure 2 .
Figure2.Electron affinity data of various species as calculated using the eReaxFF method and comparison with the experimental42,43 and DFT data.

Figure 3 .
Figure 3. Radicals used for the eReaxFF MD simulations (a) C 12 H 19 • and (b) C 14 H 23 • .(c) Potential energy profiles from the eReaxFF MD simulations on C 12 H 19• at three different temperatures, 400, 500, and 600 K. Blue and green shaded regions indicate electron localization on the polyacetylene/aliphatic and radical sites, respectively.Color scheme: black: carbon, and white: hydrogen.

Figure 4 .
Figure 4. Time-averaged electron localization around the contact point of the conjugated and aliphatic chains of the (a) C 12 H 19 ̇• and (b) C 14 H 23 • ṙadicals at three different temperatures, 400, 500, and 600 K.The blue sphere represents the electron.
(b) represents the

Figure 5 .
Figure 5. (a) The time-scale required for an electron to transfer from the polyacetylene to the radical site at different temperatures for the C 12 H 19 • and C 14 H 23 • radicals.(b) Snapshots of typical localization of an electron in the C 12 H 19 • radical during MD simulations.
(a), which gives a characteristic ET time of 11 ps and agrees very well with the eReaxFF predicted ET time of ∼11 ps at 600 K, presented in Figure 5(a).If we include the ultrafast ET part in the exponential fit, the time scale is reduced to ∼4 ps.Clearly many

Figure 6 .
Figure 6.(a) Electron injection process into C 12 H 19 • within ED-rt-TDDFT.The unpaired electron is depicted as an orange circle and is localized at the rightmost atom in the aliphatic part of the C 12 H 19 • radical.An electron is injected into the conjugated side of the radical.(b) Hirshfeld charge during the ED-rt-TDDFT at 600 K. Blue spheres represent the Hirshfeld atomic charge referenced to the neutral radical; the size of these blue spheres is proportional to the charge: (i) t = 0 ps, start of the Ehrenfest dynamics; the charge is mostly localized on the polyacetylene part of the radical.(ii) t = 0.3 ps, initial −0.1 |e| is transferred to the aliphatic chain.(iii) t = 3.0 ps, about −0.2 |e| charge is transferred the aliphatic chain.(iv) t = 4.0 ps, about −0.5 |e| charge is transferred to the aliphatic chain, and a significant portion is localized at the radical site.

Figure 7 .
Figure 7. Charge-density difference between the time-propagated density and the ground state at (a) the start of the trajectory; (b) just before the fast electron transfer; and (c) just after the fast electron transfer.Blue color depicts the excess of the negative charge ("electron"); green color depicts its depletion ("hole").

Figure 8 .
Figure 8. Hirshfeld charge on the polyacetylene and aliphatic chains of the radicals at (a) T = 600 K and (b) T = 300 K.The higher slope of the charge profile at 600 K indicates faster electron transfer rate compared to the 300 K simulation.

■
ACKNOWLEDGMENTS M.M.I., G.K., A.C.T.v.D., and E.K. acknowledge funding from a grant from the U.S. Army Research Laboratory through the Collaborative Research Alliance (CRA) for Multi-Scale Multidisciplinary Modeling of Electronic Materials (MSME).T.V. thanks the Research Board of Ghent University (BOF) for financial support.G.K. also acknowledges computational resources provided by XSEDE, Grant No. TG-DMR120073, which is supported by National Science Foundation Grant No. ACI-1053575.

Figure 9 .
Figure 9. Adiabatic electronic energies of the ground (black line) and excited (gray line) states along the 600 K trajectory around the time point where ultrafast electron transfer occurs.The inset shows an expanded version of the region of intersection of the potential energy surfaces, with the intersection occurring between t = 3.595 and 3.596 ps.
(BO) based classical force field method that allows on-the-fly bond breaking and formation during simulations.The general form of the ReaxFF energy terms are