Flexible boundary layer using exchange for embedding theories. II. QM/MM dynamics of the hydrated electron

Flexible boundary layer using exchange for embedding theories. II. QM/MM dynamics of the hydrated electron Zhuofan Shen, 2 Shaoting Peng, and William J. Glover 2, 3, a) NYU Shanghai, 1555 Century Ave., Shanghai, 200122, China NYU-ECNU Center for Computational Chemistry at NYU Shanghai, 3663 Zhongshan Road North, Shanghai, 200062, China Department of Chemistry, New York University, New York, NY, 10003, USA


I. INTRODUCTION
The hydrated electron, e − aq , an excess electron embedded in liquid water, is the quintessential system exhibiting solvent-supported electronic states, and has been used as an ideal probe of water's solvation dynamics, since the electron itself has no internal degrees of freedom that might muddy a pump-probe spectrum. Nevertheless, a reasonable question to ask is whether the excess electron is simply a spectator in water's solvation dynamics, or whether it strongly perturbs the electronic and molecular structure of the surrounding water to an extent not thermally accessible in the neat liquid.
Much of our understanding of e − aq 's dynamics has come from one-electron mixed quantum classical (MQC) simulations, which consider the quantum mechanics only of the excess electron, while treating the solvent molecules classically. [1][2][3][4][5][6][7][8][9][10] The low computational cost of MQC methods has allowed for converged sampling of the equilibrium properties and dynamics of e − aq , as well as its nonadiabatic dynamics following photoexcitation. 2,3,5,8,10 Nevertheless, an underlying approximation in all current MQC descriptions of e − aq is that water's electronic structure and thus intramolecular potential energy surface is unchanged from the neat liquid. MQC is thus inherently incapable of addressing the question of whether the excess electron is a spectator, since it is assumed to be.
In addition to the assumption of a neat-liquid water electronic structure, MQC methods rely also on approximate electron-molecule pseudopotentials to describe the interaction of the excess electron with the closed-shell water molecules. While much effort has been devoted to their development, 1,4,6-9 these pseudopotentials introduce large sources of uncertainty in the accuracy of MQC descriptions of e − aq . Unlike typical atomic pseudopotentials that replace core electrons, the electron-molecule a) Electronic mail: william.glover@nyu.edu pseudopotentials replace valence electrons, which introduces an uncontrolled error. Furthermore, the structural predictions of MQC simulations of e − aq can depend quite sensitively on the pseudopotential, 8,9,11,12 which one of us found recently to arise from a competition of entropic and enthalpic driving forces in the solvation structure of e − aq . 13 A first-principles description of e − aq is thus preferred, and several ab initio molecular dynamics (AIMD) strategies have emerged recently, including Quantum Mechanics/Molecular Mechanics (QM/MM), 14,15 periodic planewave Density Functional Theory (DFT), 16,17 periodic second-order Möller-Plesset theory (MP2), 18 and Path Integral Molecular Dynamics (PIMD) with Machine-Learned (ML) potentials. 19 Given the computational cost of ab initio methods, explorations of the equilibrium dynamics of e − aq have so far been limited, although the initial localization dynamics has received attention. 14,16,18,20 An exception is the very recent PIMD study of e − aq with ML potentials, which afforded extensive sampling due to the forcefield-like cost of ML potentials. 19 Of the AIMD methods, QM/MM is particularly efficient, since most of the solvent molecules far from the electron can be treated with a cheap MM forcefield, and only the first few solvent shells around e − aq need to be treated at a QM level. However, a problem arises since water is diffusible, and without a special treatment, the QM and MM water molecules would mix, leading to a breakdown of QM/MM separation. This limits the length of trajectory that can be simulated to a few ps. The issue of QM/MM mixing by diffusible particles is precisely addressed by our Flexible Boundary Layer using Exchange (FlexiBLE) method, introduced in the previous companion paper. 21 FlexiBLE belongs to a class of constrained QM/MM methods that maintain QM/MM separation by applying a repulsive potential on particles near the QM/MM boundary. [21][22][23][24] We showed in paper I that with a careful choice of boundary potential, QM/MM separation can be maintained while leaving en-semble averages unaffected and allowing for density fluctuations in the QM region. Unlike adaptive QM/MM approaches, [25][26][27][28][29][30][31][32][33] FlexiBLE prevents QM or MM particles from traversing the boundary, meaning long-time diffusional dynamics will not be accurately described. However, we showed in paper I that FlexiBLE preserves equilibrium dynamics in the inner QM region, at least on sub-diffusional timescales, 21 which is sufficient to simulate e − (aq) 's vibrational spectroscopy. To explore the influence of an excess electron on water's electronic structure and dynamics, we performed many-electron FlexiBLE-QM/MM AIMD simulations at the BH&HLYP-D3/6-31++G* level on a model of the condensed-phase e − (aq) . We find that the equilibrium dynamics of the electron exhibits enhanced coupling to OH stretch modes compared to an MQC description. By performing Natural Bond Orbital (NBO) analysis, we show that this coupling arises from significant occupation of OH σ * antibonding orbitals, in a similar fashion to that observed in water anion clusters. 34 This allows us to develop a minimal frontier orbital description of e − (aq) involving a cavity orbital and strong coupling to 4-5 OH σ * orbitals. e − (aq) can thus be viewed as a strong hydrogen bond acceptor.
The remainder of the paper is as follows. In section II we benchmark DFT for electron-water interactions. In section III we describe the computational details of our condensed-phase model of e − (aq) . Results are presented in section IV, where we first apply FlexiBLE to a QM/MM treatment of neat water at the BH&HLYP-D3/6-31+G* level in section IV A and then to a QM/MM description of e − (aq) in section IV B. Dynamics of e − (aq) are explored in section IV C, and the NBO analysis is discussed in section IV D. Finally, conclusions are drawn in section V.

