Widely used continuum solvation models for electronic structure calculations, including popular polarizable continuum models (PCMs), usually assume that the continuum environment is isotropic and characterized by a scalar dielectric constant, *ε*. This assumption is invalid at a liquid/vapor interface or any other anisotropic solvation environment. To address such scenarios, we introduce a more general formalism based on solution of Poisson’s equation for a spatially varying dielectric function, *ε*(**r**). Inspired by nonequilibrium versions of PCMs, we develop a similar formalism within the context of Poisson’s equation that includes the out-of-equilibrium dielectric response that accompanies a sudden change in the electron density of the solute, such as that which occurs in a vertical ionization process. A multigrid solver for Poisson’s equation is developed to accommodate the large spatial grids necessary to discretize the three-dimensional electron density. We apply this methodology to compute vertical ionization energies (VIEs) of various solutes at the air/water interface and compare them to VIEs computed in bulk water, finding only very small differences between the two environments. VIEs computed using approximately two solvation shells of explicit water molecules are in excellent agreement with experiment for F^{−}(aq), Cl^{−}(aq), neat liquid water, and the hydrated electron, although errors for Li^{+}(aq) and Na^{+}(aq) are somewhat larger. Nonequilibrium corrections modify VIEs by up to 1.2 eV, relative to models based only on the static dielectric constant, and are therefore essential to obtain agreement with experiment. Given that the experiments (liquid microjet photoelectron spectroscopy) may be more sensitive to solutes situated at the air/water interface as compared to those in bulk water, our calculations provide some confidence that these experiments can indeed be interpreted as measurements of VIEs in bulk water.

## I. INTRODUCTION

Fundamental aspects of ion solvation at the air/water interface have attracted significant attention in recent years,^{1–7} including investigations of how ion coordination motifs, concentrations, and reactivity may differ at the interface versus bulk water. At the same time, the development of liquid microjet photoelectron spectroscopy has opened the way to experimental measurements of vertical ionization energies (VIEs) of molecules in solution,^{8–11,17} as opposed to the gas-phase VIEs afforded by traditional photoelectron spectroscopy. However, interpretation of solution-phase photoelectron spectra is complicated by the possibility that the ejected electron may be scattered and/or recaptured by the liquid and thus detected with reduced kinetic energy or possibly not detected at all. As such, the microjet experiments are likely more sensitive to species solvated at the liquid/vapor interface than they are to the same species in a bulk liquid environment. The wavelengthdependent nature of the electron attenuation length^{12,13} (a measure of the likelihood that the emitted electron is recaptured) leads to solution-phase photoelectron spectra that depend on the wavelength of the photodetachment laser.^{14,15} For these reasons and others,^{16} theoretical prediction of VIEs in solution is desirable in order to facilitate the interpretation of the experiments. From a quantum chemistry point of view, one may expect a significant polarization response from the medium upon ionization of the solute, so the question arises how this effect can be incorporated in a tractable way. A continuum representation of the solvent represents one cost-effective strategy.

The most common continuum solvation models in quantum chemistry are based upon the framework of the polarizable continuum model (PCM),^{19–23} in which the continuum solvent’s electrostatic interaction with the atomistic solute is parameterized in terms of a single, scalar dielectric constant, *ε*. For accurate solvation energies, nonelectrostatic interactions must be included as well, but the electrostatic contribution can still be obtained from a PCM.^{23–27} These models are simple, efficient, and—assuming nonelectrostatic corrections are included—reasonably accurate,^{20,23–27} and are therefore widely used in quantum chemistry. As conventionally formulated, however, these models assume that the solvation environment is isotropic, as appropriate for solvation in bulk liquid but not at an interface.

There have been some attempts to modify the PCM formalism for use in anisotropic environments, including a formulation that uses a dielectric tensor in place of a scalar dielectric constant.^{28–31} This is useful, e.g., in the case of liquid crystals where the dielectric “constant” is strongly direction-dependent and therefore a diagonal tensor with different values *ε*_{xx}, *ε*_{yy}, and *ε*_{zz} might afford a reasonable description. Mennucci *et al.*^{32–34} and others^{35,36} have developed PCMs designed for liquid/vapor interfaces by modifying certain matrix elements of the PCM equations within the interfacial region. Mennucci *et al.* used a smooth switching function to interpolate between liquid and vapor values of *ε*,^{32–34} as is also done in the approach presented here. For complete generality, however, an anisotropic continuum environment should be described theoretically using Poisson’s equation,^{22} not with a scalar (or tensor) dielectric *constant* but rather with a spatially varying dielectric *function*, *ε*(**r**).

A simplified version of such a model, in which *ε*(**r**) is replaced by a set of distinct dielectric constants *ε*_{1},*ε*_{2}, … in different spatial regions, was introduced long ago by Sakurai *et al.*^{37,38} and used in semi-empirical electronic structure calculations.^{38} At its core, this model amounts to solution of Poisson’s equation in each spatial region, subject to appropriate boundary conditions. In the present work, we introduce an even more general formalism and computational algorithm in which the function *ε*(**r**) is allowed to be completely arbitrary. It is ultimately defined by the value of *ε* at each point on a discretization grid.

Other solvers for Poisson’s equation have been reported recently,^{39–43} including several for use with quantum chemistry.^{40–43} What is novel in the present work is the introduction of nonequilibrium corrections. These account for the response of the continuum solvent to a sudden change in the electron density of the solute, such as that which occurs upon (vertical) ionization.^{44–47} We have previously formulated this nonequilibrium theory for use with PCMs,^{48–50} and here we make the appropriate modifications for use with Poisson’s equation. A preliminary version of this methodology was reported in Ref. 42, but whereas that formulation was perturbative (following along the lines of our group’s previous work on PCMs^{48}), the present version includes the full solvent response. We have also made significant improvements to our grid-based Poisson solver, as compared to the one described in Ref. 42.

This work aims to evaluate the limitations of nonequilibrium anisotropic Poisson boundary conditions in quantum chemistry calculations, by comparing to aqueous-phase VIEs measured using liquid microjet photoelectron spectroscopy.^{11} Perhaps unsurprisingly, VIEs for atomic ions computed using nothing but a PCM representation of the solvent afford extremely poor agreement with experiment;^{18} hence, we will include explicit water molecules in the atomistic, quantum-mechanical (QM) region. In PCM calculations, where the solute cavity that defines the solute/continuum interface is usually constructed from atom-centered spheres,^{19–23} inclusion of a large number of explicit solvent molecules sometimes leads to erratic convergence with respect to the size of the QM region.^{51} This occurs because the dielectric medium artificially intrudes into the interstices between explicit solvent molecules, which should properly be characterized by *ε* = 1 since these are part of the QM region. We avoid such artifacts by using a single, spherical solute cavity around the entire QM region or alternatively using a novel cavity construction that is described herein. Finally, we consider whether VIEs computed in bulk water differ appreciably from those obtained at the air/water interface. This is easily addressed computationally but less trivial to interrogate experimentally, although experimental insight might be gained from angle-resolved photoelectron spectroscopy.^{52–55}

## II. NONEQUILIBRIUM POISSON FORMALISM

### A. Self-consistent equilibrium solvation

Solution of the gas-phase Poisson equation,

affords the electrostatic potential *φ*_{sol}(**r**) arising from the electronic and nuclear components of the solute’s charge density, *ρ*_{sol}(**r**). The quantity *ρ*_{sol}(**r**) is to be computed from a quantum chemistry calculation, but we make no assumptions about the level of electronic structure theory. The quantities *φ*_{sol} and *ρ*_{sol} can be partitioned into electronic and nuclear components,

The solute’s internal energy in a vacuum is

and includes the electron–electron, electron–nuclear, and nuclear–nuclear interactions. Upon immersion of the solute in continuum solvent, characterized by a spatially varying dielectric function *ε*(**r**), the solute–solvent interaction is governed by the most general form of Poisson’s equation,

where

includes an induced polarization potential, *φ*_{pol}(**r**).

For electronic structure methods using atom-centered Gaussian basis functions *g*_{μ}(**r**), the electronic contribution to the electrostatic potential is

where **P** is the one-electron density matrix and *φ*_{μν}(**r**) is the electrostatic potential generated by the shell pair *g*_{μ}(**r**) *g*_{ν}(**r**),

Nuclear charges are smeared out using Gaussian functions to avoid problems with discretizing them onto a grid. The electrostatic potential generated by these Gaussian nuclear charges is

where *Z*_{α} and **R**_{α} are the charge and position of nucleus *α* and *σ* is the standard deviation of the Gaussian, which is an input parameter to the method.

The solute’s charge density is obtained from *φ*_{sol}(**r**) according to

In this work, $\u2207^2\phi sol$(**r**) is computed using an eighth-order central finite-difference scheme, as described in Sec. III A. Unlike our original implementation of Poisson boundary conditions,^{42} which required the direct evaluation of the electron density on the real-space grid, the present implementation evaluates *φ*_{elec}(**r**) on the grid, via Eq. (2.6), and then computes *ρ*_{sol}(**r**) from Eq. (2.9). We find that the present approach is more robust with respect to changes in the grid size and spacing.

To solve the Poisson problem in Eq. (2.4), we adapt a procedure outlined in Refs. 40 and 41 for obtaining the solvent polarization response, which is characterized by the quantities *φ*_{pol}(**r**) and *ρ*_{pol}(**r**). Equation (2.4) is first recast as a vacuum-like Poisson equation,

where the total charge density is

Note carefully the difference between Eq. (2.10) and Eq. (2.1). The effects of the inhomogeneous dielectric function *ε*(**r**) are contained in the polarization charge density *ρ*_{pol}(**r**), the form of which is^{40,41}

The first term on the right side of Eq. (2.12) is the solute charge density scaled by a dielectric-dependent factor that is only non-zero outside of the atomistic region (solute cavity), where *ε* > 1. The second term *ρ*_{iter}(**r**) is a charge density induced by the inhomogeneous dielectric in regions where it transitions from *ε* = 1 near the solute molecule to a value appropriate for bulk solvent outside of the solute cavity. As the notation implies, this correction is obtained iteratively, and its value at the *k*th iteration can be expressed as^{40,41}

