Ab Initio Molecular Dynamics of Temporary Anions Using Complex Absorbing Potentials.

Dissociative electron attachment, that is, the cleavage of chemical bonds induced by low-energy electrons, is difficult to model with standard quantum-chemical methods because the involved anions are not bound but subject to autodetachment. We present here a new computational development for simulating the dynamics of temporary anions on complex-valued potential energy surfaces. The imaginary part of these surfaces describes electron loss, whereas the gradient of the real part represents the force on the nuclei. In our method, the forces are computed analytically based on Hartree-Fock theory with a complex absorbing potential. Ab initio molecular dynamics simulations for the temporary anions of dinitrogen, ethylene, chloroethane, and the five mono- to tetrachlorinated ethylenes show qualitative agreement with experiments and offer mechanistic insights into dissociative electron attachments. The results also demonstrate how our method evenhandedly deals with molecules that may undergo dissociation upon electron attachment and those which only undergo autodetachment.

Temporary anions (TAs)[1,2] are formed when a neutral molecule captures a free electron, most often one with low energy below 10 eV. Although their lifetime is in the range of femto- to milliseconds, TAs play a central role in our understanding of important chemical processes: these include DNA lesions in living cells upon exposure to ionizing radiation,[3,4] the formation of radicals and anions in the Earth’s atmosphere,[5,6] and plasma etching in the semiconductor industry.[7] TAs are also currently reshaping our understanding of the physics and chemistry in interstellar and circumstellar environments.[8] Studying TAs by means of electronic-structure methods designed for bound states, that is, states related to the discrete spectrum of the electronic Hamiltonian, is difficult and, in many cases, a hopeless venture. To begin with, temporary anions are electronic resonances, i.e., metastable states embedded in the continuum of the electronic Hamiltonian.[9−11] Second, the formation of a TA may lead to fragmentation of the molecule, which is commonly referred to as dissociative electron attachment (DEA). If DEA is possible, it competes against autodetachment, which makes the computational modeling and prediction of the dissociation channels an intricate problem in its own right.[2,12−14]

Various methods have been proposed to address the difficulties associated with the modeling of TAs and DEA; recent overviews are available in refs (10 and 14). A group of methods combine scattering theory and operator projection techniques,[15−17] for example, the approach by Domcke[18] or the Schwinger multichannel approach by Takatsuka and McKoy.[19] In these methods, an explicit treatment of the continuum is necessary. A very accurate description of TAs and DEA is offered by nonlocal resonance theory,[18,20−22] but this entails high computational cost so that polyatomic molecules must be treated in reduced dimensionality considering only a few vibrational degrees of freedom.[23] A local approximation is thus usually invoked, which is valid as long as one does not deal with very broad resonances or effects near the threshold.[18,20] However, if the potential energy surfaces (PES) are constructed by single-point energy calculations and the nuclei are treated quantum mechanically, one is still limited to very few degrees of freedom.[24−29]

It appears, thus, worthwhile to apply ab initio molecular dynamics (AIMD) simulations to TAs, and we would like to mention here the recent work by Kossoski and co-workers.[30] They used bound-state methods to compute the resonance energies and determined the corresponding lifetimes using the Schwinger multichannel method.

In this Letter, we present a new computational development, dubbed CAP-AIMD, which integrates complex absorbing potentials (CAPs)[31] into AIMD simulations.[32] In our approach, we avoid scattering calculations and obtain the energy and lifetime of a TA as well as the forces acting on the nuclei in a single computation at each time step. Our approach invokes the local approximation and may not be suited for very broad resonances and near-threshold effects. In addition, we treat the nuclei classically, which represents a further approximation, but, together with analytic gradients, reduces the computational cost substantially so that the multidimensional PES of a polyatomic molecule can be treated in full dimensionality.

As initial applications, we report CAP-AIMD results for the temporary anions of dinitrogen, ethylene, chloroethylene, the three isomers of dichloroethylene, trichloroethylene, tetrachloroethylene, and chloroethane. It is well established that the chlorinated compounds all undergo DEA, while dinitrogen and ethylene do not;[14] the CAP-AIMD approach handles both cases evenhandedly.

