Ab Initio Molecular Dynamics Simulations of Solvated Electrons in Ammonia Clusters.

We investigated excess electron solvation dynamics in (NH3)n- ammonia clusters in the n = 8-32 size range by performing finite temperature molecular dynamics simulations. In particular, we focused on three possible scenarios. The first case is designed to model electron attachment to small neutral ammonia clusters (n ≤ ∼10) that form hydrogen-bonded chains. The excess electron is bound to the clusters via dipole bound states, and persists with a VDE of ∼160 meV at 100 K for the n = 8 cluster. The coupled nuclear and electronic relaxation is fast (within ∼100 fs) and takes place predominantly by intermolecular librational motions and the intramolecular umbrella mode. The second scenario illustrates the mechanism of excess electron attachment to cold compact neutral clusters of medium size (8 ≤ n ≤ 32). The neutral clusters show increasing tendency with size to bind the excess electron on the surface of the clusters in weakly bound, diffuse, and highly delocalized states. Anionic relaxation trajectories launched from these initial states provide no indication for excess electron stabilization for sizes n < 24. Excess electrons are likely to autodetach from these clusters. The two largest investigated clusters (n = 24 and 32) can accommodate the excess electron in electronic states that are mainly localized on the surface, but they may be partly embedded in the cluster. In the last 500 fs of the simulated trajectories, the VDE of the solvated electron fluctuates around ∼200 meV for n = 24 and ∼500 meV for n = 32, consistent with the values extrapolated from the experimentally observed linear VDE-n-1/3 trend. In the third case, we prepared neutral ammonia cluster configurations, including an n = 48 cluster, that contain possible electron localization sites within the interior of the cluster. Excess electrons added to these clusters localize in cavities with high VDEs up to 1.9 eV for n = 48. The computed VDE values for larger clusters are considerably higher than the experimentally observed photoelectric threshold energy for the solvated electron in bulk ammonia (∼1.4 eV). Molecular dynamics simulations launched from these initial cavity states strongly indicate, however, that these cavity structures exist only for ∼200 fs. During the relaxation the electron leaves the cavity and becomes delocalized, while the cluster loses its initial compactness.

Introduction

The ammoniated electron is an archetypal example of solvated electrons. The first observations of the ammoniated electron are attributed to Sir Humpry Davy, and later, to Wolfgang Weyl.[1] This phenomenon, the appearance of a blue color after dropping a piece of sodium or potassium metal in liquid ammonia, is widely known and described in detail in standard university textbooks.[2] Originally, Kraus proposed that the color is due to the presence of ammoniated electrons,[3] conveniently viewed in this picture as free electrons stabilized (solvated) in liquid ammonia. Solvated electrons have also been observed in several other solvents. The bottleneck of the observation of the solvated electrons is their (usually very short) lifetimes. While the ammoniated electron persists for days, the solvated electron in water (hydrated electron) is a short living species surviving only several microseconds.[4] Since its experimental detection in 1962,[5] the hydrated electron has become the main target of the investigations, which is due to its paramount role as a highly reactive intermediate in aqueous phase chemical reactions of high practical and theoretical importance.[6] Nevertheless, the ammoniated electron has still remained a fascinating and challenging case of the solvated electron family.

The physical properties of the ammoniated electron has been thoroughly explored by experimental studies both in bulk liquid[7−11] and in finite size clusters.[12,13] The earliest theoretical studies relied on a one-electron dielectric continuum model.[14,15] This model led to a simple structural picture of the ammoniated, and in general, the solvated, electron: a localized excess electron distribution in a solvent void (cavity) surrounded by ammonia (solvent) molecules.[14] This is the cavity model. More sophisticated approaches that replace the solvent to classically treated molecular baths in one-electron path-integral and quantum molecular dynamics simulations brought about significant improvements in our understanding of the physical properties of the species.[16−19] The many-electron treatment of the ammoniated electron was first introduced in static calculations for finite size clusters at the Hartree–Fock (HF) level,[20,21] later using DFT[22,23] and ab initio calculations.[23,24] In the only many-electron based study that models the temporal behavior of the ammoniated electron, Chaban and Prezhdo simulated the dynamics of alkali and alkaline earth metals solvated in a bath of 32 ammonia molecules in a periodic simulation cell.[25]

Since the application of many electron quantum chemistry necessarily introduces severe limitations on the size of the investigated systems, the focus of interest shifted to finite-size ammonia cluster anions. Nevertheless, the investigation of smaller ammonia clusters raised several challenges, most notably, on the existence of the cavity structure. The cavity model of Ogg[14] has been consistent with, and therefore, largely supported by one-electron simulations,[16−19] and early many-electron calculations.[20,21] On the basis of DFT calculations on small clusters, Shkrob, however, suggested that the bulk ammoniated electron should be more properly considered as a multimer radical anion.[22] Later, high-level correlated ab initio calculations by Sommerfeld questioned this proposal.[24] Most recently, we also performed ab initio and DFT calculations on ammonia cluster anions, extending the investigated systems up to 32 monomers.[23] We found that long extended chains of ammonia cluster anions localize less than 1% of the excess electron density in the frontier orbitals of the nitrogen atoms, a minor effect that largely supports Sommerfeld’s view.[24] Although this argument seems to be valid for chains of dipole bound anions, one has to realize that physical properties observed for finite size systems are not necessarily transferable to the bulk.

The central challenge of the ammoniated electron structure is also strongly related to the size effect. It is still to be rationalized how the dipole-bound character of the excess electron of small cluster anions evolves with cluster size to an excess electron that is confined and solvated in the bulk of the ammonia liquid. Experimentally, photoelectron spectroscopy of ammonia cluster anions provides information on the vertical detachment energy (VDE) of the excess electron from the cluster anions. Sarkas and his co-workers measured the VDE of ammoniated electron clusters in the cluster size range of n = 41–1000 and found a linear VDE–n–1/3 relationship.[12]