II. BENCHMARKING DFT DESCRIPTIONS OF ELECTRON-WATER INTERACTIONS
DFT electronic structure has a sufficiently low computational cost that merits its use in AIMD; however, for weakly-bound electronic states one might be concerned about the accuracy of standard functionals. In particular, delocalization error is known to be quite severe for cluster models of the hydrated electron. 35,36 This motivates a benchmarking of functionals in describing electron-water interactions.
One of us previously developed a method for computing single water-electron interaction energies with varying average relative displacement. 9 The idea is to include in the electronic Hamiltonian a confining potential of the form:V where k and k z control the strength of the confinement in the x, y and z directions respectively. In addition, a linear term inẑ is added with strength C z = −0.015167 Hartree Bohr −1 that serves to localize the excess electron on one side of the molecule. The +z direction was always chosen to point from the oxygen atom to the electron. The high order of the position operators in Eq. 1 introduces a steep repulsive wall surrounding the water molecule that confines the excess electron. This is needed since a single water molecule alone does not bind an excess electron. A value of k = 1 × 10 −7 Hartree Bohr −8 confines the excess electron within a fewÅ of the water, while leaving the neutral water's electronic structure essentially unperturbed. 37 k z is then varied between 1 × 10 −7 and 1 × 10 −10 Hartree Bohr −8 , which has the effect of moving the minimum in the confining potential, and therefore the excess electron's density, to increasing displacement from the water molecule (see dashed red curves in Fig. 2 below).
With the confining potential of Eq. 1, electronic structure calculations were then performed on a water anion. Since the excess electron localizes up to severalÅngstom from the water, it is important to augment the standard atom-centered basis functions to describe highly diffuse orbitals. Following our previous work, 9 we achieved this by adding to the aug-cc-pVTZ basis 38 a cubic grid of floating s type primitive Gaussian functions with exponent 0.16 Bohr −2 . The grid contains 7 × 7 × 13 functions with a spacing of 2.5 Bohr. The electron-water interaction energy is defined as: where E anion and E neut are the total energies of the anionic and neutral systems respectively, and E 1e is the ground-state eigenvalue of a single electron in the confining potential of Eq. 1. We found previously that unrestricted MP2 results agreed well with gold-standard Coupled Cluster Singles Doubles and perturbative Triples (CCSD(T)) 39 interaction energies, 9 so we use MP2 as the reference in this work. Density functional calculations were performed in Q-Chem 5.0.2, 40 with the SG-1 Exchange-Correlation (XC) integration grid. 41 Care had to be taken due to our use of a floating Gaussian basis formed with distributed ghost helium atoms that each contributed to XC integration grid centers. As discussed in detail in Supplementary Material Section S-I, SG-1's use of pruned angular grids in atomic core regions introduces large integration errors if a ghost atom is too close to a real atom. To avoid this issue, we removed the ghost atom at the origin, and added its s primitive function to the oxygen atom's basis. We consider first the BLYP functional, 42,43 since this was used in the first planewave DFT study of e − (aq) . 16 Fig. 1 plots the electron-water interaction energy for three different orientations (indicated by the molecular graphics at the top of the figure). Panel (a) shows the interaction energies for an electron on the oxygen side of water. As expected, the interaction energy is positive, corresponding to a repulsion from the partially negative oxygen atom. However, BLYP (dashed blue curves) sig-nificantly underestimates the magnitude of repulsion by ∼0.2 eV compared to the reference MP2 results (green circles) at high confinement strengths (small e − -water distances). Similar behavior is seen in two other orientations: dipole aligned (panel (b)) and OH bond aligned (panel (c)), which show attractive interactions as expected, but with BLYP overestimating the attraction by ∼0.2 eV when the electron is close to the water. The overestimated attraction between the electron and water at short range in BLYP is suggestive of the delocalization error. Indeed, Johnson and co-workers found significant convexity error in the energy versus fractional charge of a Kevan model of the hydrated electron with BLYP, raising the concern that delocalization error is severe when this functional is applied to solvated electrons. 35,36 We confirm this is the case in Fig. 2, which plots spin densities for the electron-water system. We indeed see that for all orientations and confining strengths, BLYP (dashed blue curves) underestimates the spin density of the excess electron in the attractive well region and overestimates the spin density on the water molecule compared to MP2 (green circles), with the error most severe at large confining strengths.
Johnson and co-workers found that using a hybrid functional with an increasing fraction of exact exchange reduced delocalization error. 35,36 Herbert and co-workers found similar behavior, with the SOMO energy of e − (aq) from MQC snapshots increasing with fraction of exact exchange from PBE to PBE0 to HF. 44,45 To test whether increasing the amount of exact exchange reduces delocalization errors and overbinding in our benchmark electronwater system, we plot interaction energies for the halfand-half hybrid, BH&HLYP, 46 as solid black curves in Fig. 1. We see much improved agreement to MP2 with this functional, with errors in the interaction energy not exceeding 0.06 eV, although a small amount of delocalization error remains, as seen in the solid black curves of Fig. 2, and noted also by Johnson and co-workers. 35 However, the good agreement between BH&HLYP and MP2 interaction energies in Fig. 1, in addition to Herbert and Head-Gordon's previous observation that this functional well reproduces vertical detachment energies of water cluster anions, 47,48 suggests that the remaining delocalization error is not too severe, and this functional is useful in describing e − (aq) .