The fundamental difference between conventional AIMD and CAP-AIMD is that in the latter, the nuclei evolve on a complex PES (CPES). The CPES is obtained from a non-Hermitian Hamiltonian with complex eigenvaluescalled Siegert energies.[33] Here, ER is the energy of a resonance state, and Γ is its decay width, which is inversely proportional to the lifetime τ. For states, which are electronically bound, Γ = 0.

In our approach, the CPES is computed on the fly using the CAP method. The CAP Hamiltonian Hη is obtained from the physical Hamiltonian H as[31,34−38]where η is a real scalar parameter and W is the CAP. We choose W to be quadratic in the electronic coordinate. Mathematical and technical details of this CAP, which has historically been the most popular one, have been discussed elsewhere.[31,35,36,39] Other possible choices include the Voronoi CAP[40] and reflection-free CAPs.[41−43]

The important thing to note in the context of the present work is that W depends on a set of three parameters, , which define its onset. Hence, the CAP is parametrized by the quartet , which needs to be optimized in order to minimize the perturbation of the resonance wave function. In principle, this needs to be redone at each time step as the nuclear configuration changes. We have found, however, that determining the CAP parameters for the initial structure of the TA and keeping them fixed during the AIMD simulation until the anion becomes bound produces meaningful results. When the anion is bound, η is set to zero. Details on the procedure for determining the optimal quartet are discussed in the Supporting Information.

Although it is possible to generate the CPES using other complex-variable techniques,[9,10] we chose here the CAP method with a quadratic potential W because of the availability of analytic gradients.[44,45] Similarly, while one can compute CPES at a higher level of electronic-structure theory, we limit ourselves to Hartree–Fock (HF) theory here because of the computational demands of AIMD simulations. Density functional theory (DFT) would be an attractive alternative but has only been implemented within the local density approximation for complex-variable techniques.[46,47] Also, analytic gradients are not available for complex-variable DFT.

The low level of HF theory notwithstanding, our CAP-AIMD results show qualitative agreement with experiments and offer insights into the connection between the initial TA and the final DEA products. Questions like, for example, how the molecular orbitals (MOs) of the TA evolve in time may also be answered.

We propagate the nuclei on the CPES according to the classical Hamilton equationswhere R, P, and M are the coordinate, momentum vector, and mass of the kth nucleus, and Vnuc–nuc is the nuclear repulsion energy. Note that the force F depends on only the real part of E. Indeed, if the nuclei are treated classically and nuclear quantum effects are ignored, the force F is simply proportional to the gradient of the real part of the CPES.[48]

After running a CAP-AIMD simulation, one obtains a collection {E} of complex energies, one for each time step n. The imaginary part of E relates to the resonance width Γ through eq . The {Γ} profiles can be used to distinguish TAs which undergo only autodetachment from those which undergo dissociation. In the first case, Γ oscillates with time, while in the second case, it drops to zero at a critical point and stays so indefinitely.

Furthermore, one can estimate from Γ the probability P that the TA survives autodetachment between steps t and t. Taking Γ as the width in the interval t ≤ t < t leads to the expressionwhere Δt = t – t. This probability may be employed to estimate the “lost fraction”, f, that is, the fraction of an ensemble of TAs, all characterized by exactly the same initial conditions, statistically expected to undergo autodetachment by the time step t. In the Supporting Information, we show that

It must be noted that the summation in the exponent in eq is an approximation to an integral of Γ over time in the continuous time limit. Hence, for a good estimate of f, one needs to choose the time step Δt such that Γ0Δt < 1, where Γ0 is the resonance width at the TA’s initial geometry.[49] As we will demonstrate in the following, our computed f values are in agreement with trends in experimental DEA cross sections. Eq thus provides a way to connect the AIMD simulations to experimental observables.