We note here that varying the backing pressure of the carrier gas in the cluster preparation procedure may lead to different tendencies of the VDE with cluster size, as was observed for water and methanol clusters.[26] The different trends are usually associated with the existence of different electron binding motifs in the clusters. While there has been only one single linear tendency observed for ammonia clusters, at least three such patterns exist for water, and two exist for methanol.[26] Previous one-electron model based simulations[27−29] and ab initio molecular dynamics (AIMD) studies on water[30,31] and methanol cluster anions[32] identified two main excess electron localization modes: interior states, where the excess electron density is mainly localized within the cluster, and surface states, where the electron stabilizes mostly outside the cluster. The extra complexity in the aqueous case has been attributed to cases that lay between these two limiting geometrical arrangements.[29]

Until now, similar experiments have not been performed for ammonia cluster anions. The single observed linear relationship has been attributed to the presence of “embryonic” solvated electron states that gradually converge with size to the bulk ammoniated electron.[12] This picture opens the way to (at least) two alternative interpretations. It is either that the smallest observed clusters (n ∼ 40) already confine the electron in some kind of embedded states and the interior character is smoothly enhanced as the number of solvent molecules increase around the electron or that an initially nonembedded state (i.e., surface state) in smaller clusters develops an increasing interior character and progressively becomes an interior state in sufficiently large clusters.

In the present paper, we aim to investigate this problem by performing AIMD simulations. Previously, we reported a combined quantum mechanical and molecular dynamics study, a first step in this direction.[23] First, we examined the electron binding properties of optimized linear hydrogen-bonding chains along with several other branched configurations. These small clusters in the n = 3–8 range bind the electron weakly in dipole bound states with VDEs up to ∼200 meV. These clusters, anionic configurations at ∼0 K, represent minima on the potential energy surfaces that are unlikely to be accessible experimentally. We also generated intermediate size configurations with n = 8–32 monomers using a different sampling procedure, finite temperature classical and ab initio molecular dynamics.[23] In these simulations, neutral ammonia cluster equilibrium trajectories were generated, and selected configurations along the trajectories were tested for excess electron binding. We note that this simulation procedure can be thought of as a simplified model of the first step of an electron attachment scenario, where a zero kinetic energy electron is attached to a neutral cluster. The electron binding strength in these clusters appeared to be significantly weaker than that in water clusters[33] but is similar to that of methanol clusters of the same size.[34] Subsequent quantum mechanical calculations along the neutral trajectories indicated that initially the excess electron is always attached to the surface of the ammonia clusters for the investigated cluster sizes.[23]

In the present work, we wish to investigate the molecular dynamics of ammonia cluster anions. We launch different relaxation trajectories for various types of clusters. The investigated cluster sizes are selected from the n = 8–32 range. First, we investigate the dynamics of small, linear chains of ammonia cluster anions. The prototype of these simulations will be performed for the n = 8 cluster. The next scenario investigates the relaxation of cluster anions that are prepared from equilibrated neutral clusters by adding an extra electron to the cluster. The main questions we address here are how the initially surface-bound electrons relax on the surface and whether we can identify signs of internalization in the investigated clusters. The third group of simulations is launched from ammonia clusters that are prepared with initial (although artificial) binding sites for the extra electron in the interior of the clusters. Here, we aim to investigate the initial electron localization properties in such clusters and follow the fate of the original excess electron distribution. The results of the simulations will be compared to available experimental results.

The structure of the paper will be as follows. Section contains the description of the cluster preparation procedure of the investigated configurations, and provides the technical details of the molecular dynamics techniques employed in the present study. Section collects the results of the simulations, and analyses the resulting relaxation trajectories. In section we discuss the main messages of the study and conclude the paper.

Methods

We performed ab initio molecular dynamics simulations to model excess electron solvation in finite size ammonia clusters. Since the applied methods are similar to those employed in a previous work,[23] we outline the simulation techniques only briefly here, but provide more discussion on a few specific key details.

In the present work, we employ DFT with the BHandHLYP functional[35,36] for the electronic structure calculations in the MD simulations. This functional was successfully tested and employed in a series of studies for aqueous,[37−39] methanolated,[40] and ammoniated excess electron systems.[23] In particular, it was found[37,38] that that the higher Hartree–Fock (HF) exchange character in the BHandHLYP functional leads to more reliable description of the solvated electron systems relative to previously applied functionals, such as BLYP, B3LYP, and PBE.

In general, these findings are in line with the observations about the problems of the application of standard DFT methods to anions.[36,41] Functionals that do not include correction for the self-interaction error or for the incorrect long-range behavior are especially prone to perform poorly when applied to loosely bound electrons.[42] One possible way for the improvement is by increasing the HF exchange character of the functional, such as in BHandHLYP, that we examined in methanol[40] and ammonia cluster anions.[23] An alternative solution is to use range-separated functionals.[41,42] Our initial tests on ammonia cluster anions indicated that the long-range corrected (LC) VDEs shift unfavorably upon optimizing the range separation parameter with respect to the nonoptimized long-range corrected calculations (i.e., LC-BLYP) and the CCSD(T) numbers, a matter that needs further clarification.[23]

Although the most common DFT methods have been shown to fail to account for intermolecular dispersion interactions,[43] in the present study we choose to apply the BHandHLYP functional without the application of dispersion corrections. Our argument for this choice is based on test calculations of the interaction energy of the ammonia dimer (for details, see the Supporting Information). We found that the BHandHLYP method used with the basis sets of the present study (see below) gives a reasonable estimate of the ammonia–ammonia interaction energy that is mostly dominated by hydrogen-bonding and dipole–dipole interactions. Dispersion accounts for ∼10% of the interaction energy. We also find that the popular D3 dispersion correction[44] overestimates the interaction energy by ∼10%. Overall, we conclude that the method we applied in molecular dynamics simulations appears to be approximately as adequate for simulating ammonia clusters as the D3-corrected method. We anticipate that the present mild underestimation of the intermolecular interactions will not influence significantly the general trends of the dynamics of the solvated electron clusters.