Algorithm 1 outlines a procedure for the iterative solution of Eq. (2.10) to obtain *φ*_{tot}(**r**), *φ*_{pol}(**r**), and *ρ*_{pol}(**r**). We call this the “equilibrium” Poisson-equation solver (PEqS) method, which we now describe. The quantities *φ*_{sol}(**r**) and *ρ*_{sol}(**r**) are initialized using Eqs. (2.2a) and (2.6)–(2.9), then $\rho iter(k)$(**r**) is computed using Eq. (2.13), and Eq. (2.12) is then used to generate $\rho pol(k)$(**r**). With the total density now in hand, the total electrostatic potential is obtained via the numerical solution of Eq. (2.10) using a multigrid conjugate gradient (CG) procedure that is described in Sec. III B. This affords $\phi tot(k+1)$(**r**), and the iterative part of the charge density is then updated using Eq. (2.13). However, rather than using this directly to define $\rho iter(k+1)$(**r**), we instead use a damping procedure to stabilize the update between iterations *k* and *k* + 1,

We take *η* = 0.6 as in Refs. 40 and 41. Convergence of the solvent polarization response is achieved when the residual

falls below a threshold, *τ*_{PEqS}.

1: Initialize Δh = 0. |

2: for j = 1, 2, … do until SCF error < τ_{SCF} |

3: Diagonalize $F(j)=F0(j)+\Delta h(j)$ to obtain P^{(j)} |

4: Compute ρ_{sol}(r) and φ_{sol}(r) |

5: if j = 1 then |

6: ρ_{tot}(r) = ρ_{sol}(r) |

7: φ_{tot}(r) = φ_{sol}(r) |

8: else |

9: ρ_{tot}(r) = ρ_{sol}(r) + ρ_{pol}(r) |

10: φ_{tot}(r) = φ_{sol}(r) + φ_{pol}(r) |

11: end if |

12: Initialize $\rho iter(0)$(r) using $\phi tot(0)$(r) |

13: for k = 0, 1, … do until Δρ_{iter} < τ_{PEqS} |

14: Compute $\phi tot(k+1)$(r) |

15: Update $\rho iter(k+1)$(r), ρ_{pol}(r), and ρ_{tot}(r) |

16: end for |

17: Update Δh^{(j)} |

18: Compute E = E_{int} + G_{elst} |

19: end for |

1: Initialize Δh = 0. |

2: for j = 1, 2, … do until SCF error < τ_{SCF} |

3: Diagonalize $F(j)=F0(j)+\Delta h(j)$ to obtain P^{(j)} |

4: Compute ρ_{sol}(r) and φ_{sol}(r) |

5: if j = 1 then |

6: ρ_{tot}(r) = ρ_{sol}(r) |

7: φ_{tot}(r) = φ_{sol}(r) |

8: else |

9: ρ_{tot}(r) = ρ_{sol}(r) + ρ_{pol}(r) |

10: φ_{tot}(r) = φ_{sol}(r) + φ_{pol}(r) |

11: end if |

12: Initialize $\rho iter(0)$(r) using $\phi tot(0)$(r) |

13: for k = 0, 1, … do until Δρ_{iter} < τ_{PEqS} |

14: Compute $\phi tot(k+1)$(r) |

15: Update $\rho iter(k+1)$(r), ρ_{pol}(r), and ρ_{tot}(r) |

16: end for |

17: Update Δh^{(j)} |

18: Compute E = E_{int} + G_{elst} |

19: end for |

Operationally, the solvent polarization response is included in the QM calculations by augmenting the gas-phase Fock matrix **F**_{0} with a correction Δ**h** to its one-electron part. This correction has matrix elements

Finally, the converged solution of Eq. (2.10) affords a total energy

which consists of the solute’s internal energy *E*_{int} from the electronic structure calculation, plus the electrostatic contribution

to the solvation free energy.

### B. State-specific nonequilibrium solvation

To incorporate solvent polarization effects following vertical ionization of the solute, we have adapted the nonequilibrium solvation approach developed for PCMs,^{45,48–50,56–58} for use with three-dimensional charge densities rather than the apparent surface charges used by PCMs. In previous work, we developed a perturbative approach to correcting the solute/continuum interaction for a sudden change in the electronic state of the solute, either electronic excitation or ionization.^{48–50} The perturbative approach is tantamount to “freezing” the inertial components of a reference-state solvent reaction field such that, upon vertical ionization, this frozen reaction field governs the solvent response to the ionized state’s charge distribution, which is prevented from relaxing. In the case of electronic excitation, the perturbative nature of the correction avoids a state-switching problem that arises in the case of near-degeneracies when using a state-specific Hamiltonian,^{59} and in Ref. 42 we introduced a Poisson equation solver based on the perturbative approach.

For VIEs, the state-switching problem is not an issue and in this work we develop a Poisson equation solver based on a state-specific, nonequilibrium treatment of solvent polarization, rather than a perturbative approach. Within the state-specific approach, the final (ionized) state’s charge distribution relaxes in the presence of the slow inertial component of the reference-state reaction field as well as the fast noninertial component of its own reaction field.

For the state-specific method, the solute wave function for state |Ψ_{i}⟩ is obtained by solving the Schrödinger equation

with a state-specific Hamiltonian of the form

Subscripts indicate whether a particular quantity originates from the equilibrium reference state (*i* = 0) or else the nonequilibrium final state. (For ionization to the electronic ground state, there is only one possible final state that we will indicate by *i* = 1 in what follows.) Superscripts “slow” versus “fast” in Eq. (2.20) indicate which part of the solvent response is considered: either the slow inertial part, representing nuclear degrees of freedom (orientational and vibrational fluctuations of the solvent), or else the fast electronic part. The quantity $\u0124ivac$ is the molecular Hamiltonian that affords the solute’s vacuum internal energy *E*_{int,i} for state *i*, and the operators $V^islow$ and $V^ifast$ generate the indicated components of the solvent polarization response, i.e., $\phi pol,islow$ and $\phi pol,ifast$, where

As in previous work,^{48–50} we use the Marcus partition of the fast and slow components of the polarization response. (See Ref. 49 for a comparison to the common alternative, Pekar partitioning, with the conclusion that this choice makes essentially no difference for solvation energies.) Within this approach, the slow component of the reference-state polarization charge density, $\rho pol,0slow$(**r**), which affords $\phi pol,0slow$(**r**) according to Eq. (2.21), is computed according to^{48,49}

where *χ* = *χ*_{slow} + *χ*_{fast} is the static susceptibility, partitioned into slow and fast components,

Here, *ε*_{solv} is the static dielectric constant of the solvent and *ε*_{∞} = *n*^{2} is its optical dielectric constant, where *n* denotes the solvent’s index of refraction. The quantity *ε*_{∞} encodes the fast electronic contribution to the solvent polarization response.

To obtain the fast components of the ionized state’s polarization response, $\phi pol,1fast$(**r**) and $\rho pol,1fast$(**r**), within the Marcus partition,^{49} the Poisson equation is modified such that the total source-charge density is the ionized solute’s charge density, $\rho sol,1(r)+\rho pol,0slow$(**r**), and the dielectric function is the optical one. This modified form of Eq. (2.4) is

Here, $\phi tot,1fast$(**r**) is the *total* fast component of the ionized state’s electrostatic potential. To apply the equilibrium PEqS procedure introduced in Sec. II A, Eq. (2.24) is rewritten as

where $\rho tot,1fast$(**r**) and $\phi tot,1fast$(**r**) are to be computed self-consistently. The total nonequilibrium source charge density is

For the Marcus partitioning scheme,^{49} $\rho pol,1fast$(**r**) takes the form

where *ρ*_{iter,1}(**r**) is computed iteratively according to

Between iterations, we apply the damping procedure of Eq. (2.14).

Finally, the nonequilibrium free energy is^{48,49}

The term

arises within the Marcus partitioning scheme^{49} due to Coulomb interactions between fast and slow components of the solvent polarization response,^{19,46,48,49} in which the slow components of the reference-state response affect the fast components of the final-state response.^{45,60} The quantity *G*_{elst,1} defined in Eq. (2.29) is added to the gas-phase internal energy *E*_{int,1} to generate the electrostatic interaction energy of the final (ionized) state, *E*_{elst,1}. The state-specific nonequilibrium VIE is then evaluated as the difference between the ionized- and reference-state electrostatic energies, *G*_{elst,1} − *G*_{elst,0}.^{48,49} The result is

where

is the difference between the ionized- and reference-state internal energies.

Operationally, the gas-phase Fock matrix for the ionized state must be corrected for the solvent response using a matrix Δ**h**_{1} whose elements are

The state-specific nonequilibrium PEqS method is summarized in Algorithm 2.

1: Proceed with Algorithm 1 and save data to disk. |

2: Reference data: E_{0}, G_{elst,0}, $\phi pol,0slow$(r), and $\rho pol,0slow$(r) |

3: Initialize Δh_{1} = 0 |

4: for j = 1, 2, … do until SCF error < τ_{SCF} |

5: Diagonalize $F(j)=F0(j)+\Delta h1(j)$ to obtain P^{(j)} |

6: Compute ρ_{sol,1}(r) and φ_{sol,1}(r) |

7: if j = 1 then |

8: $\rho tot,1fast(r)=\rho sol,1(r)+\rho pol,0slow$(r) |

9: $\phi tot,1fast(r)=\phi sol,1(r)+\phi pol,0slow$(r) |

10: else |

11: $\rho tot,1fast(r)=\rho sol,1(r)+\rho pol,0slow(r)+\rho pol,1fast(r)$ |

12: $\phi tot,1fast(r)=\phi sol,1(r)+\phi pol,0slow(r)+\phi pol,1fast$(r) |

13: end if |

14: Initialize $\rho iter,1(0)$(r) using $\phi tot,1fast,(0)$(r) |

15: for k = 0, 1, … do until Δρ_{iter,1} < τ_{PEqS} |

16: Compute $\phi tot,1fast,(k+1)$(r) |

17: Update $\rho iter,1(k+1)$(r), $\rho pol,1fast$(r), and $\rho tot,1fast$(r) |

18: end for |

19: Update $\Delta h1(j)$ |

20: Compute E_{1} = E_{int,1} + G_{elst,1} |

21: end for |

22: Compute VIE = E_{1} − E_{0} |

1: Proceed with Algorithm 1 and save data to disk. |