III. COMPUTATIONAL DETAILS
We built a spherical nanodroplet model of the condensed-phase e − (aq) following a similar protocol to the companion paper. 21 A cavity for the excess electron was carved out by placing a chloride ion at the origin of a 25-Å radius droplet of 2017 water molecules created using the SolvateCap command of tleap in the Amber 2018 package. 49 The water forcefield was described by the SPC/Fw model, 50 and chloride Lennard-Jones parameters were taken from Ref. 51. A half-harmonic cap potential with a force constant of 10 kcal/mol/Å 2 prevented evaporation of water at the droplet/vacuum interface, and the chloride ion was restrained to the origin using a harmonic potential with a force constant of 1000 kcal/mol/Å 2 . Following minimization, heating, and equi-libration, a 100-ps NVT production run was performed to generate ten independent snapshots separated by 10 ps, from which the chloride ion was then removed to generate initial conditions for the e − (aq) simulations. To model e − (aq) , 32 water molecules nearest the origin of the droplet were treated at a QM level, with the remaining 1985 water molecules treated at an MM level. Since the switch from MM to QM introduced large forces on the QM water molecules, we quenched the largest components of their forces by performing 5 steps of minimization using the L-BFGS algorithm, 52 implemented in DL-FIND. 53 During dynamics, QM/MM separation was maintained with FlexiBLE implemented in a development version of Terachem, [54][55][56] and interfaced with OpenMM 7.3. 57 We used a FlexiBLE exponent parameter of α = 15Å −1 and a convergence parameter of γ = 10 −3 k B T . The QM region was described at the BH&HLYP-D3/6-31++G* level of theory, where D3 represents dispersion corrections with Becke-Johnson damping. 58 The QM and MM subsystems were coupled with the standard approach of elecrostatically embedding the QM particles in the field of the MM point charges. Mechanical embedding forces were also included via Lennard-Jones pair potentials between the QM and MM particles according to the SPC/Fw forcefield. 50 AIMD was propagated with the velocity Verlet algorithm 59 using the Niklasson time-reversible Born-Oppenheimer integrator with a 0.5-fs timestep. 60 Temperature was maintained at 298 K with the Bussi-Parrinello thermostat. 61 We generated ten trajectories of 7 ps in length. The first 2 ps of data of each trajectory was discarded as an equilibration period.
We note that one challenge arising in the simulation of e − (aq) is the occasional rapid diffusion of the electron to the QM/MM boundary. This issue was noted previously by Holden et al, 15 and we followed their solution of discarding trajectories that had the electron diffuse to the boundary, based on monitoring the outermost extent of the electron, defined as the sum of the electron's centroid distance: and its radius of gyration: where ψ SOMO is the singly-occupied molecular orbital, occupied by the excess electron. If the sum of these values exceeded 6.0Å, corresponding to the QM/MM boundary region, the trajectory was discarded. Fig. 3 shows an example of an accepted trajectory (black curve), in which the radial extent of the electron avoids the QM/MM interface. An example of a rejected trajectory (red curve) is also shown, in which the radial extent exceeded 6Å at ∼200 fs. Discarding trajectories this way would likely lead to an underestimate of the diffusion constant of e − (aq) and also limits the length of our simulations to no more than several ps; however, this is sufficiently long to study representative equilibrium dynamics of e − (aq) on sub-picosecond timescales, which governs its vibrational spectroscopy.