For the representation of the electron density, we applied the Gaussian and plane wave method (GPW) in combination with the Goedecker–Teter–Hutter pseudopotentials.[45,46] We found previously that the GPW method with the MOLOPT/TZVP basis set provides an effective way to perform MD simulations on methanol cluster anions.[47] We also noticed, however, that this basis set lacks the diffuse character that is necessary to model relatively small size cluster anions. We remedied this problem by empirically adding Pople-type diffuse functions to the original MOLOPT/TZVP basis set to improve the computed VDE values. We apply a similar procedure here for the ammonia cluster anions.

The localized basis set of the simulations is similar to that used in our previous work to generate neutral ammonia cluster trajectories,[23] where we used the MOLOPT-TZVP basis[48] augmented on the hydrogen atoms by a diffuse function of the 6-31++G(d) basis set[49,50] and an additional s orbital with a coefficient that is one-third of the previous diffuse function. Since the excess electron distribution in the anionic ammonia clusters is presumably weakly bound and diffuse, here we apply a more diffuse basis set by augmenting the previous set with further extra functions, the diffuse functions of the 6-31(++)G(d) set for the nitrogen[50] and one additional extra s orbital on the hydrogen atoms with a coefficient that is one-third of the previous diffuse function of the hydrogen atoms. The extra diffuse functions that are added to the MOLOPT-TZVP set are listed in the Supporting Information. To reduce the cost of the evaluation of the HF exchange energy, we applied the auxiliary density matrix methods (ADMM)[51] with the diffuse aug-FIT3 basis set. The energy cutoff of the plane wave part of the GPW basis was set to 480 Ry in the simulations. The applied basis set will be referred to as GPW-BS1 in the text.

The simulations were performed in cubic simulation cells with box sizes between 25 and 40 Å with the application of the Martyna–Tuckerman Poisson solver[52] to treat long-range interactions of the clusters. The molecular dynamics simulations were carried out in the canonical (NVT) ensemble using a 0.5 fs time step at T = 40 and 100 K. The temperature was regulated by the canonical sampling through velocity rescaling (CSVR) thermostat[53] using a time constant of 100 fs.

The central quantity we focus on in the present study is the vertical detachment energy of the excess electron of the anions. The time evolution of the VDE is evaluated along the molecular dynamics trajectories using the atom-centered 6-31(1+3+)G(d) basis set[38] and the aug3-cc-pVDZ basis.[38,54,55] The justification for the use of these basis sets will be given in the next section.

In the present work, we employ the CP2K program package to perform the AIMD simulations.[56] The Gaussian software suite is also used for molecular dynamics, electronic structure calculations, and VDE calculations.[57] The molecular structures and spin densities were analyzed with the GaussView program of the Gaussian09 program package.[57]

Results

Molecular Dynamics of Dipole Bound Ammonia Cluster Anions

This section is devoted to studying the dynamics of octamer anion chains that represent the electron binding motif in the smallest size clusters, dipole bound anions. Here, at the outset, our goal is two-fold. First, we aim to model the behavior of the dipole bound anionic species per se, and, at the same time, test the applied basis sets and simulation tools and develop an error estimate for our later simulation results for more extended systems.

In order to assess and exclude the possible artifacts associated with the use of the combination of periodic basis functions (i.e., plane waves) and Poisson solvers in the simulations of finite size clusters, first we performed BHandHLYP-driven MD on linear octamer anions with the atom-centered 6-31(1+3+)G(d) basis set. The VDEs were computed using the 6-31(1+3+)G(d) and aug3-cc-pVDZ sets. The combinations of this method and basis sets were previously illustrated to give a reasonable estimate of the vertical detachment energies of water,[37] methanol,[40] and ammonia cluster anions.[23] In particular, we pointed out that BHandHLYP/aug3-cc-pVDZ and BHandHLYP/6-31(1+3+)G(d) VDE values for ammonia cluster anions are only ∼10 and ∼25 meV higher, respectively, than the corresponding CCSD(T) numbers in the 100–200 meV VDE range.[23] We carried out this simulation with the Gaussian program package.[57] For comparison, we also computed the VDEs for the BHandHLYP/6-31(1+3+)G(d)-generated configurations with the GPW-BS1 basis set using the CP2K package.[56]

We started the simulation from the BHandHLYP/6-31(1+3+)G(d)-optimized neutral ammonia chain that we gradually warmed up from 0 to 100 K during a 2 ps long trajectory (see also below). Once we reached the desired temperature, an extra electron was attached to the cluster, and a 2 ps long MD trajectory was launched to monitor the temporal evolution of the VDE. Figure shows the time evolution of the VDE along the trajectory using the 6-31(1+3+)G(d) basis set. The neutral cluster initially binds the electron weakly (∼100 meV). Nuclear relaxation of the cluster leads to the enhancement of the VDE up to ∼160 meV (157 ± 4 meV, p = 0.01 significance level) during the trajectory. This VDE range is consistent with (and only slightly lower than) the reasonably converged estimate of the VDE of the optimized octamer anion.[23] While the use of a more extended atom-centered basis set (aug3-cc-pVDZ) decreases the VDE values only slightly, by ∼20 meV (15 ± 3 meV, p = 0.01 significance level), the application of the GPW-BS1 set results in significantly lower VDE values by a constant downshift of ∼0.2 eV (190 ± 5 meV, p = 0.01 significance level, see also Figure S1). We further discuss the basis set issues below.