For meaningful results, it is paramount to conduct the CAP-AIMD simulation with the desired HF state, meaning the self-consistent field (SCF) solution that corresponds to the resonance and not some discretized continuum state. Indeed, as we show in the Supporting Information, the dynamics of pseudocontinuum states are fundamentally different from those of a resonance. We use the following procedure to find the right SCF solution in the continuum: (i) at time step n = 0, we build the SCF guess from the core Hamiltonian; then, (ii) at any subsequent time step n > 0, we use the SCF solution found at the (n – 1)th step as guess. As indicated above, the CAP is turned off as soon as the anion is bound. In order to determine this, we use Koopmans’ theorem and consider the real part of the energy of the highest occupied MO (HOMO): If it is negative, we consider the anion to be bound.

To judge how reliable a CAP-AIMD simulation is, we compute at every time step n a deperturbative correction to the complex energy according to[31,39]

For a well-represented resonance, and E deviate only little from each other, whereas this is not the case for pseudocontinuum states. Thus, in general, the closer E and are, the better the simulation. Note that the deperturbed lost fraction, , may also be computed by replacing Γ in eq with .

We implemented our method in the Q-Chem program[50] making use of the implementation of CAP-HF energies[36] and analytic gradients.[44,45] The results reported below should be viewed as proof of concept and as illustration of the robustness of CAP-AIMD simulations. For these reasons, we performed no averaging over initial structures and always took the optimized geometry of the neutral molecule, determined at the HF level, as the initial geometry of the corresponding TA, unless stated otherwise. This is equivalent to a vertical electron attachment at the neutral equilibrium geometry and simulating the time evolution of the formed TA. The cc-pVTZ + 3p basis set with the diffuse p-functions placed on all atoms except hydrogen was used for all simulations and the initial geometry optimizations. Only for tetrachloroethylene, we used aug-cc-pVDZ + 3p instead. For the TAs of dinitrogen, ethylene, chloroethylene, and the three isomers of dichloroethylene, we ran about 200 trajectories each, which enables a statistical analysis. For chloroethane, trichloroethylene, and tetrachloroethylene, we ran 67, 24, and 21 trajectories, respectively. All simulations were done in the microcanonical ensemble.

The 2Π shape resonance of N2– is an example of a TA which only undergoes autodetachment but does not induce fragmentation.[51] A similar case is the TA of ethylene, for which we report results in the Supporting Information. In Figure , we report resonance widths Γ and lost fractions f for N2– from randomly chosen CAP-AIMD trajectories. In Figure , we show the real part of the CPES derived from the same simulations. Cases A and B in Figure refer to simulations started at the same initial N–N bond length of 1.067 Å with the same CAP box dimensions but different η values: η = 0.02490 au, which is the optimal η for the initial geometry, in case A and η = 0.00500 au in case B. In case C, we started the simulation at an initial N–N bond length of 1.400 Å, where N2– is bound (see Figure ) and chose η = 0.00500 au for the unbound region.

Figure 1

Lost fraction f and resonance width Γ, as well as their deperturbed counterparts and for randomly chosen CAP-AIMD trajectories of N2– simulated under different conditions. The profiles are limited to the first 40 fs after vertical electron attachment but Γ and remain periodic in the entire duration of the simulations, which is ∼140 fs. Time step is 2 au ≈0.05 fs.

Figure 2

Real part of CPES of N2– (in green) and the PES of N2 (in magenta). Time step is 2 au ≈ 0.05 fs, total number of time steps is 2000. The curve for N2 derives from a standard AIMD simulation, while the curve for N2– is from a CAP-AIMD simulation.

The resonance width has a cyclic profile in all three cases in Figure because of the vibration of the molecule. We observe a direct correlation between Γ and the N–N bond length. Segments of the Γ profile in Figure where we see a decrease (increase) correspond to stretching (shortening) of the bond below 1.20 Å, whereas segments with Γ = 0 correspond to RN–N ≳ 1.20 Å. When Γ is nonzero, we also see an increase in the lost fraction, whereas it remains constant while Γ = 0. This suggests that N2– is unbound when RN–N ≲ 1.20 Å and that the potential curves of N2 and N2– cross at RN–N ≈ 1.20 Å according to Koopmans’ theorem. But Figure shows that the potential curves cross at RN–N ≈ 1.33 Å if independent HF calculations are performed for N2 and N2–; in fact, it is known that methods that describe a TA and the neutral molecule with the same Hamiltonian yield a more consistent description.[52] (See Supporting Information, Figure 10, for a plot of the imaginary part of the CPES of N2– against N–N bond length.)