2: Reference data: E_{0}, G_{elst,0}, $\phi pol,0slow$(r), and $\rho pol,0slow$(r) |

3: Initialize Δh_{1} = 0 |

4: for j = 1, 2, … do until SCF error < τ_{SCF} |

5: Diagonalize $F(j)=F0(j)+\Delta h1(j)$ to obtain P^{(j)} |

6: Compute ρ_{sol,1}(r) and φ_{sol,1}(r) |

7: if j = 1 then |

8: $\rho tot,1fast(r)=\rho sol,1(r)+\rho pol,0slow$(r) |

9: $\phi tot,1fast(r)=\phi sol,1(r)+\phi pol,0slow$(r) |

10: else |

11: $\rho tot,1fast(r)=\rho sol,1(r)+\rho pol,0slow(r)+\rho pol,1fast(r)$ |

12: $\phi tot,1fast(r)=\phi sol,1(r)+\phi pol,0slow(r)+\phi pol,1fast$(r) |

13: end if |

14: Initialize $\rho iter,1(0)$(r) using $\phi tot,1fast,(0)$(r) |

15: for k = 0, 1, … do until Δρ_{iter,1} < τ_{PEqS} |

16: Compute $\phi tot,1fast,(k+1)$(r) |

17: Update $\rho iter,1(k+1)$(r), $\rho pol,1fast$(r), and $\rho tot,1fast$(r) |

18: end for |

19: Update $\Delta h1(j)$ |

20: Compute E_{1} = E_{int,1} + G_{elst,1} |

21: end for |

22: Compute VIE = E_{1} − E_{0} |

### C. Dielectric function

#### 1. Bulk environment

In conventional PCM calculations, the solute cavity is a two-dimensional surface constructed from a union of atom-centered spheres, possibly with additional surface elements added to smooth the seams where those spheres intersect. In any case, it is assumed that the dielectric constant changes abruptly at the cavity surface, switching from its vacuum value (*ε* = 1) inside the cavity to a value appropriate for bulk solvent (*ε*_{solv}) outside. This abrupt change in *ε* poses no problems within the PCM formalism but is problematic in the present context, where it is necessary to discretize three-dimensional space. As such, the sharp transition in *ε*(**r**) must be smoothed.^{61}

Several groups have proposed dielectric functions that are functionals of the electron density and thus conform automatically to molecular shape,^{40,62–64} analogous to using an isodensity contour to define the cavity in a PCM calculation.^{65–67} Such cavities (or dielectric functions) must be self-consistently updated at each self-consistent field (SCF) cycle. Instead, we choose the rigid cavity model of Ref. 41 that uses a product of spherically symmetric atomic switching functions *s*_{α} to smooth the discontinuous function *ε*(**r**) that is used (implicitly) in PCM calculations and is based on atom-centered spheres. The resulting dielectric function is

For generality, we have written this in terms of an arbitrary “vacuum” dielectric constant *ε*_{vac} inside the cavity. For any choice *ε*_{vac} ≠ 1, however, the electronic structure program ought properly to be modified to use a Coulomb operator $(\epsilon vacr)\u22121$ rather than *r*^{−1}. All numerical calculations presented here use *ε*_{vac} = 1.

The switching functions in Eq. (2.34) are defined as

where *d*_{α} is the radius of the atomic sphere centered at **R**_{α}. With *s*_{α} chosen in this way, the dielectric function transitions smoothly from *ε*_{vac} to *ε*_{solv} over a region whose length is ≈4Δ and is centered at a distance *d*_{α} from nucleus *α*. Following Ref. 41, we set Δ = 0.265 Å, and following standard PCM convention we take *d*_{α} = 1.2 *r*_{vdW,α},^{19,22,56} where the atomic van der Waals (vdW) radii *r*_{vdW,α} are taken from Bondi’s set,^{68} except that for hydrogen we use the updated value of 1.1 Å.^{69} As such, Eq. (2.34) mimics the dielectric function that is used implicitly in PCM calculations, except that the former is continuous everywhere.

However, this dielectric function poses a problem when explicit solvent molecules are included as part of the solute. An egregious example is the case of the hydrated electron, *e*^{−}(aq),^{70} represented in the following example as a $(H2O)28\u2212$ cluster model extracted from a condensed-phase simulation.^{71} As shown in Fig. 1, this cluster model consists of approximately two solvation shells of water molecules coordinated around an unpaired electron. Cluster models of this type, combined with a PCM to capture long-range solvation effects, have previously been used to estimate the VIE of *e*^{−}(aq),^{72} but this is potentially problematic because the vdW cavity that is conventionally used in PCM calculations places high-dielectric regions in between water molecules.

This can be seen explicitly by plotting the dielectric function *ε*(**r**) in Eq. (2.34), for a vdW cavity corresponding to the water configuration shown in Fig. 1. (We emphasize that up to a switching function to smooth the transition between *ε* = 1 and *ε* = 78, this is the dielectric function that is used, implicitly, in PCM calculations.) A two-dimensional slice through *ε*(**r**) is plotted in Fig. 2(a). Although the cavity correctly conforms to the molecular shape of the water cluster, the dielectric function is problematic in the region near the cluster’s center of mass (c.o.m.). There, blue and green regions indicate a solvent-like value of *ε* that penetrates into the region of space occupied by the unpaired electron that, as part of the solute, ought instead to experience *ε* = 1. While *e*^{−}(aq) might seem like an unusual case due to the esoteric nature of the solute, the problem is a general one, as illustrated by the dielectric function for a $F\u2212(H2O)31$ cluster that is plotted in Fig. 3. High-dielectric regions can once again be found inside of the atomistic QM region.

To address this problem, we pursue an approach used also in the context of PCMs, in which a fictitious spherical “solvent probe” is rolled along the surface of the vdW cavity (constructed from *unscaled* Bondi radii); the locus of points traced out by the center of this probe sphere defines the solvent-accessible surface (SAS).^{73} Equivalently, the SAS is simply a vdW surface constructed using radii *d*_{α} = *r*_{vdW,α} + *r*_{probe} that are equal to vdW (Bondi) radii augmented by the probe radius. For aqueous solvation, the standard choice is *r*_{probe} = 1.4 Å,^{19,73} representing half the distance to the first peak in the oxygen–oxygen radial distribution function of liquid water.^{74} Using *d*_{α} = *r*_{vdW,α} + *r*_{probe} in Eq. (2.34) successfully removes values *ε* > 1 in the interstices between water molecules, however, the resulting VIEs are quite poor and in some cases the PEqS procedure is difficult to converge. A “modified” SAS (mSAS) construction, using the reduced value *r*_{probe} = 0.7 Å, alleviates the convergence problems and affords more reasonable VIEs but does not entirely eliminate artifactual high-dielectric regions between water molecules, as shown in Fig. 2(b).

To rectify this, we introduce a “hybrid” cavity that retains the conformity to molecular shape exhibited by the vdW and SAS cavities but eliminates problematic high-dielectric regions in this *e*^{−}(aq) test case. The hybrid cavity is built upon the mSAS cavity (*r*_{probe} = 0.7 Å), adding a geometric constraint that exploits the roughly spherical nature of the solute configurations to ensure that the dielectric function assumes the value *ε* = 1 inside the solute region. To this end, we adapt a procedure from Ref. 75 that was used to characterize binding motifs of excess electrons in $(H2O)n\u2212$ clusters. The shape of each solute configuration (water cluster) is approximated as an ellipsoid centered at the cluster c.o.m. (*x*_{0}, *y*_{0}, *z*_{0}), whose surface is defined by the equation *S*(*x*, *y*, *z*) ≡ 1, where

The volume enclosed by *S*(*x*, *y*, *z*) is treated as the solvent-excluded region, and the hybrid cavity, whose dielectric function is plotted in Fig. 2(c), is constructed from a mSAS cavity by enforcing the condition that *ε*(**r**) = *ε*_{vac} if *S*(*x*, *y*, *z*) < 1. The parameter *a* is set equal to the maximum atomic-to-c.o.m. distance along the *x* axis, plus a distance *d*_{α} − 2Δ that centers the switching function in Eq. (2.35) at a distance *d*_{α} from the nucleus. This furthermore ensures that *ε*(**r**) ≈ *ε*_{vac} at a distance 2Δ from any atomic center. The parameters *b* and *c* are defined similarly, for the *y* and *z* directions.

Figure 2(c) illustrates the hybrid cavity dielectric function for $(H2O)28\u2212$ that is obtained using this procedure. This model affords a satisfactory description of the dielectric environment, in the sense that high-dielectric regions between water molecules are eliminated. As such, we use this definition of the cavity and corresponding dielectric function for all PEqS calculations reported in Sec. V.

#### 2. Interfacial environment

The dielectric function for the liquid/vapor interface is defined as in our previous work.^{42} Specifically, we interpolate *ε*(**r**) from *ε*_{solv} = 78.4 to *ε*_{vac} = 1.0 across the Gibbs dividing surface (GDS). The periodic slabs used for molecular dynamics (MD) simulations at the interface extend infinitely in the *x* and *y* directions, and the location *z*_{GDS} of the GDS is determined over the course of the simulation by computing ensemble-averaged solvent density profiles. These are computed individually, for each solute, using 0.5 Å bins along the *z* direction, and the resulting density profiles are then fit to the following functional form:^{32,76,77}