Figure 1

Time evolution of the VDE of an ammonia octamer anion following electron attachment (at t = 0 fs) to a neutral cluster at 100 K. The Born–Oppenheimer MD run was generated using the BHandHLYP functional with the 6-31(1+3+)G(d) basis set. The VDE values of the cluster anion configurations along the trajectory are computed with the 6-31(1+3+)G(d) basis set.

Interestingly, the VDE reaches the average relaxed VDE very quickly, within the first 100 fs (see Figure ). The inset of Figure also reveals that in this time regime the most important contribution of the VDE fluctuations takes place on the ∼30 fs characteristic time scale. This fact indicates that a ∼ 1000 cm–1 wavenumber nuclear motion, the umbrella mode of the ammonia molecules, is apparently strongly coupled to the energy levels of the excess electron. One may also notice a ∼ 10 fs periodic contribution that appears as a minor component in the energy relaxation and corresponds to the ∼3300 cm–1 N–H stretching vibrations. However, the most important component of the VDE fluctuations (∼300 fs) is related to the ∼100 cm–1 wavenumber librational modes associated with the intermolecular motions within the chain. These modes result that the hydrogen-bonded chains distort and get slightly kinked during the finite temperature simulation. The distortion of the chains does not affect the electron localization mode; the electron localizes to the free hydrogen bond donor end of the cluster and remains there during the dynamics. The spatial extent of the excess electron, faithfully reflected in the spin density isosurfaces, shrinks during the dynamics correlating well with the increasing VDE (Figure ). Nevertheless, the excess electron density remains very diffuse extending over and covering several ammonia units of the chain.

Figure 2

Mesh representation of the spin density isosurfaces of octamer anions generated in atom-centered basis set MD simulations at t = 0 fs (top figure, the excess electron attached to the neutral ammonia cluster) and after t = 1000 fs of relaxation. The surfaces contain 80% of the total spin density (isovalues of 0.000007 and 0.000010, respectively). The electron distributions are computed with the 6-31(1+3+)G(d) basis set.

As another basis set comparison, we performed MD simulations in similar spirit using the BHandHLYP functional with the GPW-BS1 basis set, and computed the VDE along these trajectories using the GPW-BS1 basis and the two atom-centered sets, 6-31(1+3+)G(d) and aug3-cc-pVDZ. Here, we also started from the geometry of the optimized neutral octamer (see above), heated up the structure to the two desired temperatures (40 and 100 K) without the excess electron (using AIMD simulations), attached an extra electron to the structures, and performed 1 ps long simulations at 40 and 100 K. We also performed similar simulations starting from the optimized anion at 0 K,[23] heating it up to 40 and 100 K, and running 1 ps long MD trajectories at these two temperatures. The subsequent temporal evolutions of the anions were monitored by computing the vertical detachment energies of the anions along the AIMD production trajectories. Figure S2 shows that the anion VDEs in all four simulations (computed with the GPW basis set) are around (or slightly below) zero, while single-point VDE calculations on the same configurations using the atom-centered 6-31(1+3+)G(d) and aug3-cc-pVDZ basis sets predict VDE values that fluctuate around the VDE of the optimized octamer anion chain, 150 meV. Thus, here we can also conclude that the GPW basis set MD simulations systematically underestimate the VDE values predicted by the atom-centered sets by ∼200 meV. Nevertheless, the positive VDE values of the atom-centered basis calculations suggest that the set of collected configurations along the GPW-generated MD trajectories represent stable structures. In addition, the atom-centered basis VDE values along the GPW trajectory are very similar to those obtained from the atom-centered basis set MD simulations.

Now it remains to be assumed, or better to be demonstrated, that the two simulations sample the same region of the configuration space. To this end, we compared the dipole moments of the neutral octamer cluster along the two molecular dynamics trajectories (Figures and S2a) computed with the 6-31(1+3+)G(d), the aug3-cc-pVDZ, and the GPW-BS1 basis sets. We obtained (at the p = 0.01 significance level) 12.97 ± 0.35 D and 12.53 ± 0.34 D in the atom-centered MD and the GPW-generated MD simulations with the 6-31(1+3+)G(d) set, respectively, 11.68 ± 0.52 D and 11.41 ± 0.71 D with the aug3-cc-pVDZ basis, and 12.61 ± 0.31 D and 12.16 ± 0.32 D with the GPW-BS1 basis set. Although the dipole moment significantly varies within the same trajectory depending on the applied basis, the two more reliable atom-centered basis sets predict statistically identical behavior (at the p = 0.01 significance level) in the two different MD trajectories. On the basis of this statistical similarity, we may assume in the remainder of this paper that the two simulation methods essentially sample the same set of configurations.

Since simulations using the GPW approach are much more feasible both in terms of computational resources, the maximum size of the investigated systems and the duration of the simulations than the atom-centered basis set simulations, for the larger clusters we will use only the GPW-BS1 basis AIMD simulations. We demonstrated that although these simulations are likely to underestimate the VDE of the clusters, they essentially sample the same set of configurations as the atom-centered 6-31(1+3+)G(d) or aug3-cc-pVDZ basis simulations (that provide more reliable estimate of the VDE). For this very reason, in later sections we generate the molecular dynamics trajectories with the GPW-BS1 basis set (AIMD simulations with the CP2K package)[56] and perform the VDE (and spin density) calculations with the atom-centered 6-31(1+3+)G(d) and, for testing purposes, aug3-cc-pVDZ basis sets (with the Gaussian program package).[57]