The similarity between cases A and B, which differ only in the η value, indicates some flexibility in choosing this parameter. The extent of such flexibility is, however, dependent on the molecule and basis set; this is illustrated in the Supporting Information (Figures 17 and 18 and the corresponding discussion). Overall, we have better agreement between the uncorrected Γ (in red) and the deperturbed (in dark gray) in case B. Also, the discontinuity at Γ → 0 is smaller. However, there is good agreement between the uncorrected (f) and corrected lost fractions in both cases A and B. In case C, where we start at RN–N = 1.40 Å, the lost fraction curve is different from that in cases A and B; Γ stays at zero for the first ∼4.2 fs. In that initial interval, the N–N bond compresses to about 1.2 Å. Thereafter, as the N–N bond continues to compress, the autodetachment process ensues and Γ begins to increase. Because the nuclei gain more kinetic energy in case C than they do in cases A and B, more compressed bond lengths and higher Γ values are reached. Importantly, in all three cases in Figure , f and become ∼1 within 40 fs meaning that it will be very rare to find an anion 40 fs after electron attachment, which is less than 3 vibration periods of neutral N2.

We now turn our attention to molecules which may undergo DEA. In Figure , we show results from arbitrarily chosen trajectories for the anions of chloroethylene and the three isomers of dichloroethylene. Immediately, we see a stark difference with N2–: Γ does not oscillate but falls to zero in less than ∼10 fs, during which time we also see a steep rise and stabilization in the lost fraction profile.

Figure 3

Lost fraction f and resonance width Γ, as well as their deperturbed counterparts and for randomly chosen CAP-AIMD trajectories of the anions of (A) chloroethylene, (B) cis-, (C) trans-, and (D) 1,1-dichloroethylene. The profiles are limited to the first 40 fs.

Dissociation is possible in all four cases of Figure . The predominant dissociation channel we observe is the formation of C1– and an organic radical (see Supporting Information), which is in agreement with experiments.[53] It is worth noting that all TAs become bound before dissociation takes place. The latter is marked by a stabilization of the anion’s energy, which becomes discernible when one averages ReE over many trajectories (see Supporting Information). In our simulations, the C–Cl bond cleavage is complete 100–150 fs after electron attachment.

The important differences among the molecules presented in Figure are the lost fraction profiles, which are determined by two factors (see eq ): the resonance width and the time it takes a TA to become bound. Both factors relate to the shape of the CPES. In general, the closer the initial geometry of the TA is to a bound region on the CPES, the smaller the initial resonance width Γ0. Also, the more the real part of the gradient at the initial geometry points toward the bound region, the lesser time it will take the TA to reach that region.

Since the bound region is reached very quickly in all four simulations in Figure , Γ0 is the most important factor in explaining the lost fraction profiles. For the TAs of interest here, Γ0 decreases in the following order: CCl2 = CH2 > CHCl = CH2 > trans-CHCl = CHCl > cis-CHCl = CHCl. The equilibrium lost fractions f ≡ limf decrease indeed in the same order. The statistical analysis of our simulations (see Supporting Information) yielded for values of 0.93, 0.92, 0.33, and 0.10 for CCl2=CH2, CHCl=CH2, trans-CHCl=CHCl, and cis-CHCl=CHCl, respectively.

These values allow us to estimate how effective DEA is. For example, for 1,1-dichloroethylene, an average of 7% of the initially formed TAs live long enough for dissociation to happen. For trans- and cis-dichloroethylene, this average jumps to 67% and 90%, respectively, resulting in a much larger Cl– ion yield. This agrees well with experiments,[53] where it was found that the DEA cross section for the Cl– channel is highest for cis-CHCl=CHCl, followed by trans-CHCl=CHCl, and CH2=CCl2.