Here, *ρ*_{solv} is the bulk liquid density (treated here as a fitting parameter) and *α* is a parameter such that the thickness of the liquid/vapor interface is ≈4/*α*. The hyperbolic tangent term in Eq. (2.37) is positive if *z* > *z*_{GDS} and negative if *z* < *z*_{GDS}.

Best-fit parameters *ρ*_{solv}, *α*, and *z*_{GDS} are listed in Table I for each solute considered in this work. Fitted values of *ρ*_{solv} are in reasonable agreement with the actual density of liquid water at 298 K. Fitted values of *α* demonstrate that the interfacial region for the ionic solutes is discernibly thicker than the liquid/vapor interface for neat water, an observation that is also reported in other studies of anions at interfaces.^{77,78} Taking parameters from Table I, we describe the *z*-dependence of the interfacial dielectric function as in previous work,^{32,42,76}

Figure 4 illustrates the dielectric function *ε*(**r**) for a $(H2O)24\u2212$ cluster extracted from an interfacial configuration of *e*^{−}(aq),^{79} with the c.o.m. placed at the origin.

Solute . | ρ_{solv} (g/cm^{3})
. | α (Å^{−1})
. | z_{GDS} (Å)
. |
---|---|---|---|

Li^{+} | 0.976 | 0.652 | −9.205 |

Na^{+} | 0.979 | 0.604 | −9.154 |

H_{2}O | 0.986 | 0.724 | −8.963 |

e^{−} | 1.016 | 0.668 | −1.508 |

F^{−} | 0.976 | 0.635 | −9.105 |

Cl^{−} | 0.969 | 0.626 | −9.328 |

Solute . | ρ_{solv} (g/cm^{3})
. | α (Å^{−1})
. | z_{GDS} (Å)
. |
---|---|---|---|

Li^{+} | 0.976 | 0.652 | −9.205 |

Na^{+} | 0.979 | 0.604 | −9.154 |

H_{2}O | 0.986 | 0.724 | −8.963 |

e^{−} | 1.016 | 0.668 | −1.508 |

F^{−} | 0.976 | 0.635 | −9.105 |

Cl^{−} | 0.969 | 0.626 | −9.328 |

## III. NUMERICAL SOLUTION OF POISSON’S EQUATION

Solution of the linear equations that define PCMs is often (though not always^{22,80–82}) accomplished via matrix inversion. Matrix diagonalization incurs a cost that is $O(Ngrid3)$ in floating-point operations and $O(Ngrid2)$ in memory, for *N*_{grid} discretization points, and for PCMs this is typically trivial in comparison to the cost of the QM calculation for the solute. An exception is QM/MM/PCM calculations, where the QM/MM solute is potentially quite large, and conjugate gradient (CG) procedures have been developed to handle such cases.^{22,80} For the PEqS method, however, direct inversion is prohibitively expensive from the start, as ∼10^{6} Cartesian grid points might be required to discretize three-dimensional space, with a memory cost alone that would exceed 7 Tb to store the discretized Laplacian. It is therefore essential to employ relaxation techniques such as iterative CG procedures. Here, we discuss finite-difference discretization schemes for solving Poisson’s equation on large Cartesian grids and also discuss improving the efficiency of PEqS using a multigrid method. Much of this work, including the multigrid method, is new since the pilot implementation of PEqS that was reported in Ref. 42.

### A. Finite-difference scheme

where the factor of −4*π* that ordinarily appears in Poisson’s equation is instead included in *ρ*(**r**). Writing out the Laplacian operator explicitly, this is

subject to the Dirichlet boundary condition

at the surface boundary *δ*Ω (see below).

For a uniform rectangular grid centered at the origin, with side lengths {*L*_{x}, *L*_{y}, *L*_{z}} containing {*N*_{x}, *N*_{y}, *N*_{z}} grid points (so that *N*_{grid} = *N*_{x}*N*_{y}*N*_{z}), the domain Ω is defined as

The surface boundary is defined by the collection of rectangular planes

For convenience, we assume in what follows that the grid is cubic, with equal spacing *h* in each direction. Cartesian coordinates are then mapped onto grid coordinates as *x*_{i} = −*L*_{x}/2 + *ih*, where *i* = 0, …, (*N*_{x} − 1).

The value of the electrostatic potential *φ*(*x*, *y*, *z*) at the grid point (*x*_{i}, *y*_{j}, *z*_{k}) is denoted as

with a similar notation for other discretized quantities. Expressions for the discretized first and second derivatives of *φ*_{i,j,k} are obtained using a central finite-difference scheme. A general expression for an *n*th-order approximation to the *m*th-order derivative, whose finite-difference approximation exhibits error of $O(h2n)$, is

for certain coefficients *c*_{m,p}. We use an eighth-order (*n* = 4) finite-difference approximation for the first (*m* = 1) and second (*m* = 2) derivatives. Coefficients *c*_{m,n} for this approximation are given in Table II.

### B. Multigrid approach

We employ a multigrid method to improve the efficiency of the iterative CG procedure. To facilitate the following discussion, let us rewrite Poisson’s equation [Eq. (3.1)] in a discretized form involving a matrix–vector product,

The *N*_{grid} × *N*_{grid} matrix $Lh$ contains the discretized Laplacian operator, and the vectors *φ*^{h} and *ρ*^{h} contain the discretized values $\phi i,j,k,h$ and $\u22124\pi \rho i,j,k,h$, respectively. The quantities *φ*^{h} and *ρ*^{h} are subject to appropriate boundary conditions, and the superscript *h* signifies the level (fineness) of the discretization. The CG algorithm avoids the prohibitive memory cost associated with forming $Lh$ explicitly; only its action on the vector *φ*^{h} is required.

Via Fourier analysis of the discretization errors

it has been shown that the spectrum of errors contains wavelengths *λ* that are comparable to, or larger than, the grid resolution, *h*.^{83} The CG routine efficiently eliminates discretization errors with *λ* ≃ *h* but struggles with error components for which *λ* > *h*.^{83} Iterative techniques such as CG therefore effectively smooth out short-wavelength discretization errors, but they do not perform well (or at least, efficiently) for obtaining a fully converged solution due to long-wavelength error components. To illustrate this, the CG routine was employed to compute *φ*^{h} in Eq. (3.8) for a single H_{2}O molecule placed at the center of a (15 Å)^{3} cubic Cartesian grid with a resolution *h* = 0.074 Å. Figure 5 shows the Euclidean norm of the residual error in the electrostatic potential,

as a function of iteration number. There is a rapid decrease in ∥**r**^{h}∥ during the first few iterations, but in the end more than 250 iterations are required to achieve convergence, defined as ∥**r**^{h}∥ < 10^{−5} a.u. The inability of the CG routine to eliminate the long-wavelength error components furthermore manifests as a broad and slowly decaying shoulder, evident in iterations 10–110 in Fig. 5.

Since the charge density *ρ*(**r**) computed by the electronic structure calculation is known essentially to arbitrary accuracy, at least compared to the ∼10^{−5} a.u. residual convergence error in the CG approach, we can safely assume that *ρ*(**r**) is free of discretization error upon formation of *ρ*^{h}. Supposing that Eq. (3.8) can be solved exactly for the electrostatic potential, then $\rho h=Lh\phi exact$. Thus Eq. (3.10) can be written in terms of **v**^{h} [defined in Eq. (3.9)], even when the exact solution is not known,^{83}

Despite never actually acquiring *φ*^{exact}, Eq. (3.11) provides an avenue for computing the quantity **v**^{h}, which is an integral part of the multigrid method.

The multigrid method seeks to obviate undesired computational effort spent eliminating long-wavelength error components that results in the slow convergence exhibited by the CG routine. The ultimate goal is a solution to Eq. (3.8) on a fine rectangular grid, and we refer to this as the “target” grid resolution, denoted by *h*. The schematic for a two-level “V-cycle” multigrid approach^{83} is illustrated in Fig. 6. (The nomenclature is explained below, where we introduce an alternative “W-cycle” approach as well.) The multigrid method is designed to relax the iterative solution on the target grid, where the computational cost is highest, as few times as possible. This is step 1 in the algorithm outlined in Fig. 6, and it serves to reduce short-wavelength errors. The resulting residual **r**^{h} obtained using the target grid is then *restricted* to a grid with only half as many grid points in each Cartesian direction. The resolution of this grid is denoted as *H* = 2*h*, and the iterative solution on the coarser grid serves to reduce longer-wavelength components of the error. The restriction **r**^{h} → **r**^{H} is illustrated in step 2 of Fig. 6 and is accomplished using a restriction matrix $IhH$,

In practice, $IhH$ is not formed and its action on **r**^{h} to generate **r**^{H} is expressed as

The notation (*I*, *J*, *K*) in $rI,J,KH$ is introduced to denote that a different mapping scheme for the coarse-grid coordinates is required, namely, *x*_{I} = −*L*_{x}/2 + *IH*_{x} for *I* = 0, …, (*N*_{x} − 1)/2 and *H*_{x} = 2*h*_{x}. Equation (3.13) is valid for three dimensions and shows that a coarse grid point takes its value from all neighboring points on the target grid, with a weight determined by its proximity to $rI,J,KH$.

After forming **r**^{H}, a Poisson-like equation is solved to obtain the coarse grid discretization error, **v**^{H},

This is illustrated in step 3 of Fig. 6. Relaxing **v**^{H} on the coarse grid reduces the problematic long-wavelength error components with which the primitive CG routine struggles; consequently, **v**^{H} provides a better correction to *φ*^{h}. However, **v**^{H} cannot be used directly to correct *φ*^{h} because the former is defined only on the coarser grid and must be interpolated back to the target grid. This process of inverse restriction is illustrated in step 4 of Fig. 6. Interpolation of **v**^{H} to form **v**^{h} is accomplished via the inverse iteration matrix $IHh$,

The action of $IHh$ is expressed by the following set of equations:

The interpolated discretization error **v**^{h} corrects *φ*^{h} according to

This is step 5 of Fig. 6. Convergence is achieved when $\Vert rnewh\u2212roldh\Vert $ falls below a desired threshold, at which point *φ*^{h} is the fully relaxed solution to Eq. (3.8). Otherwise the process shown in Fig. 6 repeats with step 2.

In contrast to the two-level algorithm outlined in Fig. 6, the PEqS method implemented here actually uses a four-level method that employs two additional coarse grids whose resolutions are *H* = 4*h* and *H* = 8*h*. These four-level schemes improve upon the two-level scheme by eliminating error components on multiple length scales, affording more effective corrections at each grid level and resulting in a fully relaxed, target-grid solution in fewer iterations.^{83} Schematics for “V-cycle” and “W-cycle” variants of the four-level method are provided in Fig. 7, which introduces parameters *γ*_{0}, *γ*_{1}, *γ*_{2}, and *γ*_{3} that control the maximum number of CG iterations spent at various grid levels. The strategy is to always fully relax the error vector **v**^{H=8h} at grid level 3 (the coarsest grid), so we set the parameter *γ*_{0} equal to the maximum number of allowed iterations, *γ*_{0} = 500 here. The other parameters are set to *γ*_{1} = 2, *γ*_{2} = 3, and *γ*_{3} = *γ*_{1} + *γ*_{2}, as in Ref. 83. Passing through either cycle outlined in Fig. 7, one performs CG iterations of the equation $Lnhvnh=rnh$ in order to compute **v**^{nh}, for *n* = 1, 2, 4, or 8 as appropriate. Iterations continue until the residual **r**^{nh} is reduced below threshold or until the maximum number of iterations is reached. At grid levels finer than *H* = 8*h*, this means that **v**^{nh} need not be fully relaxed at each step. (Convergence failure on the coarsest grid, however, implies the failure of the whole algorithm, though we have not found this to be a problem with the parameters described herein.) Proceeding in this way, the solution *φ*^{h} on the target grid, which is the most expensive to compute, is relaxed a total of *γ*_{1} + *γ*_{2} times, in either the V- or W-cycle approach.

In more detail, the four-level algorithms proceed as follows. At the target grid level, the residual error **r**^{h} corresponding to *φ*^{h} after *γ*_{1} CG iterations is restricted to grid level 1, signified by a downward arrow in Fig. 7. Using **r**^{H=2h} in Eq. (3.14), **v**^{H=2h} is relaxed *γ*_{1} times and is stored in memory. This process of restriction is repeated for each downward arrow until reaching the coarsest level of discretization (grid level 3), where the solution is then fully relaxed. Upward arrows in Fig. 7 signify interpolation of the discretization error from a coarser grid to a finer one using Eq. (3.15), and the resulting error vector is then used to update the residual error at grid levels 2 and 1, or the electrostatic potential on the target grid. Each of these is relaxed *γ*_{2} times. Returning to the example in Fig. 5, we note that the four-level methods require about 75 iterations (V-cycle) and 35 iterations (W-cycle) on the target grid, a remarkable improvement over the standard CG routine. We select the four-level W-cycle method for PEqS calculations due to its superior performance in this example.

## IV. COMPUTATIONAL METHODS

In this work, we compute VIEs for neat liquid water and for five aqueous ions: F^{−}, Cl^{−}, Li^{+}, Na^{+}, and the hydrated electron. Experimental VIEs are available for each of these species,^{11,15} and while we examine *e*^{−}(aq) due mainly to our group’s long-standing interest in this species,^{42,70,75,84–87} and because preliminary calculations on *e*^{−}(aq) were previously reported using a perturbative version of PEqS,^{42} the four atomic ions are selected because their simple structure should eliminate any concerns about adequate MD sampling. This is especially true given that an accurate polarizable force field (amoeba) is available for these species.^{88,89} This allows us to focus on the role of the continuum model in VIE calculations.

### A. Molecular dynamics simulations

Simulations of neat liquid water were performed with 222 water molecules in a periodic cell 18.8 Å on a side, corresponding to a density of 0.9995 g/cm^{3} at *T* = 300 K, using the amoeba force field^{88,89} as implemented in the tinker software package,^{90} v. 7.1.2. Electrostatic interactions were computed using Ewald summation with a real-space cutoff of 9.4 Å. The neat liquid water simulations were equilibrated for 1 ns with the final 500 ps extracted for further use. For simulations of the neat liquid/vapor interface, the final snapshot from the bulk water simulation is used as a starting point and the simulation box was extended to 90.0 Å in the *z* direction so that the dimensions of the simulation cell measured 18.8 Å × 18.8 Å × 90.0 Å. The resulting water “slab” was equilibrated for an additional 1 ns at *T* = 300 K.

MD simulations for the aqueous halides and alkali cations were initialized starting from the equilibrated neat liquid water simulation, replacing either the H_{2}O molecule nearest to the center of the cell (in the bulk simulations) or that nearest to the interface (in slab simulations) with an ion. The bulk simulations were then equilibrated for an additional 250 ps at *T* = 300 K followed by a 500 ps production run with snapshots stored every 5 ps. In contrast, simulations at the liquid/vapor interface were not equilibrated, and a 500 ps production run began immediately after insertion of the ion, again with a stride of 5 ps between saved snapshots. The snapshots for *e*^{−}(aq) in liquid water and at the air/water interface were taken from QM/MM simulations reported in Refs. 71 and 79. They are the same snapshots used in some of our previous studies of *e*^{−}(aq).^{42,87}

Solute configurations for subsequent PEqS and PCM calculations were generated from the stored snapshots by selecting a QM region that includes all water molecules within a sphere of radius 5.5 Å centered at the ion, or in the case of *e*^{−}(aq), centered at the centroid of the spin density. Extensive convergence tests in previous work suggest that larger QM regions are unnecessary for VIE calculations.^{42} For neutral water, the H_{2}O molecule nearest to the center of the simulation cell is chosen as the center for the bulk liquid configurations, whereas the water molecule nearest to the GDS is chosen for the liquid/vapor configurations. VIEs reported here are averages over 101 snapshots, each separated in time by 5 ps, except in the case of *e*^{−}(aq) where we use 87 snapshots, each separated in time by 100 fs.

### B. Continuum solvation models

The dielectric function, charge densities, and electrostatic potentials required for PEqS calculations were discretized on a (25 Å)^{3} Cartesian grid with spacing Δ*x* = Δ*y* = Δ*z* = 0.24 Å. The hybrid cavity model described in Sec. II C 1 is used to construct *ε*(**r**), and for the interfacial configurations the dielectric function is modified according to Eq. (2.38), with solute-specific parameters taken from Table I. Concurrent acquisition of the polarization response charge density and electrostatic potential through Eqs. (2.10) and (2.25) is accomplished using the four-level W-cycle multigrid technique (Sec. III B), with relaxation parameters *γ*_{0} = 500, *γ*_{1} = 2, *γ*_{2} = 3, and *γ*_{3} = *γ*_{1} + *γ*_{2}. The iterative charge density is updated until the residual Δ*ρ*_{iter} falls below a threshold *τ*_{PEqS} = 10^{−5} a.u.

To complement the PEqS calculations, we also compute VIEs of the bulk water configurations using a nonequilibrium version^{48–50} of IEF-PCM,^{91} the “integral equation formalism” (IEF) version of PCM.^{22,91–93} The solute cavity in these calculations is defined in one of the two ways. One approach uses a single sphere to encapsulate the QM region; this region was carved out of the solution using a radius of 5.5 Å and we set the spherical cavity radius to 7.5 Å. This sphere is then discretized using a Lebedev integration grid with 5294 points. Alternatively, we use the SAS cavity constructed from a union of atom-centered spheres with radii *r*_{vdW} + *r*_{probe}, where *r*_{vdW} is an unscaled Bondi radius^{68,69} and *r*_{probe} = 1.4 Å. (This is the standard SAS definition,^{19,73} not the modified one used in Sec. II C 1 to define the hybrid cavity.) In this case, each atomic sphere is discretized with a 302-point Lebedev grid using the switching/Gaussian algorithm.^{94}

For isotropic solvation in bulk water, we expect that the PEqS and PCM methods should afford similar results, up to discretization errors that can be made arbitrarily small,^{22,91} and neglecting charge penetration effects (i.e., volume polarization^{95,96}) arising from the part of the solute’s charge density that extends beyond the cavity. The latter effects are mitigated within the IEF-PCM framework,^{92,97,98} and through the inclusion of explicit water molecules around the anions. Unfortunately, the PEqS method in its current implementation is sensitive to the Gaussian width parameter *σ* that is used to blur the nuclear charges [Eq. (2.8)], so we exploit the expected numerical equivalence with IEF-PCM to set this parameter. Setting *σ* = 0.300 Å (F^{−}), 0.525 Å (H_{2}O), or 0.570 Å (Li^{+}) affords PEqS and IEF-PCM solvation energies that are in good agreement, for a few test configurations. The F^{−} value of *σ* was used for all three anions and the Li^{+} value was used for both cations. Values of *σ* determined in bulk water were used also for the interfacial PEqS calculations.

### C. Electronic structure calculations

The state-specific PEqS method has been implemented in a locally modified version of the q-chem program^{99} and will be released with v. 5.1. Electronic structure calculations were performed at the level of second-order Møller-Plesset perturbation theory (MP2) within the resolution-of-identity (RI) approximation.^{100} All electrons were correlated for the Li^{+}(aq) and Na^{+}(aq) calculations, whereas other calculations use a frozen core. The PEqS part of the calculation resides in the Hartree-Fock iterations and uses the Hartree-Fock electron density. (Nonequilibrium PCM results for excitation energies suggest that a fully self-consistent use of the correlated electron density has negligible effect on the results.^{50}) The SCF convergence threshold is set to *τ*_{SCF} = 10^{−5} a.u. in all calculations, with an integral and shell-pair screening threshold of 10^{−8} a.u.

We use the 6-311+G* basis set for all H_{2}O molecules, except in the case of *e*^{−}(aq) where we use 6-311++G* instead, to ensure that the interstices between the water molecules are supported by basis functions. (We have previously concluded that one set of diffuse functions on all atoms is sufficient to support a hydrated electron that occupies an excluded volume in the structure of liquid water.^{84,86}) For Li and Na, we use the cc-pCVTZ basis set, which includes functions to describe core/valence correlation, whereas for F and Cl we use aug-cc-pVTZ. In all cases, we employ the auxiliary basis set designed for either cc-pVTZ or aug-cc-pVTZ,^{101} as appropriate.

The valence photoelectron spectrum of liquid water consists of a broad absorption feature centered at 11.3 eV,^{102} attributed to ionization of a 1*b*_{1} MO localized on a single H_{2}O molecule.^{11,102–105} Experiments to determine the VIE of F^{−}(aq) are complicated by the fact that the fluoride signal is embedded in the 1*b*_{1} band of water.^{11,16} When explicit water molecules are included in the QM region, the corresponding VIE calculation is problematic as well, not only for F^{−}(aq) but also for Na^{+}(aq), Li^{+}(aq), or any solute whose VIE lies near or above that of water. In such cases, a simple calculation of the lowest-energy state of the ionized system results in ionization of water rather than the desired solute,^{18} as shown for Li^{+}(aq) in Fig. 8(a) and for F^{−}(aq) in Fig. 8(b). In these examples, the VIE is computed for a system consisting of the atomic ion surrounded by about 30 explicit water molecules and spherical PCM boundary conditions. In both cases, however, it is a water molecule that is ionized rather than the atomic ion. Figure 8(c) shows that even neat liquid water is problematic, as in this case the lowest-energy ionized system places the cation hole on an “edge” water molecule that lies near the QM/continuum boundary and does not fully participate in the hydrogen-bonding network. Computed VIEs are in very poor agreement with experiment, e.g., 8.3 eV for neat liquid water.

To ionize F^{−}, Cl^{−}, Li^{+}, or Na^{+} embedded in a water cluster, one must remove an electron from an orbital other than the HOMO of the full QM system. At the same time, one still wants to include orbital relaxation effects upon ionization yet prevent variational collapse of the wave function to the lowest-energy solution, which is the one depicted in Figs. 8(a)–8(c). The maximum overlap method^{106} (MOM) was designed precisely for this purpose and has been used, for example, to compute core-excited states (K-edge spectra) by moving an electron from a core orbital into an extra-valence (virtual) orbital and then relaxing the orbitals while searching for the maximum overlap solution.^{107}

Here, we start from an initial guess generated using the so-called fragment MO (FragMO) procedure^{108} and then use the MOM method to preserve the character of the initial orbitals during subsequent SCF iterations. The FragMO procedure first computes MOs on isolated fragments (here, either H_{2}O or an atomic ion), which affords an easy means to remove an electron from an MO associated with a particular fragment, e.g., the 1*b*_{1} orbital of a *particular* water molecule. A superposition of fragment density matrices is then used as the initial guess for the supersystem SCF calculation. The converged SCF solution obtained from this FragMO/MOM procedure contains a “hole,” which manifests as a single virtual orbital with an energy below that of the HOMO. For Li^{+}(aq), the core hole occasionally leads to accidental quasi-degeneracies amongst the Hartree-Fock eigenvalues such that the MP2 energy denominator becomes very small. For this species only, we therefore omit this “hole” orbital from the MP2 calculation of the correlation energy.

VIEs recomputed using nonequilibrium IEF-PCM in conjunction with this FragMO/MOM SCF procedure are shown in Figs. 8(d)–8(f). It is first of all clear that this approach successfully ionizes the desired species, which is a centrally located H_{2}O molecule in the case of neat liquid water. The computed VIEs are also in far better agreement with experiment, e.g., for the snapshots shown in Fig. 8 they are 11.8 eV (computed) for H_{2}O(aq) versus 11.3 eV (experiment);^{102} 11.3 eV versus 11.6 eV for F^{−}(aq);^{11} and 62.0 eV versus 60.4 eV for Li^{+}(aq).^{18} We use the FragMO/MOM SCF procedure for all systems except *e*^{−}(aq), where the VIE lies well below that of liquid water and a standard SCF procedure is adequate.

It is worth emphasizing that the use of the FragMO/MOM procedure does not force the converged, singly occupied MO to be localized on any single monomer; only the initial-guess MOs are localized and in subsequent SCF iterations the orbitals are free to delocalize, as seen, for example, in Fig. 8(e). This is likely to be important, e.g., to simulate the full valence photoelectron spectrum of liquid water. Although the lowest-energy feature is assigned to ionization of the 1*b*_{1} orbital localized on a single water molecule, the higher-energy 3*a*_{1} feature is broadened into an (unresolved) doublet,^{11,102–105} assigned to ionization from bonding and anti-bonding combinations of the 3*a*_{1} orbitals on two hydrogen-bonded H_{2}O molecules.^{105,109} In this work, we consider only the lowest VIE of each solute, but in principle it should be possible to simulate the splitting of the 3*a*_{1} feature by constructing the appropriate initial-guess orbitals on a water dimer.

## V. RESULTS

### A. VIEs in bulk water

Configurationally averaged VIEs in bulk liquid water, computed using the nonequilibrium PEqS and PCM methods, are shown in Table III along with experimental values obtained from liquid microjet experiments.^{11,15,18,102} Also listed is the average number ⟨*N*⟩ of explicit water molecules in the QM region, which is about 30 in each case, amounting to roughly two solvation shells. In previous work on *e*^{−}(aq),^{42} we observed that increasing ⟨*N*⟩ up to 90, corresponding to an increase in the radius of the QM region from 5.5 Å to 8.0 Å, changed VIEs by ≤0.1 eV, both in bulk liquid water and at the liquid/vapor interface. Amongst the solutes considered here, we anticipate that *e*^{−} is most acutely affected by the size of the QM region and thus we regard the present calculations to be adequately converged in this respect.

. | . | . | Computed VIE (eV) . | ||
---|---|---|---|---|---|

. | . | Experimental . | PEqS . | PCM . | |

Solute . | ⟨N⟩
. | VIE (eV) . | Hybrid . | Spherical . | SAS . |

Li^{+} | 30 | 60.40 ± 0.07^{a} | 61.41 ± 0.45 | 61.85 ± 0.45 | 61.59 ± 0.41 |

Na^{+} | 29 | 35.40 ± 0.04^{a} | 36.09 ± 0.43 | 36.54 ± 0.44 | 36.34 ± 0.40 |

H_{2}O | 30 | 11.31 ± 0.04^{b} | 11.53 ± 0.51 | 11.64 ± 0.51 | 11.55 ± 0.45 |

e^{−} | 30 | 3.7 ± 0.1^{c} | 3.75 ± 0.55 | 3.16 ± 0.32 | 3.18 ± 0.28 |

F^{−} | 30 | 11.58^{d} | 11.66 ± 0.36 | 11.37 ± 0.39 | 11.52 ± 0.37 |

Cl^{−} | 32 | 9.60 ± 0.07^{a} | 9.65 ± 0.37 | 9.36 ± 0.38 | 9.41 ± 0.35 |

. | . | . | Computed VIE (eV) . | ||
---|---|---|---|---|---|

. | . | Experimental . | PEqS . | PCM . | |

Solute . | ⟨N⟩
. | VIE (eV) . | Hybrid . | Spherical . | SAS . |

Li^{+} | 30 | 60.40 ± 0.07^{a} | 61.41 ± 0.45 | 61.85 ± 0.45 | 61.59 ± 0.41 |

Na^{+} | 29 | 35.40 ± 0.04^{a} | 36.09 ± 0.43 | 36.54 ± 0.44 | 36.34 ± 0.40 |

H_{2}O | 30 | 11.31 ± 0.04^{b} | 11.53 ± 0.51 | 11.64 ± 0.51 | 11.55 ± 0.45 |

e^{−} | 30 | 3.7 ± 0.1^{c} | 3.75 ± 0.55 | 3.16 ± 0.32 | 3.18 ± 0.28 |

F^{−} | 30 | 11.58^{d} | 11.66 ± 0.36 | 11.37 ± 0.39 | 11.52 ± 0.37 |

Cl^{−} | 32 | 9.60 ± 0.07^{a} | 9.65 ± 0.37 | 9.36 ± 0.38 | 9.41 ± 0.35 |

Agreement with experimental VIEs is excellent for the anionic solutes and for H_{2}O, especially for the PEqS treatment of solvation. Although we do not expect exact agreement between the PEqS and PCM calculations, primarily because the solute cavities differ but also due to the approximate manner in which IEF-PCM accounts for volume polarization,^{95–97} results from all three solvation models agree to within about 0.4 eV. (In fairness, the Gaussian widths for the PEqS nuclear charges were determined in order to match IEF-PCM solvation energies for a few snapshots.) For *e*^{−}(aq), both PCM methods underestimate the VIE by 0.5 eV, whereas the PEqS calculations are spot on, at 3.7 eV; arguably, this is the system for which PCM boundary conditions are most questionable^{110} due to the delocalized nature of the solute. For reasons that are not clear, computed VIEs for Li^{+}(aq) and Na^{+}(aq) exhibit much larger errors of 0.7–1.0 eV (PEqS) or 0.9–1.4 eV (PCM) with respect to experiment.

Uncertainties on the calculated VIEs in Table III represent one standard deviation across MD snapshots and provide an estimate of the inhomogeneous broadening arising from thermal sampling of solvent configurations. We characterize the width of the computed spectrum in terms of the full width at half maximum (FWHM = 2.355 × standard deviation), assuming a Gaussian distribution of the computed values, and comparison to experiment should provide some insight regarding the quality of the underlying MD simulations. Computed FWHMs for Li^{+} and Na^{+} are 0.9–1.1 eV, in good agreement with experimental widths of 1.11–1.24 eV.^{18} Peak width measurements are not available for F^{−}, but for Cl^{−} the measured width is 0.60 eV.^{18} Our computational uncertainties for F^{−}(aq) and Cl^{−}(aq) correspond to a FWHM of 0.8–0.9 eV, slightly larger than what is observed experimentally for Cl^{−}(aq) but in keeping with the trend that the halide anions have narrower photoelectron spectra as compared to the alkali cations. The halides also have narrower photoelectron spectra as compared to neat liquid water, with the latter at 1.45–1.47 eV (experiment^{102,104,105}) versus 1.1–1.2 eV (theory).

One note of caution is in order with regard to spectral widths. Although the favorable comparison between the computed VIEs and the experimental values demonstrates the success of the FragMO/MOM approach, the effect of the lifting of *p*-orbital degeneracy cannot easily be elucidated using this approach. However, the magnitude of the *p*-orbital splitting resulting from an asymmetric distribution of water molecules has been estimated to be rather small,^{18} viz., 0.03 eV for Na^{+}(aq), 0.11 eV for F^{−}(aq), and 0.12 eV for Cl^{−}(aq).

The magnitude of the nonequilibrium correction to the VIE warrants consideration. It might be assumed that this correction, which comes from the boundary conditions, would be rather small given the number of explicit water molecules in our calculations, but in fact this is not the case. The magnitude of the nonequilibrium correction is not available directly from the solvation models (i.e., it is not separable in the total nonequilibrium energy expression) but can be deduced by considering two strategies by which a VIE might be computed.

Employ a nonequilibrium method that uses the optical dielectric constant

*ε*_{∞}for the ionized state and the static dielectric constant*ε*_{solv}for the initial state.Perform two equilibrium solvation free energy calculations, both of which use

*ε*_{solv}, and compute the VIE as the difference between the free energies of the initial and the ionized states.

The free energy of the ionized state computed using the *equilibrium* strategy (2) includes within it the polarization effects resulting from both the slow (nuclear) and fast (electronic) contributions from the solvent. We might express this VIE as

where the notation indicates that *ε*_{solv} is the dielectric constant of merit in both calculations. This differs from *nonequilibrium* strategy (1), which accounts only for the fast component,

The difference between these two calculations provides a measure of the nonequilibrium correction to the VIE,

The average value of ΔVIE_{noneq} from each set of calculations is listed in Table IV. This correction ranges from 0.5 to 1.2 eV, and this value characterizes the error that would be made if only equilibrium solvation models were available. As such, computational strategies for *vertical* ionization energies that are based on *equilibrium* PCMs should not be trusted, although they are sometimes encountered in the literature.

. | ⟨ΔVIE_{noneq}⟩ (eV)
. | ||
---|---|---|---|

. | PEqS . | PCM . | |

Solute . | Hybrid . | Spherical . | SAS . |

Li^{+} | 0.65 | 0.53 | 0.59 |

Na^{+} | 0.65 | 0.53 | 0.59 |

H_{2}O | 0.11 | 0.55 | 0.62 |

e^{−} | 1.17 | 0.54 | 0.61 |

F^{−} | 0.88 | 0.53 | 0.58 |

Cl^{−} | 0.88 | 0.54 | 0.59 |

. | ⟨ΔVIE_{noneq}⟩ (eV)
. | ||
---|---|---|---|

. | PEqS . | PCM . | |

Solute . | Hybrid . | Spherical . | SAS . |

Li^{+} | 0.65 | 0.53 | 0.59 |

Na^{+} | 0.65 | 0.53 | 0.59 |

H_{2}O | 0.11 | 0.55 | 0.62 |

e^{−} | 1.17 | 0.54 | 0.61 |

F^{−} | 0.88 | 0.53 | 0.58 |

Cl^{−} | 0.88 | 0.54 | 0.59 |

### B. VIEs at the liquid/vapor interface

Figure 9 shows the time-dependent VIE, computed using the nonequilibrium PEqS method, for a single MD trajectory of Li^{+} and of Na^{+}, initialized at the liquid/vapor interface. Also plotted is the distance *d*_{GDS} between the ion and the instantaneous GDS.

Within 50 ps, Li^{+} moves away from the interface in favor of a more bulk-like environment, and for the remainder of the simulation its position fluctuates in the range 4.0 Å < *d*_{GDS} < 8.5 Å. Na^{+} departs the interface on a similar time scale and *d*_{GDS} then fluctuates from 6.0 to 9.0 Å for most of the simulation, except when the ion briefly drifts back to the interface around 250–300 ps, before quickly descending again intoa bulk-like solvation environment by 350 ps. Linear regressions of the VIE(*t*) data (shown in blue in Fig. 9) are nearly flat and almost indistinguishable from the average bulk VIE (shown in red), demonstrating the extent to which the VIE is insensitive to the position of the ion relative to the liquid/vapor interface.

Figure 10 presents the same data for representative trajectories of F^{−}(aq) and Cl^{−}(aq) initialized at the interface. The fluoride ion remains within 2–4 Å of the interface for the first 50 ps but then moves away, with *d*_{GDS} fluctuating from 6 to 8 Å. This behavior is similar to the cation data in Fig. 9. In contrast, the migration of Cl^{−} away from the interface is noticeably slower and occurs over 350 ps, with the ion then returning rapidly to within 2 Å of the GDS near the end of the simulation. Linear fits to the instantaneous VIE data are not quite as flat as in the case of the cations, indicating that the VIE exhibits a minor dependence on the location of the anion relative to the interface.

These interfacial simulations were initialized by replacing a water molecule with an ion in a previously equilibrated simulation of neat liquid water, so the early-time dynamics reflect the rearrangement of the solvent to accommodate the ion. We therefore attribute the slope in the linear VIE fits for F^{−} and Cl^{−}, which is not observed for Li^{+} or Na^{+} (where the slopes are essentially zero) as evidence of greater disruption of the water network when the larger and more polarizable halide ions are inserted. Following an equilibration period of roughly 100 ps, however, the interfacial VIEs fluctuate around mean values of 11.63 eV (for F^{−}) and 9.62 eV (for Cl^{−}), which are nearly identical to the average bulk values of 11.66 eV and 9.65 eV. As such, the slight difference in the early-time dynamics of the anions relative to that of the cations seems insignificant and mostly an artifact of the simulation procedure, i.e., the fact that the ion is not equilibrated at the interface at *t* = 0.

Table V compares the average VIEs and average value of *d*_{GDS} for nonequilibrium PEqS calculations of neat liquid water and *e*^{−}(aq). In contrast to the calculations for the interfacial halide anions and alkali cations (e.g., Figs. 9 and 10), where the ion starts at the interface at *t* = 0 but quickly diffuses deeper into the liquid, the simulations leading to Table V are more truly interfacial. The average VIE for liquid water that is reported in Table V is obtained by ionizing the H_{2}O molecule that is closest to the instantaneous GDS at each time step. For *e*^{−}(aq), an electron initialized at the interface remains there long enough to generate a meaningful interfacial trajectory,^{79} and for the snapshots used to compute the average *e*^{−}(aq) VIE in Table V the centroid of the spin density is no farther than 2.5–3.0 Å from the liquid surface. In contrast, halide anions and alkali cations initialized at the interface sample values of *d*_{GDS} in the range 6–8 Å even in the early-time dynamics, as can be seen from the representative trajectories in Figs. 9 and 10.

Solute . | ⟨d_{GDS}⟩ (Å)
. | ⟨VIE⟩ (eV) . |
---|---|---|

H_{2}O | 0.28 ± 0.31 | 11.61 ± 0.52 |

e^{−} | 1.82 ± 0.35 | 3.35 ± 0.46 |

Solute . | ⟨d_{GDS}⟩ (Å)
. | ⟨VIE⟩ (eV) . |
---|---|---|

H_{2}O | 0.28 ± 0.31 | 11.61 ± 0.52 |

e^{−} | 1.82 ± 0.35 | 3.35 ± 0.46 |

Average VIEs for water and for *e*^{−}(aq) reported in Table V can thus be cleanly identified as interfacial VIEs for these species, and the interfacial VIE for liquid water (11.61 ± 0.52 eV) is indistinguishable from the bulk value (11.53 ± 0.51 eV). For *e*^{−}(aq), the interfacial VIE of 3.35 ± 0.46 eV is discernibly lower than the bulk value of 3.75 ± 0.55, albeit not by much. The latter simulations, which are based on the same trajectory data as our previous ones in Ref. 42 but with a slightly better treatment of the electronic structure (including a state-specific PEqS approach rather than a perturbative one) generally support our previous conclusion that the hydrated electron at the air/water interface would be difficult to distinguish from its counterpart in bulk water using liquid microjet photoelectron spectroscopy.

## VI. DISCUSSION

Winter *et al.*^{18} have used liquid microjet photoelectron spectroscopy, in conjunction with a variety of computational strategies, to investigate the VIEs of aqueous halide anions and alkali cations. In the following discussion, we consider the QM/MM *equilibrium* PCM calculations reported in Ref. 18 alongside results from the present work. Errors relative to experiment, for the calculations reported in Ref. 18 and for the present work, are listed in Table VI.

. | Signed error vs. experiment (eV) . | ||||
---|---|---|---|---|---|

. | PEqS (noneq.)^{a}
. | PCM (noneq.)^{a}
. | PCM (eq.) . | . | |

Solute . | Hybrid . | Spherical . | SAS . | Isodensity^{b}
. | QM/MM . |

Li^{+} | 1.01 | 1.45 | 1.19 | 1.83^{c} | 6.43^{c} |

Na^{+} | 0.69 | 1.14 | 0.94 | 0.05^{c} | 4.15^{c} |

H_{2}O | 0.22 | 0.33 | 0.24 | … | … |

e^{−} | −0.0 | −0.5 | −0.5 | … | … |

F^{−} | 0.08 | −0.21 | −0.06 | −3.78^{d} | −0.39^{d} |

Cl^{−} | 0.05 | −0.24 | −0.19 | −2.75^{d} | −0.93^{d} |

. | Signed error vs. experiment (eV) . | ||||
---|---|---|---|---|---|

. | PEqS (noneq.)^{a}
. | PCM (noneq.)^{a}
. | PCM (eq.) . | . | |

Solute . | Hybrid . | Spherical . | SAS . | Isodensity^{b}
. | QM/MM . |

Li^{+} | 1.01 | 1.45 | 1.19 | 1.83^{c} | 6.43^{c} |

Na^{+} | 0.69 | 1.14 | 0.94 | 0.05^{c} | 4.15^{c} |

H_{2}O | 0.22 | 0.33 | 0.24 | … | … |

e^{−} | −0.0 | −0.5 | −0.5 | … | … |

F^{−} | 0.08 | −0.21 | −0.06 | −3.78^{d} | −0.39^{d} |

Cl^{−} | 0.05 | −0.24 | −0.19 | −2.75^{d} | −0.93^{d} |

Due to complications arising from water ionization in the presence of explicit solvent molecules, as discussed in Sec. IV C, the equilibrium PCM calculations in Ref. 18 include only the bare ion in the QM region. The equilibrium nature of the PCM used in that work affords an *adiabatic* ionization energy (AIE) because all solvent degrees of freedom are (implicitly) relaxed following ionization. (For a molecular solute, or if explicit solvent molecules are included in the QM region, the use of an equilibrium PCM without geometry optimization in the ionized state affords something in between a VIE and an AIE because the implicit solvent degrees of freedom are relaxed but the explicit nuclear degrees of freedom are not.) QM/MM calculations reported in Ref. 18 also include only the solute in the QM region, surrounded by classical point charges representing water molecules. Unlike the equilibrium PCM calculations, ionization energies computed using this approach are indeed VIEs, albeit ones that lack any electronic polarization contribution from the solvent because the point charges cannot be polarized upon ionization of the solute.

For the alkali cations, VIEs computed at the QM/MM level err by 6.4 eV (Li^{+}) and 4.2 eV (Na^{+}) and are clearly unacceptable. AIEs computed with the equilibrium PCM approach are larger than experimental VIEs, by 1.8 eV for Li^{+}(aq) but only by 0.05 eV for Na^{+}(aq). This large discrepancy in accuracy is puzzling and is likely fortuitous, but in any case the juxtaposition of calculated AIEs with experimental VIEs is not reasonable, especially for the halides where the hydration structure of X^{−} differs considerably from that of neutral X. Overall, the enormity of the errors for this approach and for the QM/MM calculations suggests that a proper description of “specific” solvent effects, by means of explicit QM water molecules, is essential, even when using a PCM. (This fact is well known, e.g., in p*K*_{a} calculations.^{111}) With respect to the nonequilibrium PCM results with explicit solvent, for which a direct comparison to experiment is appropriate, errors for both Li^{+}(aq) and Na^{+}(aq) are 0.94–1.45 eV, which we consider to be surprisingly large given the rather simple electronic structure of these solutes.

MP2/PEqS calculations including ≈30 explicit QM water molecules represent our best attempt for these systems, yet these calculations still overestimate the cation VIEs by 1.0 eV (Li^{+}) and 0.7 eV (Na^{+}). (That the error for Na^{+} is more comparable to that for Li^{+}, as compared to the calculations reported in Ref. 18, suggests that the accuracy of the equilibrium PCM result for Na^{+} is indeed fortuitous.) The reasons behind this remaining error remain a topic for further study; a more thorough examination of cavity construction is probably warranted at the very least.

As compared to the alkali cations, where all methods considered here and in Ref. 18 overestimate the experimental VIE, both positive and negative errors are observed for the halide ions, perhaps simply because the errors are closer to zero. (Exceptions are the AIEs computed with an equilibrium PCM, which are much too small.) The QM/MM results are much more accurate than they were for the cations, but this seems fortuitous given that an anionic QM solute likely suffers more from overpolarization by the point charges than does a cation solute due to the more diffuse nature of the wave function. Notably, the equilibrium PCM results for the anions are far less accurate than those for the cations, but the relative accuracy for cations represents a form of error cancellation as suggested in Ref. 18. Namely, for the cations, only minor reorientation of the solvent molecules is required to accommodate the resulting divalent ion, due to the pre-existing favorable alignment of the solvent dipoles in the monovalent state, but ionization of the anions results in a charge-neutral solute and thus significant reorientation of the solvent. As a result, one expects the AIE to be much smaller than the VIE for the anions, but more similar to the VIE for the cations.

Considering *e*^{−}(aq), we note that the experimental VIE of 3.7 ± 0.1 eV^{15} that is provided in Table III represents an upward revision of many previously reported VIEs in the range 3.3–3.4 eV.^{14,112–116} The newer value has been called the “genuine” binding energy of *e*^{−}(aq),^{15} as it includes corrections for scattering of the ejected electron that lead to wavelength dependence in the photoelectron spectrum.^{14,15} MP2/PEqS calculations also afford a VIE of 3.7 eV, which is rather remarkable given the complexity of this species though not out of line with the accuracy that we obtain for the halide anions. (Our calculation also represents an upward revision of the MP2/6-31++G* value that we reported previously, based on a perturbative version of the nonequilibrium PEqS method.^{42}) Assuming that the new experimental VIE withstands further scrutiny, then the present MP2/PEqS calculation would seem to affirm the simulation procedures^{71,79} used to generate the MD snapshots for this species. It also adds to the ongoing debate regarding the detailed structure of the aqueous electron,^{70,71,117–124} and in particular the question of whether this species is cavity-forming or not; the MD snapshots used here correspond to a cavity-forming electron, as can be seen in Fig. 1. Finally, we note that amongst the solutes examined here *e*^{−}(aq) affords the largest discrepancy between PEqS and PCM values of the VIE. This may represent a failure of IEF-PCM to adequately describe volume polarization by this highly delocalized non-classical solute.^{110}

Finally we consider the 1*b*_{1} state of liquid water. Errors in the VIE, obtained with nonequilibrium methods, are 0.2–0.3 eV. Crucial to the success of this method is the FragMO/MOM SCF technique, in order to ionize a central H_{2}O molecule not too near the QM/continuum boundary. This approach will require some care if the full photoelectron spectrum is desired, to select the appropriate orbitals for ionization, and may be hopeless in cases where the MO picture of photoionization completely breaks down.^{125} Unfortunately, alternative electronic structure techniques that afford ionization energies directly,^{126–130} and can deal with “shake-up” ionization events,^{125,126,130} are considerably more expensive. Alternatively, a simple technique to correct the Kohn-Sham eigenvalue spectrum has recently been shown to afford valence photoelectron spectra that compare well with experiment^{131} and is no more expensive than standard density-functional calculations. In our estimation, however, this “potential adjustors” technique^{131} is unlikely to work across the entire (core + valence) photoelectron spectrum.

## VII. CONCLUSIONS

We have presented a detailed description of the theory and implementation of the state-specific nonequilibrium PEqS method and its application to compute aqueous-phase VIEs. In contrast to PCMs, which are the *de facto* implicit solvation models in electronic structure calculations, PEqS calculations require discretization of three-dimensional space and not simply a two-dimensional cavity surface. This makes PEqS considerably more expensive than PCM calculations despite the efficient multigrid approach described here. Computational expense notwithstanding, the PEqS approach has the advantage that it treats volume polarization (charge leakage outside of the QM region) *exactly*, up to discretization errors. That said, among the solutes considered here this seems to matter only for *e*^{−}(aq). For more “classical” solutes, nonequilibrium PCM calculations with the solvent accessible surface construction^{73} are within 0.2 eV of PEqS results, thus validating the more affordable PCM approach in bulk solution.

A more important advantage of the PEqS approach is that it is naturally applicable to arbitrary (and therefore anisotropic) dielectric environments, defined by a dielectric function *ε*(**r**). This function could be defined based on the electron density,^{40,62–67} but here we adopt an approach analogous to PCM calculations and define a surface to delineate the boundary between the atomistic QM region and its continuum environment. The usual PCM prescription using atom-centered vdW spheres, however, proves to be problematic when explicit water molecules are included in the QM region, leading to unphysical high-dielectric regions between these explicit solvent molecules. Oddly, this problem is not often discussed in the quantum chemistry literature although a similar problem in biomolecular simulation has been widely discussed,^{132–136} where in the context of Poisson-Boltzmann electrostatics calculations the vdW cavity construction may leave high-dielectric regions in the hydrophobic interior of a protein. In this work, we introduced a “hybrid” cavity model that avoids this problem.

Here, we used the PEqS approach, in conjunction with QM regions containing ≈30 explicit water molecules, to compute VIEs for neat liquid water as well as F^{−}(aq), Cl^{−}(aq), Na^{+}(aq), Li^{+}(aq), and *e*^{−}(aq), both in bulk liquid water and at the air/water interface. Ionization energies for most of these systems lie below (or are similar to) that of liquid water itself, and a naïve calculation of the lowest-energy ionized state thus results in ionization of H_{2}O rather than the solute of interest. We circumvent this problem by means of a fragment-based initial guess combined with the maximum overlap method.^{106}

We find that nonequilibrium corrections to VIEs, which are missing from continuum models based only on the static dielectric constant, amount to 0.5–1.2 eV for each system investigated in this work. VIEs computed at the MP2/PEqS (noneq.) level for liquid water, F^{−}(aq), Cl^{−}(aq), and *e*^{−}(aq) agree with experimental results to within 0.2 eV, with slightly larger errors when nonequilibrium PCMs are used as a substitute for PEqS. For reasons that remain unclear, however, errors for alkali cations are larger, e.g., 1.0 eV for Li^{+}(aq) at the MP2/PEqS (noneq.) level. Consistent with our previous work on *e*^{−}(aq),^{42} there is very little difference between VIEs computed at the air/water interface versus those in bulk water, for any of the solutes considered here. For liquid water, the same conclusion has recently been reported based on *G*_{0}*W*_{0} calculations.^{137}

All of the QM calculations in this work were performed at the MP2 level, but the PEqS method works equally well in the context of density functional theory, and also for other correlated wave functions, if the Hartree-Fock density is used in Poisson’s equation. A nonequilibrium treatment of vertical excitation energies within the PEqS framework could be accomplished by adapting PCM algorithms described previously in the context of time-dependent density functional theory^{48,49} and the algebraic diagrammatic construction.^{50} Finally, it should be possible to adapt this methodology to develop nonequilibrium versions of density-dependent dielectric solvation models in which *ε*(**r**) is a functional of *ρ*(**r**).^{40,62–64,138} This would eliminate some of the arbitrariness in construction of the QM/continuum boundary. We hope to address this, along with the sensitivity of PEqS results to the Gaussian smearing of the nuclear charges, in future work.

## ACKNOWLEDGMENTS

This work was supported by National Science Foundation Grant Nos. CHE-1300603 and CHE-1665322. Calculations were performed at the Ohio Supercomputer Center under Project No. PAA-0003.^{139} J.M.H. is a fellow of the Alexander von Humboldt Foundation and serves on the Board of Directors of Q-Chem, Inc.

## REFERENCES

In 2016, the experimental value for the VIE of F^{−}(aq) was revised upward from 8.7 eV^{17} to 11.6 eV,^{11} with much earlier theoretical calculations^{18} serving to validate the revised value. The F^{−}(aq) VIE is difficult to discern as it is embedded in the 1*b*_{1} band of the photoelectron spectrum of liquid water.^{11}

^{−}

Reference 72 considers a four-coordinate, tetrahedral $(H2O)4\u2212$ model for *e*^{−}(aq), using nonequilibrium IEF-PCM to represent bulk water. The distance between each of the four O–H moieties and the center of the tetrahedron was selected partly based on comparing the calculated VIE, which varies dramatically as a function of this distance coordinate, to an experimental value taken to be 3.42–3.45 eV. According to the data provided in Ref. 72, a much smaller tetrahedron would be required in order to accommodate the most recent measurement, 3.7 eV.^{15} However, a smaller tetrahedron affords results that are inconsistent with other experimental data for *e*^{−}(aq), according to the calculations in Ref. 72. In view of the sizable disparity (0.5 eV) between PEqS and IEF-PCM values of the hydrated electron’s VIE, as reported in the present work, this may point to the inadequacy of PCMs for this species.

_{a}calculations with continuum solvents, alternative protocols to thermodynamic cycles