We are aware that these atom-centered basis sets may also be of limited accuracy for weakly bound anions.[42] Nevertheless, in our previous study for small ammonia cluster anions we observed the following:[23] (i) The aug3-cc-pVDZ set showed reasonably fast convergence (at ∼0.1 eV VDE) with respect to the number of diffuse shells in CCSD(T) calculations for increasing size ammonia cluster anions. This suggested that application of the aug3-cc-pVDZ set is of sufficient precision for the cluster sizes of the present study. (ii) BHandHLYP/aug3-cc-pVDZ and BHandHLYP/6-31(1+3+)G(d) VDE values are only ∼10 and ∼30 meV higher, respectively, than the CCSD(T) numbers in the 100–200 meV VDE range. This reasonably suggests that critically used BHandHLYP/aug3-cc-pVDZ and BHandHLYP/6-31(1+3+)G(d) VDE values can be useful at characterizing the energetics of ammonia cluster anions in the examined size (and VDE) range, that is, for n = 8–32 monomers.

Molecular Dynamics of Intermediate-Size Surface-State Ammonia Cluster Anions

Here we investigate the relaxation of excess electrons added to compact neutral (NH3) clusters in the n = 8–32 size range that were prepared in our previous study.[23] The procedure begins with the preparation of neutral ammonia clusters using the semiempirical AM1[58] driven molecular dynamics procedure of approximately 1 ns length (for details see ref (23)). We selected 10 neutral configurations from the last 100 ps equilibrium part of a T = 100 K equilibrium trajectory separated by 10 ps time segments for each cluster size. The 10 selected configurations were used as the starting points of additional 2.5 ps long AIMD simulations for the neutral ammonia clusters.[23]

First, we analyzed the electron-binding affinities of the trajectories along the last 2.5 ps long neutral trajectories. We added an excess electron to the neutral configurations at every 50 fs along one of the 10 trajectories and computed the VDE of the anion. The analysis of the VDE for these 51 configurations for each cluster size confirmed our previous finding[23] that the initial electron binding affinities of the clusters increase with increasing cluster size. The two largest examined clusters, the n = 24 and 32 clusters, bind the electron with average binding energies of 60 and 170 meV, respectively. Figure (and Figure S3) show the initial electron binding energies along the neutral cluster trajectories. The significant electron binding affinity of the neutral n = 32 clusters (Figure ) that reaches as high as 290 meV strongly suggests that neutral clusters of approximately at or above this size may effectively capture low kinetic energy free electrons in experiments.

Figure 3

Computed VDE values of (NH3)32– anions following electron attachment to neutral clusters (or, equivalently, electron affinities of (NH3)32 neutral clusters) sampled along the molecular dynamics trajectory of the neutral species at 100 K.[23] The VDE values of the cluster anion configurations are computed with the 6-31(1+3+)G(d) basis set.

After the collisions of the neutral clusters with slow electrons, several scenarios can take place. For the smallest clusters, the electrons will most likely be scattered or only temporarily attached to the clusters. These unstable species may suffer subsequent autodissociation to smaller clusters or/and the extra electron may be autodetached from the cluster. The increasing affinity to bind an excess electron in growing size neutral clusters (as reflected by the increasing initial electron binding energies with size) indicates that the formation of stable anionic species becomes more likely in larger ammonia clusters. In general, positive VDE values and localized spin densities of the newly formed anions are good indicators that electron attachment to the neutral clusters is likely to take place. Electron detachment, however, is signaled by nearly zero (or negative) VDE values and highly delocalized spin densities of the cluster anions. We note here that even in the case of unfavorable electron–cluster interaction the use of finite size basis sets and/or finite size simulation cells necessarily keeps the excess electron density in the vicinity of the molecular species. This artificial constraint may lead to negative VDE values in these unstable situations.

To examine the molecular level events ensuing an electron–neutral cluster encounter, we performed MD simulations of clusters that just captured an excess electron. In practice, we modeled electron attachment by simply switching the charge of selected neutral clusters along their neutral trajectories from 0 to −1. We then launched molecular dynamics simulations starting from these now negatively charged ammonia clusters. Thus, the initial cluster anions of these simulations can be thought of as the result of a hypothetical encounter of a neutral cluster and a zero kinetic energy electron. We note that similar routes were followed in previous computer experiments, MD simulations of electron attachment to neutral methanol clusters.[32,47] One should be aware, however, that this type of computer experiments represent a simplified scenario of the electron attachment and the subsequent relaxation. In real experiments, cluster anions are likely to be formed and stabilized by evaporative cooling.[59] Nevertheless, we anticipate that the present simulations provide important information on the elementary microscopic details of the relaxation mechanisms.

We performed 10 MD simulations for each of the n = 8, 12, and 16 cluster sizes starting from the last configurations of their corresponding neutral MD trajectories. For the n = 8 and 12 clusters, we did not observe the formation of stable cluster anions. Clusters of this size tend to disintegrate with highly delocalized spin densities. Similar general observations can be made for most of the n = 16 trajectories. Nevertheless, the trajectory launched from the configuration that binds the electron most strongly (by 76 meV) of the initial geometries, shows definite signs of temporary stabilization. Figures S4 follows the time evolution of the VDE for this n = 16 cluster anion. The VDE climbs up temporarily to around 150 meV, before falling back to zero, then rising and falling again. Figure S5 depicts the cluster configurations and the spin density distributions at the beginning (t = 0 fs) and at the end (t = 1000 fs) of the dynamics. The cluster clearly disintegrates by the end of the 1 ps dynamics with a highly delocalized spin density, indicating the instability of the anion. This instability is also reflected in the tendency of the VDE values to fluctuate near zero.

