eFF, a method to simulate large scale excited electron dynamics

Hydrogen molecule

A hydrogen molecule in eFF exists as two protons held together by an electron pair centered halfway between the protons. The bond length is 1.47 bohr, near the exact value of 1.40 bohr [19].

Squeezing the protons together, we find that the potential energy curve has a proper inner wall. Pulling apart the protons, we find that that at a certain point, each electron hops onto a separate proton, and thereafter the electrons sit on top of the nuclei instead of in the middle of the bond. In a Hartree-Fock wavefunction, this transition between open and closed-shell forms occurs at 2.3 bohr; in eFF, we observe the transition at 2.1 bohr.

The potential energy falls off a bit too sharply as the nuclei are separated. Like density functional theory, Hartree-Fock, and other single-determinant methods, eFF misses the resonance stabilization provided by exchanging electrons in the breaking bond.

Figure 8: Potential energy surface for $ \mathrm {H_{2}}$ dissociation.

The overall dissociation energy of the bond is 67 kcal/mol. In comparison, the bond energy computed using unrestricted Hartree-Fock (6-311g** basis) is 83 kcal/mol, and the exact bond energy is 109 kcal/mol - both the Hartree-Fock and the eFF calculations do not account for electron correlation, the instantaneous Coulomb avoidance of opposite spin electrons that reduces their repulsion.

Hydrogen reactions

eFF can distinguish between allowed and forbidden reactions involving hydrogen atoms and molecules. For the prototype allowed hydrogen exchange reaction $ \mathrm {H_{2}} + \mathrm {H} \rightarrow \mathrm {H} + \mathrm {H_{2}}$, eFF correctly finds that the transition state is colinear with a symmetric saddle point (at $ r = $ 1.04 $ \mathrm{\AA}$ eFF vs 0.95 $ \mathrm{\AA}$ quantum Monte Carlo/exact [20]).

The eFF energy barrier is 42 kcal/mol, which is too high relative to the exact energy barrier of 9.7 kcal/mol. This barrier is reduced to 9.2 kcal/mol in eFF3, a later version of the force field, whose starting point is parameterized against the uniform electron gas. However, we find that the original eFF model is still the most suitable for describing correctly the thermodynamics of warm dense hydrogen under a variety of temperatures and pressures.