IV. RESULTS AND DISCUSSION
A. Liquid water structure Before considering e − (aq) , we first verified our simulation setup by performing AIMD on pure liquid water at the BH&HLYP-D3/6-31+G* level. To our knowledge, the structure of water at the BH&HLYP level has not been previously published. Using the same protocol discussed above, but without a chloride ion, we generated a spherical nanodroplet of 2028 waters. We fixed the oxygen atom of one QM water molecule at the origin, and allowed the other 31 QM and 1996 MM water molecules to move freely. FlexiBLE was used to maintain QM/MM separation, and 110-ps of NVT dynamics was propagated.
The Radial Distribution Functions (RDF) for liquid water relative to the central QM oxygen atom are shown in Fig. 4. Compared to the experimental data of Ref. 62, we see that the overall solvent structure of liquid water is well reproduced, albeit with some noticeable overstructuring in the first solvent shell compared to experiment. It is well known that such overstructuring arises from the neglect of Nuclear Quantum Effects (NQE), [63][64][65] and thus we view the RDFs in Fig. 4 as indicating that BH&HLYP-D3/6-31+G* provides a faithful representation of liquid water's potential energy surface. Furthermore, we see that FlexiBLE has correctly maintained QM/MM separation (indicated by solid red and blue filled areas), without introducing artefacts in the total RDF.

B. Hydrated electron structure
Having verified that our simulation protocol gives a good description of liquid water's structure, we turn next to the simulated properties of e − (aq) .  (aq) is thus essentially the same as found in previous QM/MM simulations at the self-interaction corrected BLYP-D level, 14 and consistent with the overall picture provided by cavity-forming MQC simulations, 4,7 and recent PIMD simulations with ML potentials, 19 although the latter's prediction of quasi-stable dual cavity structures is not observed here, perhaps due to the neglect of NQE.