The tendency of the stabilization during the dynamics following excess electron attachment to neutral clusters continues and culminates in the n = 24 and 32 trajectories. Here we launched MD simulation from only a single, randomly selected, neutral configuration for both cluster sizes. The VDE evolution of the n = 24 MD trajectory starts from a configuration that binds the electron initially with a VDE value of only 8 meV. The VDE rises to ∼400 meV within the first few hundred femtoseconds, then falls back to and stabilizes around 200 meV by the end of the trajectory (Figure , top). A similar, but more enhanced, VDE tendency can be traced for the n = 32 anion. Here the initial ∼100 meV VDE shoots up as high as 600 meV and ends up at 400 meV after 1 ps dynamics (Figure , bottom). In the last 500 fs of the two simulated trajectories, the VDE of the solvated electron fluctuates in the 210 ± 30 meV interval for n = 24 and around 500 ± 70 meV for n = 32 (p = 0.01 significance level).

Figure 4

Time evolution of the VDE of a) (NH3)24– and b) (NH3)32– ammonia cluster anions following electron attachment (at t = 0 fs) to neutral clusters at 100 K. The MD run was generated using the BHandHLYP functional with the GPW-BS1 basis set. The VDE values of the cluster anion configurations along the trajectory are computed with the 6-31(1+3+)G(d) basis set.

The temporal change of the excess electron distribution during the dynamics parallels the VDE trends. The spin density at the outset of the dynamics appears to be very diffuse essentially engulfing the cluster framework as illustrated in Figure (top). The huge enclosed volume containing 80% of the total spin density clearly suggests extensive delocalization and the lack of presence of any specific electron binding sites in the cluster. Nuclear relaxation quickly (within ∼200 fs) acts to accommodate the extra electron effectively. In the one-electron picture, the excess electron’s wave function localizes and becomes relatively strongly bound on the cluster. This localization is clearly manifest in the shrinking volume of the spin density distribution (Figure , bottom). Here, the nuclear configuration of the cluster also undergoes a significant transformation, losing its initial compactness and forming a distinct cleavage on the surface of the cluster that apparently acts as a binding site for the extra electron. We also emphasize that due to the available limited simulation time we are unable to rule out the possibility of further transformations of the surface clusters during the course of the dynamics. Nevertheless, the observation of surface type binding patterns on the ammonia clusters is indicative and provides us important information on the qualitative magnitude of the VDE for this type of cluster anions.

Figure 5

Mesh representation of the spin density isosurfaces of (NH3)32– generated in AIMD simulations at t = 0 fs (top figure, the excess electron attached to the neutral ammonia cluster) and after t = 1000 fs of relaxation (bottom figure). The surfaces contain 80% of the total spin density (isovalues of 0.000006 and 0.000014, respectively). The electron distributions are computed with the 6-31(1+3+)G(d) basis set.

Molecular Dynamics of Intermediate-Size Interior-State Ammonia Cluster Anions

At the end of the previous section, we observed the formation and stabilization of negatively charged ammonia clusters after the attachment of an extra electron to the neutral clusters with n = 24 and 32. These cluster anions initially bind the electron in weakly bound, diffuse states, relax within the 1 ps time frame of our simulations, and bind the extra electrons on the surface of the clusters in more localized states with VDE as high as 600 meV for the n = 32 size. One may, however, wonder whether the localization of the extra electron ends in a localized surface state or progresses further with possible penetration of the electron density in the interior of the cluster (internalization). Such internalization was observed in one-electron simulations of water cluster anions.[28] This process, however, is anticipated to take place on a time scale that is unavailable for our MD simulations. The internalization problem may be approached alternatively by launching simulations from preformed interior (cavity) states of the excess electron. Following the dynamics of such clusters, one may directly detect whether (a) these clusters persist for an appreciable amount of time or (b) transition from cavity states to surface states is observable in the simulations. We examine these possibilities in this section.

The preparation of the preformed cavity states follow a similar protocol that we described previously for the preparation of neutral ammonia clusters.[23] Here the same procedure is applied, but initially we place an extra chloride ion in the center of the simulation box surrounded by randomly placed and oriented ammonia molecules. To contract the cluster, we perform AM1-based molecular dynamics simulations at around 0 K with systematic quenching of the kinetic energy of the molecules on the approximate time scale of ∼1 ns. Once the (NH3)Cl– clusters are sufficiently contracted, they are gradually heated up to 100 K. At this temperature, we remove the chloride ion from these clusters, leaving a solvent void behind, and add a negative charge to the clusters. This type of molecular structure provides an energetically favorable binding site, a cavity, for the localization of the excess electron. Here we present the results of AIMD simulations that are started from these anions with preformed cavities. We note that we prepared clusters in the n = 8–48 range, but we were able to perform the actual simulations only up to n = 32.

The initial structures for these simulations are shown in Figure . Clearly, part of the excess electron spin density is located in the preformed cavity in each case. To characterize the spatial extent of the excess electron more quantitatively, we determined the volume surrounded by the 80% isovalue surface and calculated its radius assuming a perfectly spherical shape. The effective radii of the excess electron are estimated to be 5.6 Å (n = 8), 5.1 Å (n = 12), 4.5 Å (n = 16), 4.3 Å (n = 24), 3.0 Å (n = 32), and 2.6 Å (n = 48). For the smallest clusters (n = 8 and 12), significant portion of the spin density appears outside the cluster. As the clusters grow, the excess electron tends to shrink, to get gradually more localized. The two largest clusters (n = 32 and 48) very effectively confine the electron in the cavity.

Figure 6

Mesh representation of the spin density isosurfaces for typical n = 8–48 ammonia cluster anions with configurations that contain preformed cavities to accommodate the excess electron. (a) n = 8; (b) n = 12; (c) n = 16; (d) n = 24; (e) n = 32; (f) n = 48. The surfaces contain 80% of the total spin density with isovalues of 0.00004, 0.00005, 0.00006, 0.00007, 0.00012, and 0.00030, respectively. The electron distributions are computed with the 6-31(1+3+)G(d) basis set.