The ethylene-derived TAs under discussion here are known to originate from the capture of an electron into the π* orbital of the C=C bond, giving rise to a 2Π anion state.[54,55] According to electron transmission spectroscopy (ETS), their DEAs proceed, however, via a 2Σ state, which suggests a 2Π → 2Σ transition mediated by an out-of-plane motion during the C–Cl bond elongation.[54−57] As discussed in the Supporting Information, our CAP-AIMD simulations confirm that.

The 2Π → 2Σ transition is also evident from the changing character of the HOMO along the trajectory; Figure illustrates this for the chloroethylene anion. These orbital plots were generated along the trajectory from Figure A. It is seen that the HOMO is of π* character at t = 0, while we have mixed π*−σ* character at 48 fs. Notably, this plot is very similar to that of the coupled-cluster Dyson orbital computed at the minimum-energy crossing point between neutral and anionic chloroethylene.[57] Already at 72 fs, the HOMO is largely localized on the Cl atom and the Mulliken charge of this atom is already about −0.9 au

Figure 4

Time evolution of the HOMO of chloroethylene anion. Isovalue is 0.002.

A full description of the 2Π → 2Σ transition would require the consideration of nonadiabatic effects, which are neglected in the present work. However, we would like to point out that crossings between two anionic resonances, so-called exceptional points, have been analyzed based on a two-dimensional model[58,59] as well as using complex-variable electronic-structure calculations.[59,60]

In the Supporting Information, we also report CAP-AIMD results for the anions of trichloroethylene and tetrachloroethylene. Here too, we observe that the 2Π → 2Σ transition is mediated by an out-of-plane motion. For trichloroethylene, we observe three dissociation channels: the three isomers of dichloroethylene together with Cl–, and Cl + C2Cl– + HCl. All four channels have been observed in catalyzed reductive dechlorination;[61−63] however, while formation of the trans-dichloroethylene radical dominates in our simulations, all the experiments find cis-dichloroethylene radical to be dominant. This discrepancy may be explained by equilibration between the two radicals which is known to happen but it could also be related to stereo- and regioselectivity of the catalysts (e.g., vitamin B12) used in the experiments.[64] We note that our simulations do not take into account any environment and are performed at a low level of theory, i.e., HF. For tetrachloroethylene anion, we observe only the formation of trichloroethylene radical and Cl–, which is consistent with experiments.[61−63] The experiments also find subsequent dechlorination of the trichloroethylene anion formed from the radical yielding the same DEA products discussed in the last paragraph.

To test the CAP-AIMD method for DEA involving only a σ* resonance, we studied the TA of chloroethane. This is also reported in the Supporting Information; we find an equilibrium lost fraction very close to 1 indicating a small DEA cross-section consistent with experiment[65] and previous theoretical results.[30]

It is worth noting that we did not observe dissociation without a CAP. Such simulations do not capture electron attachment but describe pseudocontinuum states. The notable exception is tetrachloroethylene, whose anion is almost bound at the neutral equilibrium structure making it possible to simulate DEA without a CAP. As shown in the Supporting Information, dissociation is also suppressed when too high or too low CAP strengths are used. Finally, we note that we observed in a small number of trajectories artificial imaginary energies long after the anion has become stable against electron loss. We ignored these artificial imaginary energies in the analysis of the results.

In summary, CAP-AIMD offers a robust method to simulate the dynamics of TAs, independent of whether they undergo dissociation or not. The nuclei are propagated on a complex potential energy surface, yielding a complex energy at each time step, where the imaginary part relates to the decay width of the TA. In addition, the fraction f of TAs lost due to autodetachment can be obtained from CAP-AIMD simulations. The equilibrium value provides a measure for the efficiency of DEA. We conducted CAP-AIMD simulations for the anions of chloroethane and different chlorinated ethylenes; the observed trends in the DEA yield are consistent with experimental results. A remarkable commonality among these simulations is that all anions become stable toward autodetachment after less than 10 fs. We conclude by noting that an obvious challenge for future work is to improve the description of the electronic structure, that is, replacing HF by DFT in CAP-AIMD simulations. The inclusion of nonadiabatic effects may be a further task worth pursuing.