C. Hydrated electron dynamics
We turn next to the equilibrium dynamics of e − (aq) , which we measure through the time autocorrelation function of its band gap fluctuations: where the angled brackets indicate a classical ensemble average and the fluctuations in the band gap are defined as: and E i is the energy of state i. For our DFT simulations, we take the band gap to be the SOMO-SOMO+1 gap. For the MQC simulations, the band gap is the excitation energy between the lowest two eigenvalues of the one-electron Hamiltonian using the Turi-Borgis electronwater pseudopotential, 4 as described in paper I. 21 Representative orbital plots of the SOMO and SOMO+1 are shown in Fig. 6(a) and (b) respectively, and are consistent with the s-like and p-like states of the particle-in-a-spherical-box model of e − (aq) originally put forward from MQC simulations, 1 and later confirmed by numerous many-electron calculations. 16,17,66,67 The density of states of the SOMO and SOMO+1 are shown in Fig. 6(c). The average BH&HLYP SOMO energy of -2.6 eV is compatible with the experimental binding energy of 3.7 eV, 68,69 given the limitations of Koopman's theorem. A slightly lower average SOMO energy of -3 eV was observed in HF QM/MM simulations; 15,45 however, we caution against a direct comparison of absolute SOMO energies between different studies, due to differences in solvation structures, QM region sizes, basis sets, MM models, and treatments of long-range interactions. From Fig. 6, the SOMO to SOMO+1 gap is seen to be ∼3 eV. This overestimates the experimental absorption peak of 1.7 eV, 70,71 likely as a result of the 50% fraction of exact exchange in the BH&HLYP functional. Timedependent density functional theory would provide a better description of electronic excitation energies of e − (aq) , 67 however at significant additional computational expense. For the purposes of revealing the influence of the excess electron on water's vibrational frequencies, our interest is in the dynamic fluctuations of the band gap, rather than the absolute magnitude of the gap, supporting our use of the orbital energy gap.
As discussed in paper I, C EE reports couplings between solvent motions and the energy levels of the excess electron, and FlexiBLE was seen to reproduce C EE from full-system calculations at the MQC level, indicating the bias potential did not measurably influence solvent motions that couple to the excess electron. 21 To directly reveal the electron-coupled solvent motions, we take the Fourier transform of C EE (t), and display the result in Fig. 7(a). We focus on the frequency range of 1000-4000 E (eV) cm −1 , which reports couplings to the intramolecular vibrations of water, and allows comparison to experimental Resonance Raman (RR) spectroscopy. 72 As Fig. 7(a) shows, the energy gap fluctuations of e − (aq) at the FlexiBLE QM/MM level (solid black curve) are dominated by two features: a broad peak centered at 3530 cm −1 , and a pair of peaks at 1500 and 1750 cm −1 . The 3530 cm −1 peak is assigned to water's OH stretch, while the 1750 cm −1 peak corresponds to water's bend. These numbers are somewhat higher than the experimental e − (aq) RR peak frequencies of 1610 and ∼3200 cm −1 ; 72 however, applying the recommended vibrational frequency scaling factor of 0.9374 for BH&HLYP 73 brings the theoretical frequencies into better agreement with experiment. We will address the 1500 cm −1 peak further below.
Comparing QM/MM to MQC results (solid black and dashed red curves of Fig. 7(a)  The spectral density of pure water with the SPC/Fw forcefield, shown as the dot dashed blue curve in Fig. 7(b), reveals the origin of the peak at 1500 cm −1 seen in the gap correlation function of e − (aq) at the Flexi-BLE QM/MM level (solid black curve in Fig. 7(a)): this peak arises from coupling between the excess electron and the MM SPC/Fw water molecules' bend. SPC/Fw was parameterized to reproduce structural, thermodynamic, and kinetic, rather than vibrational, properties of liquid water, 75 and the bend frequency is evidently underesti-mated in this model. Although our use of the SPC/Fw model results in a spurious secondary water bend peak in e − (aq) 's energy gap fluctuations, we view this as serendipitous, as it allows us to easily separate the coupling of the excess electron to QM and MM water bend vibrations, which we see to be dominated by the former.
In addition to the redshifts between QM/MM and MQC results, another difference is that the total amplitude of coupling to the OH stretch at the QM/MM level is significantly increased compared to MQC. Integrating the modulus squared of the Fourier transform between 2500 and 4000 cm −1 reveals that the coupling to the OH stretch is 8.6 times stronger in the QM/MM model compared to MQC. On the other hand, the coupling to the bend is comparable: between 1616 and 2500 cm −1 , the MQC integral is 1.9 times larger than QM/MM. The ratio of MQC:QM/MM bend intensity reduces to 1.1 if the second QM/MM peak is included by extending the lower integration limit to 1350 cm −1 .
The observation of electron induced vibrational redshifts and a significant enhancement of coupling to water's OH stretch compared to MQC results suggests that the excess electron interacts with water much more intimately than simple electrostatics would suggest. The vibrational redshifts arise from occupation of σ * OH orbitals, 44 the energy of which is strongly modulated by the OH stretch mode. It is thus reasonable to expect that σ * OH occupation also explains the enhanced coupling between the electron and stretch modes. To explore this further, in the next section we seek to quantify the occupation of σ * OH orbitals, using NBO analysis.

D. NBO analysis
We performed NBO analysis at the BH&HLYP/6-311++G** level on e − (aq) snapshots separated by 100 fs extracted from our BH&HLYP-D3 QM/MM trajectories using NBO 7.0.8, 76 interfaced with Q-Chem 5.0.2. 40 Special care is needed in the NBO analysis of solventsupported states, such as e − (aq) , since the default approach is to optimize spin NBOs to each have approximately unit occupation. For e − (aq) , this resulted in a single water close to the electron having a Lewis structure of H OH -; i.e., with a broken OH bond, and with the excess charge assigned to the resulting hydroxide ion. Since we do not observe OH bond breaking in our AIMD simulations, we view this Lewis structure as unphysical. This behavior was avoided by specifying in the input a Lewis structure that leaves water molecules intact.
A second issue in the NBO analysis of e − (aq) relates to the diffuse nature of the excess electron's orbital and its occupation of a solvent cavity. Using an exclusively atom-centered basis augmented with diffuse functions (6-311++G**), we found significant mixing in of Rydberg-like character in the σ * OH NBOs of the water molecules closest to the electron (see Supplementary Ma-terial Fig. S2), resulting in artificially high σ * OH occupations of >0.3 e. We found this issue was circumvented by adding a ghost atom with lithium's 6-311++G** basis functions (chosen as the 2nd-row element with the smallest diffuse exponent) at the centroid position of the SOMO. We then assigned the excess electron to an NBO formed largely from the lithium ghost atom's outer valence and diffuse set, with smaller contributions from Rydberg natural atomic orbitals of the water atoms. We identify this ghost atom NBO as an s-like cavity orbital, Cav s . Representative σ * OH and Cav s NBOs are shown in Fig. 8. Convergence of NBO properties with basis set and QM region size is explored in Section S-II, where we show that our chosen 6-311++G** basis provides comparable results to aug-cc-pVQZ, but at significantly lower computational cost, allowing the full QM region of 32 water molecules to be included.
We start by considering average NBO populations of the Cav s and σ * OH orbitals after ordering α-spin NBOs by population for each snapshot, with the highest populated antibonding orbital labelled as σ * 1 . The result is shown in Table I. As expected, the Cav s orbital has the highest population; however, σ * OH orbitals contribute an appreciable ∼0.3 e in total. In fact, the Cav s and five σ * OH orbitals sum to a population of 0.97 e, and thus together describe the excess electron's SOMO orbital almost in its entirety. The remaining amplitude of the SOMO is made up of contributions from lower-populated σ * OH orbitals and water Rydberg orbitals. This analysis is in qualitative agreement with the picture of excess electron binding in cavity-containing water cluster anions. In particular, Herbert and Jacobson noted that only ∼0.5 e of excess charge density is localized in the cavity of the idealized (H 2 O) − 24 cluster, with the remainder penetrating into the second solvation shell and beyond. 77 Using a QM/MM model, Uhlig, Marsalek, and Jungwirth (UMJ) confirmed this picture extends to the condensed-phase e − (aq) and assigned 0.41 e of the excess charge to the cavity region, 0.24 e to a region overlapping with the first two solvation shells water molecules, and 0.35 e to a diffuse tail extending beyond the second solvation shell. 14 We assign a larger contribution to a cavity state than UMJ, since the Cav s NBO also contains the diffuse tails that UMJ assigned separately. Our analysis provides additional quantitative insight into how the excess electron is able to have appreciable amplitude in the first two solvation shells, which is by populating the σ * OH orbitals.
To further explore the magnitudes of σ * OH occupations in e − (aq) , we found it revealing to project the distribution of occupations onto the e − -oxygen distance and e − -O-H angle, where the excess electron's position is taken to be the centroid of the SOMO. This is shown as the false color map in Fig. 9. Rather surprisingly, we see σ * OH occupations on a single OH bond reach 0.2 e. This high occupation occurs exclusively when a water molecule is close to the electron (e − -O distance < 2Å) and with the OH bond aligned to the electron (e − -O-H angle < 20°). by the σ * OH orbitals having greatest amplitude on the hydrogen atom side, which then due to phase matching, constructively interfere with the cavity orbital when they point towards the cavity (see Fig. 8). This is further confirmed by plotting transition densities between the Cav s and the closest σ * OH orbitals, which are seen to point from the hydrogen atoms towards the center of the cavity (see Fig. S4).
Although a relatively rare event (see Fig. S3), σ * OH occupations of 0.2 e are quite surprisingly high, and exceed values found even in H 2 O-OHcomplexes (0.17) 78 and water cluster anion isomers exhibiting the double acceptor (AA) motif (0.16 e). 34 The latter are associated with water stretch vibrational redshifts of >300 cm −1 . 79 It stands to reason then that similarly large, or larger, vibrational redshifts occur in the condensedphase e − (aq) . The overall 200 cm −1 redshift we observe in the OH stretch of e − (aq) thus appears to result from the smaller average σ * OH populations of ∼0.05-0.1 e seen in the first-shell water molecules. However, the spectral envelope of the OH stretch is seen to extend down to below 3000 cm −1 (see Fig. 7(a)) and we hypothesize that these strongly red-shifted frequencies are due to the water molecule closest to the electron when it has a high σ * OH occupation. This idea will be explored in more detail in a future publication.
The observed large σ * OH occupations of 0.20 e might seem surprising considering the energy levels involved. The hydrated electron state has an energy of −3.7 eV relative to vacuum, 68,69 while the σ * OH orbitals are several eV above the vacuum level. 80 For the hydrated electron state to have significant σ * OH occupation then suggests a very large Hamiltonian coupling (several eV) of the excess electron to the σ * OH orbitals. However, significant populations of σ * OH orbitals also observed in halide ion water clusters, 48 which have donor orbitals even lower in energy than e − (aq) , suggesting that large Hamiltonian couplings between anions and water σ * OH orbitals are not uncommon. NBO analysis allows for a determination of Hamiltonian couplings by computing matrix elements of the Fock operator in the NBO basis. We averaged the matrix elements over our FlexiBLE QM/MM trajectories. 81 The results are presented in Table II which shows matrix elements between the Cav s orbital and the five most populated σ * OH orbitals. For simplicity, matrix elements coupling the σ * OH orbitals with each other are omitted. Since these values are found to be sensitive to the choice of density functional (see Supplementary Material section S-III), we caution against a quantitative physical interpretation of the absolute values of the model Hamiltonian; however, relative trends are preserved across functionals, resulting in σ * OH populations that do not depend sensitively on the functional (see Fig. S7). We consider first the diagonal elements of the model Hamiltonian in Table II, where we observe that all orbitals are above the vacuum level. This makes physical sense for the σ * OH orbitals, but might at first seem surprising for the cavity orbital. However, this cavity orbital is by construction orthogonal to all water-centered occupied and virtual NBOs and is thus rather confined in the cavity region and interstitial regions between the water molecules, raising its kinetic energy.
Considering the off-diagonal elements of the Hamiltonian in Table II, we indeed find that the Cav s and σ * OH orbitals are coupled by up to several eV, with the largest average coupling of 4.49 eV found for σ * 1 , which is also the closest OH bond to the electron. Similarly large couplings are seen for σ * 2 through σ * 4 , but by σ * 5 , the coupling has dropped to 0.33 eV. These values are consistent with the inner solvent shell of e − (aq) comprising four OH-aligned water molecules. The couplings may be viewed rather like Charge-Transfer (CT) integrals between the donor Cav s orbital and acceptor σ * OH orbitals. In this CT context, the large magnitude of the couplings can be understood as resulting from two aspects of e − (aq) . First, the diffuse nature of the Cav s orbital provides a high degree of spatial overlap with the acceptor σ * OH of the first solvent shell, which explains why the coupling drops off for the fifth and further σ * OH orbitals. Second, the transition density connecting the Cav s and σ * OH orbitals is localized near the partially positive H atoms (see Fig. S4), which results in an attractive Coulombic coupling.
Diagonalizing the Hamiltonian in Table II yields a lowest eigenvalue of -1.08 eV, which while higher than the true average SOMO energy of -2.6 eV, shows that our simple model captures most of the stabilization of the cavity orbital by strongly coupling to four σ * OH orbitals (see Fig. S5 for convergence of the eigenvalue with number of included σ * OH orbitals and Fig. S6 for its correlation with the SOMO energy). 82 The energy levels of the NBOs and model Hamiltonian eigenvalues are summarized in the MO diagram presented in Fig. 10. We notice an interesting stabilization of the most populated σ * 1 orbital's energy of ∼0.9 eV relative to the other σ * OH orbitals, that all have average energies of ∼14.2 eV. We found this stabilization arises due to a lengthening of the closest OH bond to the electron, which had an average length of 0.999±0.002Å compared to the remaining QM OH bonds which had an average length of 0.9761±0.0003Å. σ * 2 is stabilized by ∼0.3 eV in a similar manner, but to a lesser degree. Similar OH bond-length extensions were observed in the zero-Kelvin minimal cluster with polarizable continuum model put forward by Sevilla and co-workers. 83 Thus, the partial occupation of water σ * OH orbitals is seen to noticeably weaken the coordinating OH bonds, which lowers the σ * OH orbital energies and further enhances their coupling to the excess electron. This physics is missed entirely by MQC models, and provides an explanation for the significant differences in e − (aq) 's dynamics we observe at the MQC and QM/MM levels.

V. CONCLUSIONS
In this paper, we applied the FlexiBLE embedding approach, developed in the preceding companion paper, 21 to perform QM/MM dynamics of the hydrated electron at the BH&HLYP-D3/6-31++G* level. We found agreement with the structural picture of e − (aq) put forward by UMJ, 14 that the excess electron occupies a cavity, but has a complex character involving significant over-  Table II. lap with first-shell coordinating water molecules. Using NBO analysis, we provided a quantitative rationale for this complex behavior: the excess electron has ∼0.3 e occupation of σ * OH orbitals in total brought about by large electronic couplings of a cavity state and the antibonding orbitals. This leads to an enhanced dynamical coupling between the excess electron's energy levels and water stretch vibrations compared to a MQC model. Furthermore, the weakening of OH bonds by σ * OH occupation leads to stretch and bend vibrations that are redshifted compared to pure water, in quantitative agreement with experiment. 72 Somewhat surprisingly, the σ * OH occupation of a single water molecule can reach 0.2 e, which leads to a noticeable extension of that OH bond and a reduction of its antibonding orbital energy by ∼1 eV.
To answer the question posed in the introduction, we conclude that e − (aq) is not a simple spectator species, but rather is intimately coupled to, and strongly modulates, the motions of the first solvation shell waters. Previous MQC dynamics studies, which inherently miss this physics, emphasized the role of water translations and librations in the solvation dynamics of e − (aq) . Our results motivate a deeper consideration of intramolecular water motions in the dynamics of e − (aq) . Overall, we believe that this successful first application demonstrates that Flexi-BLE opens the door to studying the dynamics of many other condensed-phase species, for which a QM description of the environment is necessary.

SUPPLEMENTARY MATERIAL
See the supplementary material for an explanation of DFT grid errors when using floating Gaussian basis functions, and further details on our NBO analysis and model Hamiltonian.