The structure of the cavity confined electrons also needs a short comment. Figure shows the spin density of the cavity state electron in the n = 48 cluster (see also Figure f) surrounded by the 15 nearest ammonia molecules forming the cavity. All the other molecules of the cluster have been removed from the figure for clarity. Here we again show the spin density isosurface that contains 80% of the total spin density. Previous studies showed that excess electron localization is negligible on the nitrogen atoms of dipole bound ammonia cluster anions.[23,24] The situation is clearly different here. Sizable excess spin density appears on the nitrogen atoms that surround the cavity electron. This picture appears to be qualitatively similar to that of Shkrob,[22] who proposed that the structure of the ammoniated electron may be better understood as a solvent stabilized multimer radical anion with a significant part of the excess electron density appearing in nitrogen frontier orbitals of the molecules that form the solvent cavity.

Figure 7

Mesh representation of the spin density isosurface for a typical n = 48 ammonia cluster anion that contains a preformed cavity to accommodate the excess electron (see also Figure f). The figure shows only the 15 nearest ammonia molecules that surround the cavity. The isosurface contains 80% of the total spin density with an isovalue of 0.00030. The electron distribution is computed with the 6-31(1+3+)G(d) basis set.

Molecular dynamics simulations from these starting configurations in the n = 8–32 range, however, reveal that the excess electron does not remain confined in these solvent voids. The molecular structure breaks, and concomitantly, the electron moves very quickly to the outside of the cluster. The temporal evolution of the VDE and the spin densities confirm this observation. Figure shows the progress of the VDE during the dynamics for two selected clusters, n = 16 and 32. At the outset the excess electron is significantly more strongly bound to the cluster than in surface bound states. The initial binding energies exhibit an increasing tendency with size, starting from 350 meV (n = 8) to 590 meV (n = 12), 820 meV (n = 16), and 1.32 eV (n = 24), and ultimately reaching 1.63 eV in the n = 32 cluster anion. Although we were not able to launch AIMD simulation from the n = 48 cavity containing cluster, we computed its initial VDE as 1.90 eV.

Figure 8

Time evolution of the VDE of (a) (NH3)16– and (b) (NH3)32– ammonia cluster anions that were prepared with a preformed cavity at 100 K. The MD run was generated using the BHandHLYP functional with the GPW-BS1 basis set. The VDE values of the cluster anion configurations along the trajectory are computed with the 6-31(1+3+)G(d) basis set.

The coupled nuclear and electronic dynamics, however, very quickly break the cavity, and push the electron outside the molecular frame. This process causes sharply decreasing VDE in time, melting down to about 0.1–0.2 eV by the end of the 1 ps long trajectories (Figure ). All investigated clusters with preformed cavities exhibit similar behavior. It may also be instructive to illustrate the stages of this relaxation by following the transformation of the spin density during the dynamics. Figure shows for the n = 32 cluster that the initial configuration (t = 0 fs) tightly binds the electron in the cavity. At the next illustrated instant (t = 125 fs), the electron still remains mainly in the cavity but distinctly grows in size. The VDE at this instant, although significantly diminished, is still high, ∼860 meV. By t = 250 fs, the VDE drops further to 280 meV, while the cavity opens up. By the end of the trajectory, the structure breaks to smaller structures, and the VDE remains at ∼0.3 eV. We note that for the n = 32 case we also repeated the simulation at a much lower temperature, 40 K, resulting in the same general trends and conclusions.

Figure 9

Mesh representation of the spin density isosurfaces of (NH3)32– generated in AIMD simulations. The simulation starts from a cluster structure with a preformed cavity (t = 0 fs). The figure also shows spin densities at t = 125, 250, and 950 fs. The surfaces contain 80% of the total spin density (isovalues of 0.00012, 0.00004, 0.000010, and 0.000008, respectively). The electron distributions are computed with the 6-31(1+3+)G(d) basis set.

Discussion and Conclusions

We performed ab initio molecular dynamics simulations on various ammonia cluster anions. The simulations focused on three distinct scenarios. First, we investigated the dynamics of small, dipole-bound anions of ammonia chains. The simulations suggest that these systems, particularly the longer (n ∼ 8) chains, may persist and be experimentally observable, especially at cryogenic temperatures. We note that in this size regime it is relatively straightforward to observe and identify couplings of the nuclear motions and the electronic degrees of freedom of the system. In the present case, we qualitatively assigned two main nuclear mode contributions to the VDE fluctuations, the intermolecular librational motions and the umbrella motions of the ammonia monomers, presumably at the electron binding end of the chain.

The second scenario involves addition of an extra electron to preformed low-temperature intermediate-size neutral nonlinear ammonia clusters. The general observation is that the smallest size such cluster anions (n = 8 and 12) appear to be unable to bind and stabilize an excess electron. These cluster anions easily dissociate or suffer electron detachment. Larger compact neutral clusters tend to only weakly bind the electron at the outset, although the strength of the interaction mildly increases with cluster size. The attachment of the excess electron to the neutral clusters takes place in very diffuse, delocalized orbitals. These states are concentrated mainly outside the molecular frame and engulf the molecular cluster. Nuclear dynamics may open the path to energy relaxation and localization of the excess electron. The first signs of such processes are observed for n ≥ 24 clusters. Within the limited length of our simulations we found the VDE to relax to ∼200 meV for n = 24 and ∼500 meV for n = 32 clusters. Nevertheless, from the present short simulations, it is impossible to anticipate what happens with the cluster anions at significantly longer time scales. One possibility is that the initially surface-bound diffuse, and delocalized excess electron stabilizes on the surface, then gradually gets internalized, moves to the interior of the cluster, and further stabilizes there. This possibility leads us to our third scenario of our studies.