Figure 9: Potential energy surface for hydrogen exchange reactions: (a) the allowed reaction $ \mathrm {H_{2}} + \mathrm {H} \rightarrow \mathrm {H} + \mathrm {H_{2}}$ and (b) the forbidden reaction $ \mathrm {H_{2}} + \mathrm {H_{2}'} \rightarrow \mathrm {H H'} + \mathrm {H H'}$.

For the four-center forbidden reaction $ \mathrm {H_{2}} + \mathrm {H_{2}'} \rightarrow \mathrm {H H'} + \mathrm {H H'}$, we search for a square saddle point, and with this constraint find a barrier of 132 kcal/mol, comparable to the 147 kcal/mol MRD-CI/exact barrier [21]. Thus to make the two hydrogens react, one needs to exert more energy than is needed to fragment the hydrogen molecules into its constituent atoms. This makes the reaction a forbidden one.

Hydrocarbons and diamond

eFF can describe a variety of bridged, fused-cyclic, and strained carbon skeletons, with largely correct carbon-carbon distances (Figure 10). The worst discrepancies in bond distances involve quaternary carbons in compounds like $ \mathrm{^{t}Bu-^{t}Bu}$ (1.708 $ \mathrm{\AA}$ vs 1.592 $ \mathrm{\AA}$ DFT) and diamond (1.681 $ \mathrm{\AA}$ vs 1.545 $ \mathrm{\AA}$ DFT).

Figure 10: eFF structures of several hydrocarbons and diamond, comparing eFF bond lengths with B3LYP/6-311g** values in parentheses.

Usually, eFF electrons separate into into core electrons that are centered on nuclei with $ Z>2$, and valence electrons that are located between nuclei. The valence electrons arrange themselves into nearest neighbor packings with a maximum of four electron pairs around each core (``octet rule''). When carbon has a full octet of electrons, they arrange themselves into a tetrahedral $ sp^{3}$ packing.

For diamond and hydrocarbons, eFF leads to a shell structure similar to density functional theory (Figure 11). The valence electrons in CC bonds are localized at the bond center, while the valence electrons in CH bonds are shifted off-center ($ \sim$0.62 bohr) toward the hydrogen, properly reflecting their attraction to the bare proton and their repulsion from the carbon 1s core electrons.

Figure 11: Comparison of electron densities between eFF and density functional theory (B3LYP/6-311g**) along a (a) carbon-hydrogen bond in methane, and a (b) carbon-carbon bond in ethane. In ethane the higher eFF density at the midpoint is confined to a small volume; the overall charge difference in the bonding region is only 0.13 electrons.

Over a test set of small hydrocarbons, valence electron ionization energies are within 0.1 eV of experimentally corrected [22] localized Hartree-Fock values for CC bonds (eFF 16.8 eV vs 16.7 eV) but are 2 eV too low for CH bonds (eFF 13.9 eV vs 16.0 eV). Core electrons are underbound (eFF 236 eV vs 290 eV for C(1s)) because electrons are represented by Gaussian functions rather than the optimal Slater form near the nucleus.

In eFF, the energy required to cleave a carbon-carbon bond into two radicals is too high, e.g. D( $ \mathrm{H_{3}C-CH_{3}}$) = 163.5 kcal/mol versus 93.4 kcal/mol DFT, and D( $ \mathrm{(CH_{3})_{3}C-CH_{3}}$) = 121.4 kcal/mol versus 86.0 kcal/mol DFT.

This error arises from differences in how well the eFF wavefunctions represent the true electron density in the molecule versus separated fragments. Carbon-carbon sigma bonds have a valence electron density that is concentrated in the region between the nuclei; hence the eFF bond-centered representation is a good one. On the other hand, methyl radical is poorly represented and has a too-high energy because eFF does not have the proper $ p$ functions to describe the radical electron. Hence ethane is overbound.

The eFF2 force field, which we outline later, provides an improved description of $ p$-like electrons. With eFF2, the ethane bond dissociation energy is 140 kcal/mol. However the eFF2 method is not currently as universally applicable as the original eFF method.

Steric interactions in hydrocarbons

eFF does well at capturing energy differences between multiple conformers of the same hydrocarbon. In traditional force fields, these interactions are handled via separate angle, dihedral, and noncovalent interaction terms, but in eFF they arise instead from a proper consideration of electrostatics and the Pauli principle.

Table 1 shows energy differences between staggered versus eclipsed ethane (CH/CH bond repulsion); gauche versus trans butane (CC/CC bond repulsion); chair versus boat cyclohexane (multiple gauche butane interactions); and axial versus equatorial methyl substituents on cyclohexane (diaxial interactions).

Table 1: Energy differences between conformers examined. $ \mathrm {^{*}}$Gauche butane is not a local minimum, and is constrained at $ \mathrm {60^{o}}$.
      $ \Delta E$ (kcal/mol)
system energy of relative to eFF DFT
ethane eclipsed staggered 2.1 2.7
butane gauche trans 1.6* 0.9
cyclohexane twist-boat chair 4.7 6.3
1,3-dimet.-cyclhex. ax-ax eq-eq 5.8 5.9
  ax-eq eq-eq 2.7 2.1

The geometry of cyclic alkanes is well-described (Figure 12), e.g. the curved bonds in cyclopropane (angle between bonding electrons of $ \mathrm{98^{o}}$ versus $ \mathrm{110^{o}}$ GVB [23]), the ring pucker of cyclopentane (dihedral $ \mathrm{21.5^{o}}$ versus $ \mathrm{33.2^{o}}$ DFT), and the twist-boat versus chair forms of cyclohexane (dihedral $ \mathrm{57.7^{o}}$ versus $ \mathrm{56.6^{o}}$ DFT for the chair, and $ 34.0^{o}$ versus $ 32.3^{o}$ DFT for the twist-boat).

Figure 12: Cyclic alkanes, comparison of eFF and B3LYP/6-311g** (parentheses) bond lengths in Angstroms and dihedral angles in degrees.

These results suggest that eFF (1) accurately describes the magnitude of steric repulsions between electrons on different atoms and (2) has the correct energetics for bending valence electrons away from a tetrahedral arrangement.

Ionic and multicenter bonds

The elements lithium, beryllium, and boron often participate in ionic bonds, since they are electropositive, and/or electron deficient multicenter bonds, since they lack enough electrons to complete a full octet. eFF describes lithium hydride, beryllium dihydride, and diborane well (Figure 13).

Figure 13: Lithium, beryllium and boron hydrides containing ionic and/or electron-deficient multicenter bonds. Bond lengths are in Angstroms, DFT values are in parentheses.

Lithium hydride has a bond length that is slightly too long (1.69 $ \mathrm{\AA}$ eFF versus 1.59 $ \mathrm{\AA}$ DFT). The energy needed to dissociate it into lithium and hydrogen atoms is very near the DFT value (58.1 kcal/mol versus 58.2 kcal/mol DFT).

Beryllium dihydride is linear with a bond length of 1.41 $ \mathrm{\AA}$ eFF versus 1.33 $ \mathrm{\AA}$ DFT. The energy needed to break $ \mathrm{BeH_{2}}$ into BeH and H is near the DFT value (113.0 kcal/mol versus 97.5 kcal/mol DFT); however, the energy of breaking BeH into Be and H atoms is too high (109.6 kcal/mol versus 57.7 kcal/mol DFT).

The boron compound $ \mathrm{BH_{3}}$ dimerizes into the borane $ \mathrm{B_{2}H_{6}}$ via two three-center two-electron bonds. eFF describes nearly all aspects of the $ \mathrm{BH_{3}}$ and $ \mathrm{B_{2}H_{6}}$ geometries correctly (see Figure 13). However, the dimerization energy is too low (27.6 kcal/mol versus 39.0 kcal/mol exact).

Uniform electron gas

eFF can be judged by its ability to reproduce the energetics of a uniform electron gas, the limiting case of electrons in a uniform neutralizing background charge. This system is defined by a single parameter, the Wigner-Seitz radius $ r_{s}$, related to the electron density $ \rho = (\frac{4}{3} \pi r_{s}^{3})^{-1}$.

eFF describes the uniform electron gas as a close-packed lattice of localized electrons analogous to a Wigner crystal packing [24], only at much higher electron density. Several lattices were considered and found to be similar in energy ($ <0.006$ hartree/electron), suggesting that the uniform electron gas in eFF has a fluxional character.

Figure 14: Potential energy of a uniform electron gas with respect to the density parameter $ r_{s}$.

Figure 14 shows the energetics of the eFF NaCl-packed uniform electron gas compared to Hartree-Fock energies and exact energies taken from accurate quantum Monte Carlo calculations [25,26]. We find that the well depth of the eFF and Hartree-Fock curves are similar ($ \sim$0.05 hartree/electron), but that the eFF curve is stretched out so that the minimum energy occurs at $ r_{s}\sim6$ bohr compared to 4.8 bohr for Hartree-Fock.

Lithium and beryllium solids

In eFF, metallic electrons are represented by diffuse electrons that nestle into the interstices between the ions. Interstitial sites in lithium are occupied by single electrons with random up or down spins (Figure 15), while analogous sites in beryllium are occupied by electron pairs (Figure 16). In the case of beryllium, the electron density distribution predicted by eFF agrees with those observed by X-ray crystallography [27] and in two separate periodic DFT studies [28,29].

Figure 15: Bonding in lithium bulk solid and its equation of state, DFT curve from [30]

Figure 16: Bonding in beryllium bulk solid and its equation of state.

eFF geometries and elastic constants of the two solids agree well with experimental values (Table 2), but cohesive energies are too high, probably because the $ 2s$ electrons of Li and Be atoms are not well represented in eFF.

Table 2: Lithium and beryllium parameters, comparison of eFF and experimental values.
    eFF expt
Lithium Lattice constant ( $ \mathrm{\AA}$) 4.42 4.40
  Bulk modulus 12.2 13.0
  Cohesive energy (kcal/mol) 60 38
Beryllium Lattice constant a ( $ \mathrm{\AA}$) 2.43 2.29
  Lattice constant c ( $ \mathrm{\AA}$) 3.72 3.59
  Bulk modulus 121.9 110-127
  Poisson ratio 0.092 0.032
  Cohesive energy (kcal/mol) 138 77

There is precedent to describing metallic bonding using localized electrons, i.e. the interstitial electron model developed by Li and Goddard [31]. In this model, valence electrons sit between crystal lattice sites, and interact as classical particles via effective potentials. Phonon dispersion relations in a variety of fcc metals were predicted and shown to agree well with experiment.

Dendritic lithium

Lithium metal is widely used as an electrode material for rechargeable batteries, due to its high electropositivity and low weight. During battery operation, lithium at the negative electrode oxidizes to become $ \mathrm{Li^{+}}$, which dissolves into the electrolyte. This process is reversed during the recharge cycle, but the replated lithium adds unevenly to the electrode surface, and as the battery is charged and discharged, dendrites grow from the surface, which can cause short-circuiting in the battery and lead to explosions.

Dendritic growth can be limited by encasing the lithium in large heterogeneous structures, such as intercalating solids or gels. To evaluate such structures and make them lighter and safer, it is useful to have a model that can describe both bulk and large-scale dendritic forms of lithium [32]. We demonstrate the use of eFF for this task.

For small lithium clusters, ab initio generalized valence bond (GVB) studies [33] show that as the number of atoms increases, the clusters with the most stability change from linear and zig-zag chains ($ N < 8$) to planar and three-dimensional clusters containing tetrahedral packings ( $ 8 \leq N \leq 20$) to bulk solids.

Figure 17 shows several eFF $ \mathrm{Li_{32}}$ structures obtained by minimizing different initial conditions, with a similar positioning of interstitial electrons as found in the GVB calculations.

Figure 17: Nonperiodic lithium, representative local minima with energies per atom given relative to periodic structure.

Figure 18 shows a larger heterogenous periodic structure obtained by minimizing a lithium bulk solid that was initially expanded twice along every dimension.

Figure 18: Electrons in dendritic lithium, showing interstitial electrons occupying spaces in between linear chains, two-dimensional meshes, and three-dimensional bulk-like structures. Electrons are shown reduced to 0.2 times their original size.