Instead of attempting to directly observe excess electron internalization in ammonia clusters, we prepared ammonia clusters with preformed cavities and launched cluster anion simulations using these configurations. The results of these simulations are very clear. Although the electron is bound very strongly within these clusters initially, the cluster anions themselves appear to be unstable. Apparently, the electron quickly (within ∼200 fs) leaves the cavity and the interior of the cluster, and it becomes delocalized outside the molecular frame. At the same time, the cavity breaks, and the cluster opens up and may break into smaller fragments. The observation of this sequence of events strongly suggests that in the investigated size regime, surface stabilization is the likely mechanism of electron solvation.

As a corollary, and for comparison, Figure shows selected characteristic VDE values of the present study and those inferred from the linear relationship for the n–1/3 size dependence of the experimental VDEs found by Sarkas et al.[12] Although our examined cluster sizes (except for the only n = 48 cavity cluster) are smaller than the smallest size, experimentally observed ammonia cluster anions (n = 41), similar to our previous study, we can extrapolate the linear VDE−n−1/3 relationship of the experiment to the simulated sizes.[23] These extrapolated VDE values are estimated to be 290, 410, 490, and 590 meV for n = 16, 24, 32, and 48 clusters. Figure compares these values to (a) the strongest initial binding energies of the excess electron attached to preformed neutral clusters, (b) to the strongest VDE found during the relaxation of surface stabilized electrons during their dynamics, after they were attached to preformed neutral clusters, and (c) the initial VDE values of interior-state electrons that are confined in preformed cavities in the interior of ammonia clusters. The implications of the dynamical behavior and the VDE data can be summarized in the following points: (i) The preformed interior state clusters appear to be unstable. While the cavity dissolves, the clusters break to smaller clusters, and the electron spontaneously leaves the interior of the cluster. The initial binding energies of the excess electrons are way too high compared to those of the computed surface state data and the experimental data points. Nevertheless, we point out the close similarity of our interaction energies and those simulated by Barnett et al.[19] and Marchi et al.[17] for interior state clusters, and the sharp disagreement with the experimental photoelectric threshold energy of 1.3–1.4 eV for the bulk.[12] (ii) The smallest size compact neutral clusters are unable to bind an excess electron. (iii) Excess electrons can be bound to neutral ammonia clusters for sizes n ≥ 24 via diffuse and delocalized states that surround essentially the whole cluster. The initial attachment energies are relatively weak, even with respect to the experimentally observed data. (iv) Larger neutrals clusters that bind the excess electron initially may be able to stabilize the excess electron via coupled nuclear and electronic relaxation in more localized states that are concentrated mainly outside the molecular frame. The VDE values of such relaxed structures are comparable to the experimental data extrapolated for this size range and to early one-electron quantum-path-integral simulations for surface state bound species.[17,19]

Figure 10

VDE values for various size ammonia cluster anions. Red stars: experimentally determined VDE for the n = 41 cluster anion. Red squares: data points extrapolated from the experimentally determined linear VDE–n–1/3 relationship. Green circles: strongest initial binding energies of the excess electron attached to preformed neutral clusters at 100 K. Blue triangles: strongest VDE values found during the relaxation of surface stabilized electrons during their dynamics at 100 K, after they were attached to preformed neutral clusters. Open black triangles: initial VDE values of interior state electrons that are confined in a preformed cavity in the interior of the ammonia clusters at 100 K. VDE values are computed at the BHandHLYP/6-31(1+3+)G(d) level.

In summary, we conclude that the above observations are not consistent with excess electron binding in solvent cavities at around the investigated cluster size regime. The data also seem to suggest that the smallest experimentally observed clusters may be more prone to localize the excess electron predominantly outside the cluster, although part of the electron density of these states may also penetrate in the interior of the cluster. Partial and continuously increasing embedment of the excess electron in increasing size clusters may provide a reasonable explanation for the smoothly varying VDE tendency with cluster size. In fact, hydrated electron simulations indicated the possible existence of intermediately embedded binding cases that lay between the two limiting localization modes, interior states and surface states.[29] The relatively low experimental photoelectric threshold energy of the bulk ammoniated electron, however, still remains to be understood and explained. The satisfactory explanation of the phenomenon may require large-scale bulk-ammoniated electron AIMD simulations in the future.

Before concluding the paper, we pause and summarize the limitations of our simulation approach. It is clear that the present many-electron molecular dynamics simulations represent an improvement over the previous one-electron approaches.[17,19] Nevertheless, one has to be aware that numerous technical and theoretical approximations are still applied, for example, in the electronic structure method, the choice of the basis set and the various simulation machineries outlined in the Methods section. To be more specific, we list four aspects here: (1) We point out the choice of the electronic structure method (DFT) with a specific functional (BHandHLYP). (2) Within the DFT methodology, we neglected the dispersion correction between ammonia molecules. The ammonia–ammonia interaction energy is dominated mostly by hydrogen-bonding and dipole–dipole interactions. We estimated the dispersion to account for ∼10% of the intermolecular interaction. However, we found that the application of the D3 dispersion correction overestimates the interaction strength by ∼10% (see the Supporting Information). (3) On the basis of considerations on the computational resources, the size of the investigated systems, and the duration of the simulations, we generated our MD trajectories using a specific basis set (GPW-BS1). Although this basis set somewhat underestimates the strength of the electron binding, we demonstrated that the sampling of the configurations can be assumed to be similar to that of a somewhat more strongly binding (and more reliable) atom-centered basis set (6-31(1+3+)G(d)). Single-point VDE and spin density calculations along the trajectories were then performed using this latter basis set. (4) Beyond these, we also have to note the limited sampling of our simulations both in terms of the number of generated trajectories and the length of the trajectories.

Despite the approximations and technical shortcomings, we expect that the trends of the simulations remain valid and believe that the present simulations provide a solid, semiquantitative glimpse in the molecular details underlying electron solvation in ammonia clusters.