Quantum chemistry in arbitrary dielectric environments: Theory and implementation of nonequilibrium Poisson boundary conditions and application to compute vertical ionization energies at the air/water interface

Coons, Marc P.; Herbert, John M.

doi:10.1063/1.5023916

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,¹⁶ 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,²² 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 ε₁,ε₂, … in different spatial regions, was introduced long ago by Sakurai et al.^37,38 and used in semi-empirical electronic structure calculations.³⁸ 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⁴⁸), 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.¹¹ Perhaps unsurprisingly, VIEs for atomic ions computed using nothing but a PCM representation of the solvent afford extremely poor agreement with experiment;¹⁸ 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.⁵¹ 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,

{\hat{\nabla}}^{2} φ_{s o l} (r) = - 4 π ρ_{s o l} (r),

(2.1)

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,

φ_{s o l} (r) = φ_{n u c} (r) + φ_{e l e c} (r),

(2.2a)

ρ_{s o l} (r) = ρ_{n u c} (r) + ρ_{e l e c} (r) .

(2.2b)

The solute’s internal energy in a vacuum is

E_{i n t} = \frac{1}{2} \int d r φ_{s o l} (r) ρ_{s o l} (r)

(2.3)

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,

\hat{\nabla} \cdot [ε (r) \hat{\nabla} φ_{t o t} (r)] = - 4 π ρ_{s o l} (r),

(2.4)

where

φ_{t o t} (r) = φ_{s o l} (r) + φ_{p o l} (r)

(2.5)

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

φ_{e l e c} (r) = - \sum_{μ ν} P_{μ ν} φ_{μ ν} (r),

(2.6)

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

φ_{μ ν} (r) = \int d r^{'} \frac{g_{μ} (r) g_{ν} (r^{'})}{| r - r^{'} |} .

(2.7)

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

φ_{n u c} (r) = \sum_{α}^{a t o m s} \frac{Z_{α}}{| r - R_{α} |} erf (\frac{| r - R_{α} |}{\sqrt{2} σ}),

(2.8)

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

ρ_{s o l} (r) = - \frac{1}{4 π} {\hat{\nabla}}^{2} φ_{s o l} (r) .

(2.9)

In this work, ${\hat{\nabla}}^{2} φ_{s o l}$ (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,⁴² 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,

{\hat{\nabla}}^{2} φ_{t o t} (r) = - 4 π ρ_{t o t} (r),

(2.10)

where the total charge density is

ρ_{t o t} (r) = ρ_{s o l} (r) + ρ_{p o l} (r) .

(2.11)

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

ρ_{p o l} (r) = [\frac{1 - ε (r)}{ε (r)}] ρ_{s o l} (r) + ρ_{i t e r} (r) .

(2.12)

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 kth iteration can be expressed as^40,41

ρ_{i t e r}^{(k)} (r) = \frac{1}{4 π} [\hat{\nabla} \ln ε (r)] \cdot [\hat{\nabla} φ_{t o t}^{(k)} (r)] .

(2.13)

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 $ρ_{i t e r}^{(k)}$ (r) is computed using Eq. (2.13), and Eq. (2.12) is then used to generate $ρ_{p o l}^{(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 $φ_{t o t}^{(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 $ρ_{i t e r}^{(k + 1)}$ (r), we instead use a damping procedure to stabilize the update between iterations k and k + 1,

\begin{matrix} ρ_{i t e r}^{(k + 1)} (r) & = & \frac{η}{4 π} [\hat{\nabla} \ln ε (r)] \cdot [\hat{\nabla} φ_{t o t}^{(k + 1)} (r)] \\ + (1 - η) ρ_{i t e r}^{(k)} (r) . \end{matrix}

(2.14)

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

Δ ρ_{i t e r} = ∥) ρ_{i t e r}^{(k + 1)} (r) - ρ_{i t e r}^{(k)} (r) ∥)

(2.15)

falls below a threshold, τ_PEqS.

Algorithm 1.

Equilibrium PEqS method.

1: Initialize Δh = 0.

2: for j = 1, 2, … do until SCF error < τ_SCF

3: Diagonalize

F^{(j)} = F_{0}^{(j)} + Δ 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

ρ_{i t e r}^{(0)}

(r) using

φ_{t o t}^{(0)}

(r)

13: for k = 0, 1, … do until Δρ_iter < τ_PEqS

14: Compute

φ_{t o t}^{(k + 1)}

(r)

15: Update

ρ_{i t e r}^{(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₀ with a correction Δh to its one-electron part. This correction has matrix elements

Δ h_{μ ν} = - \int d r φ_{μ ν} (r) ρ_{p o l} (r) .

(2.16)

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

E = E_{i n t} + \frac{1}{2} \int d r φ_{p o l} (r) ρ_{s o l} (r),

(2.17)

which consists of the solute’s internal energy E_int from the electronic structure calculation, plus the electrostatic contribution

G_{e l s t} = \frac{1}{2} \int d r φ_{s o l} (r) ρ_{p o l} (r)

(2.18)

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,⁵⁹ 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

Ĥ_{i}^{S S} | Ψ_{i} ⟩) = E_{i}^{S S} | Ψ_{i} ⟩)

(2.19)

with a state-specific Hamiltonian of the form

Ĥ_{i}^{S S} = Ĥ_{i}^{v a c} + {\hat{V}}_{0}^{s l o w} + {\hat{V}}_{i}^{f a s t} .

(2.20)

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 $Ĥ_{i}^{v a c}$ is the molecular Hamiltonian that affords the solute’s vacuum internal energy E_int,i for state i, and the operators ${\hat{V}}_{i}^{s l o w}$ and ${\hat{V}}_{i}^{f a s t}$ generate the indicated components of the solvent polarization response, i.e., $φ_{p o l, i}^{s l o w}$ and $φ_{p o l, i}^{f a s t}$ ⁠, where

φ_{p o l, i}^{s l o w / f a s t} (r) = \int d r^{'} \frac{ρ_{p o l, i}^{s l o w / f a s t} (r^{'})}{| r - r^{'} |} .

(2.21)

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, $ρ_{p o l, 0}^{s l o w}$ (r), which affords $φ_{p o l, 0}^{s l o w}$ (r) according to Eq. (2.21), is computed according to^48,49

ρ_{p o l, 0}^{s l o w} (r) = (\frac{χ_{s l o w}}{χ}) ρ_{p o l, 0} (r),

(2.22)

where χ = χ_slow + χ_fast is the static susceptibility, partitioned into slow and fast components,

χ_{s l o w} = \frac{ε_{s o l v} - ε_{\infty}}{4 π},

(2.23a)

χ_{f a s t} = \frac{ε_{\infty} - 1}{4 π} .

(2.23b)

Here, ε_solv is the static dielectric constant of the solvent and ε_∞ = n² 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, $φ_{p o l, 1}^{f a s t}$ (r) and $ρ_{p o l, 1}^{f a s t}$ (r), within the Marcus partition,⁴⁹ the Poisson equation is modified such that the total source-charge density is the ionized solute’s charge density, $ρ_{s o l, 1} (r) + ρ_{p o l, 0}^{s l o w}$ (r), and the dielectric function is the optical one. This modified form of Eq. (2.4) is

\hat{\nabla} \cdot [ε_{\infty} (r) \hat{\nabla} φ_{t o t, 1}^{f a s t} (r)] = - 4 π [ρ_{s o l, 1} (r) + ρ_{p o l, 0}^{s l o w} (r)] .

(2.24)

Here, $φ_{t o t, 1}^{f a s t}$ (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

{\hat{\nabla}}^{2} φ_{t o t, 1}^{f a s t} (r) = - 4 π ρ_{t o t, 1}^{f a s t} (r),

(2.25)

where $ρ_{t o t, 1}^{f a s t}$ (r) and $φ_{t o t, 1}^{f a s t}$ (r) are to be computed self-consistently. The total nonequilibrium source charge density is

ρ_{t o t, 1}^{f a s t} (r) = ρ_{s o l, 1} (r) + ρ_{p o l, 0}^{s l o w} (r) + ρ_{p o l, 1}^{f a s t} (r) .

(2.26)

For the Marcus partitioning scheme,⁴⁹ $ρ_{p o l, 1}^{f a s t}$ (r) takes the form

ρ_{p o l, 1}^{f a s t} (r) = (\frac{1 - ε_{\infty} (r)}{ε_{\infty} (r)}) [ρ_{s o l, 1} (r) + ρ_{p o l, 0}^{s l o w} (r)] + ρ_{i t e r, 1} (r),

(2.27)

where ρ_iter,1(r) is computed iteratively according to

ρ_{i t e r, 1}^{(k)} (r) = \frac{1}{4 π} [\hat{\nabla} \ln ε_{\infty} (r)] \cdot [\hat{\nabla} φ_{t o t, 1}^{f a s t, (k)} (r)] .

(2.28)

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

Finally, the nonequilibrium free energy is^48,49

\begin{matrix} G_{e l s t, 1} & = & W_{0,1} + \int d r φ_{s o l, 1} (r) ρ_{p o l, 0}^{s l o w} (r) \\ + \frac{1}{2} \int d r φ_{s o l, 1} (r) ρ_{p o l, 1}^{f a s t} (r) \\ - \frac{1}{2} \int d r φ_{s o l, 0} (r) ρ_{p o l, 0}^{s l o w} (r) . \end{matrix}

(2.29)

The term

W_{0,1} = \frac{1}{2} \int d r φ_{p o l, 0}^{s l o w} (r) [ρ_{p o l, 1}^{f a s t} (r) - ρ_{p o l, 0}^{f a s t} (r)]

(2.30)

arises within the Marcus partitioning scheme⁴⁹ 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

\begin{matrix} VIE & = & Δ E_{1,0} + W_{0,1} + \frac{1}{2} \int d r φ_{s o l, 1} (r) ρ_{p o l, 1}^{f a s t} (r) \\ + \int d r [φ_{s o l, 1} (r) - φ_{s o l, 0} (r)] ρ_{p o l, 0}^{s l o w} (r) \\ - \frac{1}{2} \int d r φ_{s o l, 0} (r) ρ_{p o l, 0}^{f a s t} (r), \end{matrix}

(2.31)

where

Δ E_{1,0} = E_{i n t, 1} - E_{i n t, 0}

(2.32)

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₁ whose elements are

Δ h_{μ ν, 1} = - \int d r φ_{μ ν} (r) [ρ_{p o l, 0}^{s l o w} (r) + ρ_{p o l, 1}^{f a s t} (r)] .

(2.33)

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

Algorithm 2.

Nonequilibrium PEqS method.

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

2: Reference data: E₀, G_elst,0,

φ_{p o l, 0}^{s l o w}

(r), and

ρ_{p o l, 0}^{s l o w}

(r)

3: Initialize Δh₁ = 0

4: for j = 1, 2, … do until SCF error < τ_SCF

5: Diagonalize

F^{(j)} = F_{0}^{(j)} + Δ h_{1}^{(j)}

to obtain P^(j)

6: Compute ρ_sol,1(r) and φ_sol,1(r)

7: if j = 1 then

8:

ρ_{t o t, 1}^{f a s t} (r) = ρ_{s o l, 1} (r) + ρ_{p o l, 0}^{s l o w}

(r)

9:

φ_{t o t, 1}^{f a s t} (r) = φ_{s o l, 1} (r) + φ_{p o l, 0}^{s l o w}

(r)

10: else

11:

ρ_{t o t, 1}^{f a s t} (r) = ρ_{s o l, 1} (r) + ρ_{p o l, 0}^{s l o w} (r) + ρ_{p o l, 1}^{f a s t} (r)

12:

φ_{t o t, 1}^{f a s t} (r) = φ_{s o l, 1} (r) + φ_{p o l, 0}^{s l o w} (r) + φ_{p o l, 1}^{f a s t}

(r)

13: end if

14: Initialize

ρ_{i t e r, 1}^{(0)}

(r) using

φ_{t o t, 1}^{f a s t, (0)}

(r)

15: for k = 0, 1, … do until Δρ_iter,1 < τ_PEqS

16: Compute

φ_{t o t, 1}^{f a s t, (k + 1)}

(r)

17: Update

ρ_{i t e r, 1}^{(k + 1)}

(r),

ρ_{p o l, 1}^{f a s t}

(r), and

ρ_{t o t, 1}^{f a s t}

(r)

18: end for

19: Update

Δ h_{1}^{(j)}

20: Compute E₁ = E_int,1 + G_elst,1

21: end for

22: Compute VIE = E₁ − E₀

C. Dielectric function

This section describes the construction of the dielectric function ε(r) that appears in Eqs. (2.4) and (2.24), for both bulk solvation and the liquid/vapor interface.

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.⁶¹

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

ε (r) = ε_{v a c} + (ε_{s o l v} - ε_{v a c}) \prod_{α}^{a t o m s} s_{α} (d_{α}, Δ; | r - R_{α} |) .

(2.34)

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 ${(ε_{v a c} r)}^{- 1}$ rather than r⁻¹. All numerical calculations presented here use ε_vac = 1.

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

s_{α} (d_{α}, Δ; | r - R_{α} |) = \frac{1}{2} [1 + e r f (\frac{| r - R_{α} | - d_{α}}{Δ})],

(2.35)

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,⁶⁸ except that for hydrogen we use the updated value of 1.1 Å.⁶⁹ 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),⁷⁰ represented in the following example as a ${(H_{2} O)}_{28}^{-}$ cluster model extracted from a condensed-phase simulation.⁷¹ 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),⁷² but this is potentially problematic because the vdW cavity that is conventionally used in PCM calculations places high-dielectric regions in between water molecules.

FIG. 1.

View large Download slide

Singly occupied molecular orbital of a ${(H_{2} O)}_{28}^{-}$ cluster model of e⁻(aq). The opaque and translucent isosurfaces encapsulate 50% and 95% of the probability density, respectively. Positive values of the orbital are shown in blue, and the very small negative regions are shown in green. The latter are only visible in the 95% isoprobability contour.

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^{-} {(H_{2} O)}_{31}$ cluster that is plotted in Fig. 3. High-dielectric regions can once again be found inside of the atomistic QM region.

FIG. 2.

View large Download slide

Two-dimensional slices through the function ε(r), for the ${(H_{2} O)}_{28}^{-}$ cluster that is shown in Fig. 1, which was extracted from a simulation of e⁻(aq). The dielectric function is constructed using either (a) the vdW solute cavity, Eq. (2.34) with parameters d_α set to scaled Bondi radii; (b) a “modified” SAS construction, created by setting d_α = r_vdW,α + 0.7 Å, which differs from the usual SAS choice, r_probe = 1.4 Å; or (c) a “hybrid” cavity, which is described in the context of Eq. (2.36). Each panel plots ε(r) in the xz plane that contains the cluster center of mass (c.o.m.). The dielectric function transitions smoothly from ε = 1 inside the cavity to ε = 78.39 outside.

FIG. 3.

View large Download slide

Two-dimensional slices through the function ε(r) for a cluster $F^{-} {(H_{2} O)}_{31}$ extracted from a simulation of F⁻(aq). The vdW cavity is used to construct the dielectric function, so this plot is analogous to that shown in Fig. 2(a) for the hydrated electron.

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).⁷³ 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.⁷⁴ 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 ${(H_{2} O)}_{n}^{-}$ clusters. The shape of each solute configuration (water cluster) is approximated as an ellipsoid centered at the cluster c.o.m. (x₀, y₀, z₀), whose surface is defined by the equation S(x, y, z) ≡ 1, where

S (x, y, z) = \frac{{(x - x_{0})}^{2}}{a^{2}} + \frac{{(y - y_{0})}^{2}}{b^{2}} + \frac{{(z - z_{0})}^{2}}{c^{2}} .

(2.36)

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 ${(H_{2} O)}_{28}^{-}$ 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.⁴² 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

ρ (z) = \frac{1}{2} ρ_{solv} {1 \pm \tanh [α (z - z_{G D S})]} .

(2.37)

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

ε (z) = \frac{1}{2} ε_{s o l v} {1 \pm \tanh [α (z - z_{G D S})]} .

(2.38)

Figure 4 illustrates the dielectric function ε(r) for a ${(H_{2} O)}_{24}^{-}$ cluster extracted from an interfacial configuration of e⁻(aq),⁷⁹ with the c.o.m. placed at the origin.

TABLE I.

Parameters for Eq. (2.37), obtained by fitting ensemble-averaged solvent density profiles from MD simulations.

Solute	ρ_solv (g/cm³)	α (Å⁻¹)	z_GDS (Å)
Li⁺	0.976	0.652	−9.205
Na⁺	0.979	0.604	−9.154
H₂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³)	α (Å⁻¹)	z_GDS (Å)
Li⁺	0.976	0.652	−9.205
Na⁺	0.979	0.604	−9.154
H₂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

FIG. 4.

View large Download slide

Two-dimensional slice through ε(r) for a ${(H_{2} O)}_{24}^{-}$ cluster representing e⁻(aq) at the liquid/vapor interface. A hybrid cavity is first constructed, as described in the discussion surrounding Eq. (2.36), and then Eq. (2.38) is used to interpolate the dielectric from ε_solv → ε_vac across the GDS, which is indicated by the dashed black line (z_GDS = −1.508 Å). Out of 24 explicit H₂O molecules inside the cavity, only 2–3 lie above the GDS.

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 (N_{g r i d}^{3})$ in floating-point operations and $O (N_{g r i d}^{2})$ 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⁶ 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

For the discussion that follows, let us rewrite Eqs. (2.10) and (2.25) in a generic form

{\hat{\nabla}}^{2} φ (r) = ρ (r),

(3.1)

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

\frac{\partial^{2} φ}{\partial x^{2}} + \frac{\partial^{2} φ}{\partial y^{2}} + \frac{\partial^{2} φ}{\partial z^{2}} = ρ (x, y, z),

(3.2)

subject to the Dirichlet boundary condition

φ (x, y, z) = 0, \forall (x, y, z) \in δ Ω

(3.3)

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_xN_yN_z), the domain Ω is defined as

\begin{matrix} Ω = [) (- L_{x} / 2 < x < L_{x} / 2), (- L_{y} / 2 < y < L_{y} / 2), \\ (- L_{z} / 2 < z < L_{z} / 2)]) . \end{matrix}

(3.4)

The surface boundary is defined by the collection of rectangular planes

δ Ω = {[\pm L_{x}, y, z], [x, \pm L_{y}, z], [x, y, \pm L_{z}]} .

(3.5)

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

φ_{i, j, k} ≔ φ (x_{i}, y_{j}, z_{k}),

(3.6)

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 nth-order approximation to the mth-order derivative, whose finite-difference approximation exhibits error of $O (h^{2 n})$ ⁠, is

\frac{\partial^{m} φ_{i, j, k}}{\partial x^{m}} = \sum_{p = - n}^{n} c_{m, p} (\frac{φ_{i + p, j, k}}{h^{m}}),

(3.7)

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.

TABLE II.

Central finite-difference coefficients c_m,±p for the discretized first (m = 1) and second (m = 2) derivatives in Eq. (3.7). These coefficients afford an eighth-order approximation scheme whose accuracy is $O (h^{8})$ ⁠.

	c_m,0	c_m,±1	c_m,±2	c_m,±3	c_m,±4
m = 1	0	±4/5	∓1/5	±4/105	∓1/280
m = 2	−205/72	8/5	−1/5	8/315	−1/560

	c_m,0	c_m,±1	c_m,±2	c_m,±3	c_m,±4
m = 1	0	±4/5	∓1/5	±4/105	∓1/280
m = 2	−205/72	8/5	−1/5	8/315	−1/560

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,

L^{h} φ^{h} = ρ^{h} .

(3.8)

The N_grid × N_grid matrix $L^{h}$ contains the discretized Laplacian operator, and the vectors φ^h and ρ^h contain the discretized values $φ_{i, j, k,}^{h}$ and $- 4 π ρ_{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 $L^{h}$ explicitly; only its action on the vector φ^h is required.

Via Fourier analysis of the discretization errors

v_{i, j, k}^{h} = φ_{i, j, k}^{e x a c t} - φ_{i, j, k}^{h},

(3.9)

it has been shown that the spectrum of errors contains wavelengths λ that are comparable to, or larger than, the grid resolution, h.⁸³ The CG routine efficiently eliminates discretization errors with λ ≃ h but struggles with error components for which λ > h.⁸³ 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₂O molecule placed at the center of a (15 Å)³ cubic Cartesian grid with a resolution h = 0.074 Å. Figure 5 shows the Euclidean norm of the residual error in the electrostatic potential,

r^{h} = ρ^{h} - L^{h} φ^{h},

(3.10)

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⁻⁵ 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.

FIG. 5.

View large Download slide

Comparison of the performance of several iterative schemes for solving Eq. (3.8) to an accuracy τ_PEqS = 10⁻⁵ a.u. (indicated by the horizontal black line), for a single water molecule whose c.o.m. resides at the center of a (15 Å)³ Cartesian grid with h = 0.074 Å. The slow decay of the residual error for the standard CG routine is indicative of λ > h, and over 250 iterations are required for convergence. The multigrid methods achieve convergence more rapidly, with the W-cycle implementation performing best. Both multigrid algorithms exhibit an exponential decrease of r^h with respect to the iteration number, in stark contrast to results from the standard CG method.

Since the charge density ρ(r) computed by the electronic structure calculation is known essentially to arbitrary accuracy, at least compared to the ∼10⁻⁵ 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 $ρ^{h} = L^{h} φ^{e x a c t}$ ⁠. 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,⁸³

\begin{align} r^{h} & = L^{h} φ^{e x a c t} - L^{h} φ^{h} \\ = L^{h} (φ^{e x a c t} - φ^{h}) \\ = L^{h} v^{h} . \end{align}

(3.11)

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⁸³ 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 = 2h, 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 $I_{h}^{H}$ ⁠,

r^{H} = I_{h}^{H} r^{h} .

(3.12)

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

\begin{matrix} r_{I, J, K}^{H} & = & \frac{1}{8} r_{i, j, k}^{h} + \frac{1}{16} (r_{i + 1, j, k}^{h} + r_{i - 1, j, k}^{h} + r_{i, j + 1, k}^{h} + r_{i, j - 1, k}^{h} + r_{i, j, k + 1}^{h} + r_{i, j, k - 1}^{h}) + \frac{1}{32} (r_{i + 1, j + 1, k}^{h} + r_{i + 1, j - 1, k}^{h} + r_{i - 1, j + 1, k}^{h} + r_{i - 1, j - 1, k}^{h}) \\ + \frac{1}{32} (r_{i + 1, j, k + 1}^{h} + r_{i + 1, j, k - 1}^{h} + r_{i - 1, j, k + 1}^{h} + r_{i - 1, j, k - 1}^{h}) + \frac{1}{32} (r_{i, j + 1, k + 1}^{h} + r_{i, j + 1, k - 1}^{h} + r_{i, j - 1, k + 1}^{h} + r_{i, j - 1, k - 1}^{h}) \\ + \frac{1}{64} (r_{i + 1, j + 1, k + 1}^{h} + r_{i + 1, j + 1, k - 1}^{h} + r_{i + 1, j - 1, k + 1}^{h} + r_{i + 1, j - 1, k - 1}^{h}) + \frac{1}{64} (r_{i - 1, j + 1, k + 1}^{h} + r_{i - 1, j + 1, k - 1}^{h} + r_{i - 1, j - 1, k + 1}^{h} + r_{i - 1, j - 1, k - 1}^{h}) . \end{matrix}

(3.13)

The notation (I, J, K) in $r_{I, J, K}^{H}$ 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 = 2h_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 $r_{I, J, K}^{H}$ ⁠.

FIG. 6.

View large Download slide

Illustration of a two-level V-cycle multigrid algorithm applied to solve Eq. (3.8). The input is a source charge density ρ^h and the output is φ^h, both of which are discretized on a target grid whose resolution is h. Steps 1 and 2 show the formation of the residual error r^h on the target grid and subsequently its restriction to a coarser grid whose resolution is H = 2h. On the coarser grid, r^H is used in a CG routine to compute a relaxed residual vector r^H in step 3. The relaxed residual is then interpolated back to the target grid to form v^h in step 4. In step 5, v^h is used to correct the solution on the target grid, and this process is repeated (starting from step 2) until convergence.

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

L^{H} v^{H} = r^{H} .

(3.14)

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 $I_{H}^{h}$ ⁠,

v^{h} = I_{H}^{h} v^{H} .

(3.15)

The action of $I_{H}^{h}$ is expressed by the following set of equations:

v_{i, j, k}^{h} = v_{I, J, K}^{H},

(3.16a)

v_{i + 1, j, k}^{h} = \frac{1}{2} (v_{I, J, K}^{H} + v_{I + 1, J, K}^{H}),

(3.16b)

v_{i, j + 1, k}^{h} = \frac{1}{2} (v_{I, J, K}^{H} + v_{I, J + 1, K}^{H}),

(3.16c)

v_{i, j, k + 1}^{h} = \frac{1}{2} (v_{I, J, K}^{H} + v_{I, J, K + 1}^{H}),

(3.16d)

v_{i + 1, j + 1, k}^{h} = \frac{1}{4} (v_{I, J, K}^{H} + v_{I + 1, J, K}^{H} + v_{I, J + 1, K}^{H} + v_{I + 1, J + 1, K}^{H}),

(3.16e)

v_{i + 1, j, k + 1}^{h} = \frac{1}{4} (v_{I, J, K}^{H} + v_{I + 1, J, K}^{H} + v_{I, J, K + 1}^{H} + v_{I + 1, J, K + 1}^{H}),

(3.16f)

v_{i, j + 1, k + 1}^{h} = \frac{1}{4} (v_{I, J, K}^{H} + v_{I, J + 1, K}^{H} + v_{I, J, K + 1}^{H} + v_{I, J + 1, K + 1}^{H}),

(3.16g)

v_{i + 1, j + 1, k + 1}^{h} = \frac{1}{8} (v_{I, J, K}^{H} + v_{I + 1, J, K}^{H} + v_{I, J + 1, K}^{H} + v_{I, J, K + 1}^{H}) + \frac{1}{8} () v_{I + 1, J + 1, K}^{H} + v_{I + 1, J, K + 1}^{H} + v_{I, J + 1, K + 1}^{H} + v_{I + 1, J + 1, K + 1}^{H})) .

(3.16h)

The interpolated discretization error v^h corrects φ^h according to

φ^{h} \to φ^{h} + v^{h} .

(3.17)

This is step 5 of Fig. 6. Convergence is achieved when $‖ r_{n e w}^{h} - r_{o l d}^{h} ‖$ 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 = 4h and H = 8h. 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.⁸³ Schematics for “V-cycle” and “W-cycle” variants of the four-level method are provided in Fig. 7, which introduces parameters γ₀, γ₁, γ₂, and γ₃ 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 γ₀ equal to the maximum number of allowed iterations, γ₀ = 500 here. The other parameters are set to γ₁ = 2, γ₂ = 3, and γ₃ = γ₁ + γ₂, as in Ref. 83. Passing through either cycle outlined in Fig. 7, one performs CG iterations of the equation $L^{n h} v^{n h} = r^{n h}$ 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 = 8h, 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 γ₁ + γ₂ times, in either the V- or W-cycle approach.

FIG. 7.

View large Download slide

Flow diagrams illustrating four-level multigrid algorithms of either the (a) V-cycle or (b) W-cycle variety. Either approach uses a target grid of resolution h and three coarser grids of resolutions H = 2h, 4h, and 8h. Downward arrows represent the restriction of the residual error vector from a finer to a coarser grid. After relaxing γ₀ and γ₁ times, the discretization error is then interpolated from coarser to finer grids, indicated by an upward arrow. The interpolated discretization error is further relaxed on the finer grids γ₂ or γ₃ times. Convergence of the solution on the target grid is then tested, and the process is repeated until the desired accuracy is achieved.

In more detail, the four-level algorithms proceed as follows. At the target grid level, the residual error r^h corresponding to φ^h after γ₁ 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 γ₁ 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 γ₂ 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,⁴² 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³ at T = 300 K, using the amoeba force field^88,89 as implemented in the tinker software package,⁹⁰ 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₂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.⁴² For neutral water, the H₂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 Å)³ 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 γ₀ = 500, γ₁ = 2, γ₂ = 3, and γ₃ = γ₁ + γ₂. The iterative charge density is updated until the residual Δρ_iter falls below a threshold τ_PEqS = 10⁻⁵ 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,⁹¹ 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.⁹⁴

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₂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⁹⁹ 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.¹⁰⁰ 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.⁵⁰) The SCF convergence threshold is set to τ_SCF = 10⁻⁵ a.u. in all calculations, with an integral and shell-pair screening threshold of 10⁻⁸ a.u.

We use the 6-311+G* basis set for all H₂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,¹⁰¹ as appropriate.

The valence photoelectron spectrum of liquid water consists of a broad absorption feature centered at 11.3 eV,¹⁰² attributed to ionization of a 1b₁ MO localized on a single H₂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 1b₁ 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,¹⁸ 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.

FIG. 8.

View large Download slide

Spin densities ρ_spin = ρ_α − ρ_β following ionization of cluster configurations representing (a) Li⁺(aq), (b) F⁻(aq), and (c) neat liquid water. (Essentially 100% of ρ_spin is encapsulated within each surface.) In each case, a standard SCF calculation of the ionized ground state results in a hole that is localized on a water molecule near the continuum boundary and disconnected from the hydrogen-bonding network. Computed VIEs are in poor agreement with experiment. Panels (d)–(f) show the spin densities obtained from the FragMO/MOM SCF approach, which ionizes the atomic ion in (d) and (e), and a central water molecule in (f). In these cases, computed VIEs are in reasonable agreement with measured values. All VIEs were computed using the nonequilibrium IEF-PCM with a spherical cavity.

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¹⁰⁶ (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.¹⁰⁷

Here, we start from an initial guess generated using the so-called fragment MO (FragMO) procedure¹⁰⁸ 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₂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 1b₁ 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₂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₂O(aq) versus 11.3 eV (experiment);¹⁰² 11.3 eV versus 11.6 eV for F⁻(aq);¹¹ and 62.0 eV versus 60.4 eV for Li⁺(aq).¹⁸ 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 1b₁ orbital localized on a single water molecule, the higher-energy 3a₁ feature is broadened into an (unresolved) doublet,^11,102–105 assigned to ionization from bonding and anti-bonding combinations of the 3a₁ orbitals on two hydrogen-bonded H₂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 3a₁ 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),⁴² 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.

TABLE III.

Average VIEs computed with nonequilibrium PEqS and PCM methods at the (RI)MP2 level using triple-ζ basis sets as described in Sec. IV C. Computed VIEs are averages over MD snapshots, and uncertainties represent one standard deviation. Experimental error bars, which come from the references indicated, represent uncertainty in the peak position and are not peak widths.

			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₂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₂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

^a

Reference 18.

^b

Reference 102.

^c

Reference 15.

^d

Reference 11.

Agreement with experimental VIEs is excellent for the anionic solutes and for H₂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¹¹⁰ 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.¹⁸ Peak width measurements are not available for F⁻, but for Cl⁻ the measured width is 0.60 eV.¹⁸ 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,¹⁸ 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

{VIE}_{e q} = E_{fi n a l} (ε_{s o l v}) - E_{i n i t i a l} (ε_{s o l v}),

(5.1)

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,

{VIE}_{n o n e q} = E_{fi n a l} (ε_{\infty}) - E_{i n i t i a l} (ε_{s o l v}) .

(5.2)

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

\begin{aligned} Δ {VIE}_{n o n e q} & = {VIE}_{n o n e q} - {VIE}_{e q} \\ = E_{fi n a l} (ε_{\infty}) - E_{fi n a l} (ε_{s o l v}) . \end{aligned}

(5.3)

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.

TABLE IV.

Nonequilibrium solvation correction to the VIE [Eq. (5.3)], averaged over MD snapshots.

	⟨ΔVIE_noneq⟩ (eV)
	PEqS	PCM
Solute	Hybrid	Spherical	SAS
Li⁺	0.65	0.53	0.59
Na⁺	0.65	0.53	0.59
H₂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₂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.

FIG. 9.

View large Download slide

VIEs computed using the nonequilibrium PEqS method [Eq. (2.31)] along a single MD trajectory for (a) Li⁺(aq) and (b) Na⁺(aq), along with the distance d_GDS from the Gibbs dividing surface that defines the interface, again for (c) Li⁺(aq) and (d) Na⁺(aq). Also shown on the VIE plots is a best-fit line to the data (in blue) and the average bulk VIE (in red).

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.

FIG. 10.

View large Download slide

VIEs computed using the nonequilibrium PEqS method [Eq. (2.31)] along a single MD trajectory for (a) F⁻(aq) and (b) Cl⁻(aq), along with the distance d_GDS from the Gibbs dividing surface that defines the interface, again for (c) F⁻(aq) and (d) Cl⁻(aq). Also shown on the VIE plots is a best-fit line to the data (in blue) and the average bulk VIE (in red).

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₂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,⁷⁹ 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.

TABLE V.

Average distance ⟨d_GDS⟩ between the c.o.m. of the solute and the GDS, along with the corresponding average VIE, for configurations extracted from a liquid/vapor MD simulation. Uncertainties reflect one standard deviation.

Solute	⟨d_GDS⟩ (Å)	⟨VIE⟩ (eV)
H₂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.¹⁸ 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.

TABLE VI.

Errors in computed VIEs, relative to experimental values from Table III. Positive values indicate that the theoretical result is larger than the experimental VIE.

	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₂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₂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

^a

MP2 results from this work.

^b

Solute cavity determined as an isocontour of the electron density.

^c

CCSD(T)/cc-pV5Z results from Ref. 18.

^d

CCSD(T)/aug-cc-pVTZ results from Ref. 18.

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 pK_a calculations.¹¹¹) 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¹⁵ 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),¹⁵ 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.⁴²) 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.¹¹⁰

Finally we consider the 1b₁ 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₂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.¹²⁵ 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¹³¹ and is no more expensive than standard density-functional calculations. In our estimation, however, this “potential adjustors” technique¹³¹ 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⁷³ 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₂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.¹⁰⁶

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),⁴² 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₀W₀ calculations.¹³⁷

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.⁵⁰ 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.¹³⁹ J.M.H. is a fellow of the Alexander von Humboldt Foundation and serves on the Board of Directors of Q-Chem, Inc.

REFERENCES

1.

B.

Garrett

, “

Ions at the air/water interface

,”

Science

303

,

1146

(

2004

).

https://doi.org/10.1126/science.1089801

Google Scholar

Crossref

PubMed

2.

P. B.

Petersen

and

R. J.

Saykally

, “

On the nature of ions at the liquid water surface

,”

Annu. Rev. Phys. Chem.

57

,

333

(

2006

).

https://doi.org/10.1146/annurev.physchem.57.032905.104609

Google Scholar

Crossref

PubMed

3.

P.

Jungwirth

and

D. J.

Tobias

, “

Specific ion effects at the air/water interface

,”

Chem. Rev.

106

,

1259

(

2006

).

https://doi.org/10.1021/cr0403741

Google Scholar

Crossref

PubMed

4.

P.

Jungwirth

, “

Ions at aqueous interfaces

,”

Faraday Discuss.

141

,

9

(

2008

).

https://doi.org/10.1039/b816684f

Google Scholar

Crossref

5.

P.

Jungwirth

and

B.

Winter

, “

Ions at aqueous interfaces: From water surface to hydrated proteins

,”

Annu. Rev. Phys. Chem.

59

,

343

(

2008

).

https://doi.org/10.1146/annurev.physchem.59.032607.093749

Google Scholar

Crossref

PubMed

6.

P. B.

Petersen

and

R. J.

Saykally

, “

Is the liquid water surface basic or acidic? Macroscopic vs. molecular-scale investigations

,”

Chem. Phys. Lett.

458

,

255

(

2008

).

https://doi.org/10.1016/j.cplett.2008.04.010

Google Scholar

Crossref

7.

D. J.

Tobias

,

A. C.

Stern

,

M. D.

Baer

,

Y.

Levin

, and

C. J.

Mundy

, “

Simulation and theory of ions at atmospherically relevant aqueous liquid–air interfaces

,”

Annu. Rev. Phys. Chem.

64

,

339

(

2013

).

https://doi.org/10.1146/annurev-physchem-040412-110049

Google Scholar

Crossref

PubMed

8.

B.

Winter

and

M.

Faubel

, “

Photoemission from liquid aqueous solutions

,”

Chem. Rev.

106

,

1176

(

2006

).

https://doi.org/10.1021/cr040381p

Google Scholar

Crossref

PubMed

9.

R.

Seidel

,

S.

Thürmer

, and

B.

Winter

, “

Photoelectron spectroscopy meets aqueous solution: Studies from a vacuum liquid microjet

,”

J. Phys. Chem. Lett.

2

,

633

(

2011

).

https://doi.org/10.1021/jz101636y

Google Scholar

Crossref

10.

M.

Faubel

,

K. R.

Siefermann

,

Y.

Liu

, and

B.

Abel

, “

Ultrafast soft x-ray photoelectron spectroscopy at liquid water microjets

,”

Acc. Chem. Res.

45

,

120

(

2012

).

https://doi.org/10.1021/ar200154w

Google Scholar

Crossref

PubMed

11.

R.

Seidel

,

B.

Winter

, and

S.

Bradforth

, “

Valence electronic structure of aqueous solutions: Insights from photoelectron spectroscopy

,”

Annu. Rev. Phys. Chem.

67

,

283

(

2016

).

https://doi.org/10.1146/annurev-physchem-040513-103715

Google Scholar

Crossref

PubMed

12.

N.

Ottoson

,

M.

Faubel

,

S. E.

Bradforth

,

P.

Jungwirth

, and

B.

Winter

, “

Photoelectron spectroscopy of liquid water and aqueous solution: Electron effective attenuation lengths and emission-angle anisotropy

,”

J. Electron Spectrosc. Relat. Phenom.

177

,

60

(

2010

).

https://doi.org/10.1016/j.elspec.2009.08.007

Google Scholar

Crossref

13.

Y.-I.

Suzuki

,

K.

Nishizawa

,

N.

Kurahashi

, and

T.

Suzuki

, “

Effective attenuation length of an electron in liquid water between 10 and 600 eV

,”

Phys. Rev. E

90

,

010302R

(

2014

).

https://doi.org/10.1103/physreve.90.010302

Google Scholar

Crossref

14.

Y.

Yamamoto

,

S.

Karashima

,

S.

Adachi

, and

T.

Suzuki

, “

Wavelength dependence of UV photoemission from solvated electrons in bulk water, methanol, and ethanol

,”

J. Phys. Chem. A

120

,

1153

(

2016

).

https://doi.org/10.1021/acs.jpca.5b09601

Google Scholar

Crossref

PubMed

15.

D.

Luckhaus

,

Y.

Yamamoto

,

T.

Suzuki

, and

R.

Signorell

, “

Genuine binding energy of the hydrated electron

,”

Sci. Adv.

3

,

e1603224

(

2017

).

https://doi.org/10.1126/sciadv.1603224

Google Scholar

Crossref

PubMed

16.

In 2016, the experimental value for the VIE of F⁻(aq) was revised upward from 8.7 eV¹⁷ to 11.6 eV,¹¹ with much earlier theoretical calculations¹⁸ serving to validate the revised value. The F⁻(aq) VIE is difficult to discern as it is embedded in the 1b₁ band of the photoelectron spectrum of liquid water.¹¹

17.

M.

Faubel

, “

Photoelectron spectroscopy at liquid surfaces

,” in

Photoionization and Photodetachment

, Volume 10A of Advanced Series in Physical Chemistry, edited by

C. Y.

Ng

(

World Scientific

,

2000

), pp.

634

–

690

.

Google Scholar

Crossref

18.

B.

Winter

,

R.

Weber

,

I. V.

Hertel

,

M.

Faubel

,

P.

Jungwirth

,

E. C.

Brown

, and

S. E.

Bradforth

, “

Electron binding energies of aqueous alkali and halide ions: EUV photoelectron spectroscopy of liquid solutions and combined ab initio and molecular dynamics calculations

,”

J. Am. Chem. Soc.

127

,

7203

(

2005

).

https://doi.org/10.1021/ja042908l

Google Scholar

Crossref

PubMed

19.

J.

Tomasi

,

B.

Mennucci

, and

R.

Cammi

, “

Quantum mechanical continuum solvation models

,”

Chem. Rev.

105

,

2999

(

2005

).

https://doi.org/10.1021/cr9904009

Google Scholar

Crossref

PubMed

20.

C. J.

Cramer

and

D. G.

Truhlar

, “

A universal approach to solvation modeling

,”

Acc. Chem. Res.

41

,

760

(

2008

).

https://doi.org/10.1021/ar800019z

Google Scholar

Crossref

PubMed

21.

B.

Mennucci

, “

Polarizable continuum model

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

2

,

386

(

2012

).

https://doi.org/10.1002/wcms.1086

Google Scholar

Crossref

22.

J. M.

Herbert

and

A. W.

Lange

, “

Polarizable continuum models for (bio)molecular electrostatics: Basic theory and recent developments for macromolecules and simulations

,” in

Many-Body Effects and Electrostatics in Biomolecules

, edited by

Q.

Cui

,

P.

Ren

, and

M.

Meuwly

(

Pan Stanford

,

2016

), Chap. 11, pp.

363

–

416

.

Google Scholar

23.

A.

Klamt

, “

The COSMO and COSMO-RS solvation models

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

8

,

e1338

(

2018

).

https://doi.org/10.1002/wcms.1338

Google Scholar

Crossref

24.

A. V.

Marenich

,

C. J.

Cramer

, and

D. G.

Truhlar

, “

Universal solvation model based on solute electron density and on a continuum model of the solvent defined by the bulk dielectric constant and atomic surface tensions

,”

J. Phys. Chem. B

113

,

6378

(

2009

).

https://doi.org/10.1021/jp810292n

Google Scholar

Crossref

PubMed

25.

A.

Klamt

,

B.

Mennucci

,

J.

Tomasi

,

V.

Barone

,

C.

Curutchet

,

M.

Orozco

, and

F. J.

Luque

, “

On the performance of continuum solvation methods. A comment on ‘Universal approaches to solvation modeling

,’”

Acc. Chem. Res.

42

,

489

(

2009

).

https://doi.org/10.1021/ar800187p

Google Scholar

Crossref

PubMed

26.

A.

Pomogaeva

and

D. M.

Chipman

, “

Hydration energy from a composite method for implicit representation of the solvent

,”

J. Chem. Theory Comput.

10

,

211

(

2014

).

https://doi.org/10.1021/ct400894j

Google Scholar

Crossref

PubMed

27.

Z.-Q.

You

and

J. M.

Herbert

, “

Reparameterization of an accurate, few-parameter implicit solvation model for quantum chemistry: Composite method for implicit representation of solvent, CMIRS v. 1.1

,”

J. Chem. Theory Comput.

12

,

4338

(

2016

).

https://doi.org/10.1021/acs.jctc.6b00644

Google Scholar

Crossref

PubMed

28.

B.

Mennucci

,

E.

Cancés

, and

J.

Tomasi

, “

Evaluation of solvent effects in isotropic and anisotropic dielectrics and in ionic solutions with a unified integral equation method: Theoretical bases, computational implementation, and numerical applications

,”

J. Phys. Chem. B

101

,

10506

(

1997

).

https://doi.org/10.1021/jp971959k

Google Scholar

Crossref

29.

E.

Cancés

,

B.

Mennucci

, and

J.

Tomasi

, “

A new integral equation formalism for the polarizable continuum model: Theoretical background and applications to isotropic and anisotropic dielectrics

,”

J. Chem. Phys.

107

,

3032

(

1997

).

https://doi.org/10.1063/1.474659

Google Scholar

Crossref

30.

E.

Cancès

and

B.

Mennucci

, “

New applications of integral equations methods for solvation continuum models: Ionic solutions and liquid crystals

,”

J. Math. Chem.

23

,

309

(

1998

).

https://doi.org/10.1023/a:1019133611148

Google Scholar

Crossref

31.

B.

Mennucci

and

R.

Cammi

, “

Ab initio model to predict NMR shielding tensors for solutes in liquid crystals

,”

Int. J. Quantum Chem.

93

,

121

(

2003

).

https://doi.org/10.1002/qua.10541

Google Scholar

Crossref

32.

L.

Frediani

,

B.

Mennucci

, and

R.

Cammi

, “

Quantum-mechanical continuum solvation study of the polarizability of halides at the water/air interface

,”

J. Phys. Chem. B

108

,

13796

(

2004

).

https://doi.org/10.1021/jp048414s

Google Scholar

Crossref

33.

L.

Bondesson

,

L.

Frediani

,

H.

Ågren

, and

B.

Mennucci

, “

Solvation of

.

N_{3}^{-}

at the water surface: The polarizable continuum model approach

,”

J. Phys. Chem. B

110

,

11361

(

2006

).

https://doi.org/10.1021/jp060794p

Google Scholar

Crossref

PubMed

34.

K.

Mozgawa

,

B.

Mennucci

, and

L.

Frediani

, “

Solvation at surfaces and interfaces: A quantum-mechanical/continuum approach including nonelectrostatic contributions

,”

J. Phys. Chem. C

118

,

4715

(

2014

).

https://doi.org/10.1021/jp4117276

Google Scholar

Crossref

35.

D.

Si

and

H.

Li

, “

Heterogeneous conductorlike solvation model

,”

J. Chem. Phys.

131

,

044123

(

2009

).

https://doi.org/10.1063/1.3187527

Google Scholar

Crossref

PubMed

36.

J.-B.

Wang

,

J.-Y.

Ma

, and

X.-Y.

Li

, “

Polarizable continuum model associated with the self-consistent-reaction field for molecular adsorbates at the interface

,”

Phys. Chem. Chem. Phys.

12

,

207

(

2010

).

https://doi.org/10.1039/b914652k

Google Scholar

Crossref

PubMed

37.

H.

Hoshi

,

M.

Sakurai

,

Y.

Inoue

, and

R.

Chûjô

, “

Medium effects on the molecular electronic structure. I. The formulation of a theory for the estimation of a molecular electronic structure surrounded by an anisotropic medium

,”

J. Chem. Phys.

87

,

1107

(

1987

).

https://doi.org/10.1063/1.453343

Google Scholar

Crossref

38.

H.

Hoshi

,

M.

Sakurai

,

Y.

Inoue

, and

R.

Chûjô

, “

Medium effects on the molecular electronic structure: Part 2. The application of the theory of medium effects in the framework of the CNDO and INDO methods

,”

J. Mol. Struct.: THEOCHEM

180

,

267

(

1988

).

https://doi.org/10.1016/0166-1280(88)80094-0

Google Scholar

Crossref

39.

A. H.

Boschitsch

and

M. O.

Fenley

, “

A fast and robust Poisson–Boltzmann solver based on adaptive Cartesian grids

,”

J. Chem. Theory Comput.

7

,

1524

(

2011

).

https://doi.org/10.1021/ct1006983

Google Scholar

Crossref

PubMed

40.

O.

Andreussi

,

I.

Dabo

, and

N.

Marzari

, “

Revised self-consistent continuum solvation in electronic-structure calculations

,”

J. Chem. Phys.

136

,

064102

(

2012

).

https://doi.org/10.1063/1.3676407

Google Scholar

Crossref

PubMed

41.

G.

Fisicaro

,

L.

Genovese

,

O.

Andreussi

,

N.

Marzari

, and

S.

Goedecker

, “

A generalized Poisson and Poisson-Boltzmann solver for electrostatic environments

,”

J. Chem. Phys.

144

,

014103

(

2016

).

https://doi.org/10.1063/1.4939125

Google Scholar

Crossref

PubMed

42.

M. P.

Coons

,

Z.-Q.

You

, and

J. M.

Herbert

, “

The hydrated electron at the surface of neat liquid water appears to be indistinguishable from the bulk species

,”

J. Am. Chem. Soc.

138

,

10879

(

2016

).

https://doi.org/10.1021/jacs.6b06715

Google Scholar

Crossref

PubMed

43.

J. C.

Womack

,

L.

Anton

,

J.

Dziedzic

,

P. J.

Hasnip

,

M. I. J.

Probert

, and

C.-K.

Skylaris

, “

DL_MG: A parallel multigrid Poisson and Poisson–Boltzmann solver for electronic structure calculations in vacuum and solution

,”

J. Chem. Theory Comput.

14

,

1412

(

2018

).

https://doi.org/10.1021/acs.jctc.7b01274

Google Scholar

Crossref

PubMed

44.

M. A.

Aguilar

,

F. J.

Olivares del Valle

, and

J.

Tomasi

, “

Nonequilibrium solvation: An ab initio quantum-mechanical method in the continuum cavity model approximation

,”

J. Chem. Phys.

98

,

7375

(

1993

).

https://doi.org/10.1063/1.464728

Google Scholar

Crossref

45.

R.

Cammi

and

J.

Tomasi

, “

Nonequilibrium solvation theory for the polarizable continuum model: A new formulation at the SCF level with application to the case of the frequency-dependent linear electric response function

,”

Int. J. Quantum Chem., Symp.

29

,

465

(

1995

).

Google Scholar

46.

D. M.

Chipman

, “

Vertical electronic excitation with a dielectric continuum model of solvation including volume polarization. I. Theory

,”

J. Chem. Phys.

131

,

014103

(

2009

).

https://doi.org/10.1063/1.3157464

Google Scholar

Crossref

PubMed

47.

A. V.

Marenich

,

C. J.

Cramer

,

D. G.

Truhlar

,

C. A.

Guido

,

B.

Mennucci

,

G.

Scalmani

, and

M. J.

Frisch

, “

Practical computation of electronic excitation in solution: Vertical excitation model

,”

Chem. Sci.

2

,

2143

(

2011

).

https://doi.org/10.1039/c1sc00313e

Google Scholar

Crossref

48.

J.-M.

Mewes

,

Z.-Q.

You

,

M.

Wormit

,

T.

Kriesche

,

J. M.

Herbert

, and

A.

Dreuw

, “

Experimental benchmark data and systematic evaluation of two a posteriori, polarizable-continuum corrections for vertical excitation energies in solution

,”

J. Phys. Chem. A

119

,

5446

(

2015

).

https://doi.org/10.1021/jp511163y

Google Scholar

Crossref

PubMed

49.

Z.-Q.

You

,

J.-M.

Mewes

,

A.

Dreuw

, and

J. M.

Herbert

, “

Comparison of the Marcus and Pekar partitions in the context of non-equilibrium, polarizable-continuum reaction-field solvation models

,”

J. Chem. Phys.

143

,

204104

(

2015

).

https://doi.org/10.1063/1.4936357

Google Scholar

Crossref

PubMed

50.

J.-M.

Mewes

,

J. M.

Herbert

, and

A.

Dreuw

, “

On the accuracy of the state-specific polarizable continuum model for the description of correlated ground and excited states in solution

,”

Phys. Chem. Chem. Phys.

19

,

1644

(

2017

).

https://doi.org/10.1039/c6cp05986d

Google Scholar

Crossref

PubMed

51.

M. R.

Provorse

,

T.

Peev

,

C.

Xiong

, and

C. M.

Isborn

, “

Convergence of excitation energies in mixed quantum and classical solvent: Comparison of continuum and point charge models

,”

J. Phys. Chem. B

120

,

12148

(

2016

).

https://doi.org/10.1021/acs.jpcb.6b09176

Google Scholar

Crossref

PubMed

52.

Y.

Yamamoto

,

Y.-I.

Suzuki

,

G.

Tomasello

,

T.

Horio

,

S.

Karashima

,

R.

Mitríc

, and

T.

Suzuki

, “

Time- and angle-resolved photoemission spectroscopy of hydrated electrons near a liquid water surface

,”

Phys. Rev. Lett.

112

,

187603

(

2014

).

https://doi.org/10.1103/physrevlett.112.187603

Google Scholar

Crossref

PubMed

53.

A. H. C.

West

,

B. L.

Yoder

,

D.

Luckhaus

,

C.-M.

Saak

,

M.

Doppelbauer

, and

R.

Signorell

, “

Angle-resolved photoemission of solvated electrons in sodium-doped clusters

,”

J. Phys. Chem. Lett.

6

,

1487

(

2015

).

https://doi.org/10.1021/acs.jpclett.5b00477

Google Scholar

Crossref

PubMed

54.

S.

Hartweg

,

B. L.

Yoder

,

G. A.

Garcia

,

L.

Nahon

, and

R.

Signorell

, “

Size-resolved photoelectron anisotropy of gas phase water clusters and predictions for liquid water

,”

Phys. Rev. Lett.

118

,

103402

(

2017

).

https://doi.org/10.1103/physrevlett.118.103402

Google Scholar

Crossref

PubMed

55.

J.

Nishitani

,

C. W.

West

, and

T.

Suzuki

, “

Angle-resolved photoemission spectroscopy of liquid water at 29.5 eV

,”

Struct. Dyn.

4

,

044014

(

2017

).

https://doi.org/10.1063/1.4979857

Google Scholar

Crossref

PubMed

56.

J.

Tomasi

and

M.

Persico

, “

Molecular interactions in solution: An overview of methods based on continuous distributions of the solvent

,”

Chem. Rev.

94

,

2027

(

1994

).

https://doi.org/10.1021/cr00031a013

Google Scholar

Crossref

57.

M.

Cossi

and

V.

Barone

, “

Separation between fast and slow polarizations in continuum solvation models

,”

J. Phys. Chem. A

104

,

10614

(

2000

).

https://doi.org/10.1021/jp000997s

Google Scholar

Crossref

58.

R.

Improta

,

V.

Barone

,

G.

Scalmani

, and

M. J.

Frisch

, “

A state-specific polarizable continuum model time dependent density functional method for excited state calculations in solution

,”

J. Chem. Phys.

125

,

054103

(

2006

).

https://doi.org/10.1063/1.2222364

Google Scholar

Crossref

PubMed

59.

L. D.

Jacobson

and

J. M.

Herbert

, “

A simple algorithm for determining orthogonal, self-consistent excited-state wave functions for a state-specific Hamiltonian: Application to the optical spectrum of the aqueous electron

,”

J. Chem. Theory Comput.

7

,

2085

(

2011

).

https://doi.org/10.1021/ct200265t

Google Scholar

Crossref

PubMed

60.

R. A.

Marcus

, “

On the theory of oxidation-reduction reactions involving electron transfer. I

,”

J. Chem. Phys.

24

,

966

(

1956

).

https://doi.org/10.1063/1.1742723

Google Scholar

Crossref

61.

J. A.

Grant

,

B. T.

Pickup

, and

A.

Nicholls

, “

A smooth permittivity function for Poisson–Boltzmann solvation methods

,”

J. Comput. Chem.

22

,

608

(

2001

).

https://doi.org/10.1002/jcc.1032

Google Scholar

Crossref

62.

J.-L.

Fattebert

and

F.

Gygi

, “

Density functional theory for efficieint ab initio molecular dynamics simulations in solution

,”

J. Comput. Chem.

23

,

662

(

2002

).

https://doi.org/10.1002/jcc.10069

Google Scholar

Crossref

PubMed

63.

J.-L.

Fattebert

and

F.

Gygi

, “

First-principles molecular dynamics simulations in a continuum solvent

,”

Int. J. Quantum Chem.

93

,

139

(

2003

).

https://doi.org/10.1002/qua.10548

Google Scholar

Crossref

64.

J.

Dziedzic

,

H. H.

Helal

,

C.-K.

Skylaris

,

A. A.

Mostofi

, and

M. C.

Payne

, “

Minimal parameter implicit solvent model for ab initio electronic-structure calculations

,”

Europhys. Lett.

95

,

43001

(

2011

).

https://doi.org/10.1209/0295-5075/95/43001

Google Scholar

Crossref

65.

J. B.

Foresman

,

T. A.

Keith

,

K. B.

Wiberg

,

J.

Snoonian

, and

M. J.

Frisch

, “

Solvent effects. 5. Influence of cavity shape, truncation of electrostatics, and electron correlation on ab initio reaction field calculations

,”

J. Phys. Chem.

100

,

16098

(

1996

).

https://doi.org/10.1021/jp960488j

Google Scholar

Crossref

66.

C.-G.

Zhan

and

D. M.

Chipman

, “

Cavity size in reaction field theory

,”

J. Chem. Phys.

109

,

10543

(

1998

).

https://doi.org/10.1063/1.477755

Google Scholar

Crossref

67.

F.

Chen

and

D. M.

Chipman

, “

Boundary element methods for dielectric cavity construction and integration

,”

J. Chem. Phys.

119

,

10289

(

2003

).

https://doi.org/10.1063/1.1615232

Google Scholar

Crossref

68.

A.

Bondi

, “

van der Waals volumes and radii

,”

J. Phys. Chem.

68

,

441

(

1964

).

https://doi.org/10.1021/j100785a001

Google Scholar

Crossref

69.

R. S.

Rowland

and

R.

Taylor

, “

Intermolecular nonbonded contact distances in organic crystal structures: Comparison with distances expected from van der Waals radii

,”

J. Phys. Chem.

100

,

7384

(

1996

).

https://doi.org/10.1021/jp953141+

Google Scholar

Crossref

70.

J. M.

Herbert

and

M. P.

Coons

, “

The hydrated electron

,”

Annu. Rev. Phys. Chem.

68

,

447

(

2017

).

https://doi.org/10.1146/annurev-physchem-052516-050816

Google Scholar

Crossref

PubMed

71.

F.

Uhlig

,

O.

Marsalek

, and

P.

Jungwirth

, “

Unraveling the complex nature of the hydrated electron

,”

J. Phys. Chem. Lett.

3

,

3071

(

2012

).

https://doi.org/10.1021/jz301449f

Google Scholar

Crossref

PubMed

72.

A.

Kumar

,

J. A.

Walker

,

D. M.

Bartels

, and

M. D.

Sevilla

, “

A simple ab initio model for the hydrated electron that matches experiment

,”

J. Phys. Chem. A

119

,

9148

(

2015

).

https://doi.org/10.1021/acs.jpca.5b04721

Google Scholar

Crossref

PubMed

73.

Chemoinformatics: A Textbook

, edited by

J.

Gasteiger

and

T.

Engel

(

Wiley-VCH

,

Weinheim

,

2003

).

Google Scholar

Crossref

74.

D. H.

Brookes

and

T.

Head-Gordon

, “

Family of oxygen–oxygen radial distribution functions for water

,”

J. Phys. Chem. Lett.

6

,

2938

(

2015

).

https://doi.org/10.1021/acs.jpclett.5b01066

Google Scholar

Crossref

PubMed

75.

L. D.

Jacobson

and

J. M.

Herbert

, “

Theoretical characterization of four distinct isomer types in hydrated-electron clusters, and proposed assignments for photoelectron spectra of water cluster anions

,”

J. Am. Chem. Soc.

133

,

19889

(

2011

).

https://doi.org/10.1021/ja208024p

Google Scholar

Crossref

PubMed

76.

L.

Frediani

,

R.

Cammi

,

S.

Corni

, and

J.

Tomasi

, “

A polarizable continuum model for molecules at diffuse interfaces

,”

J. Chem. Phys.

120

,

3893

(

2004

).

https://doi.org/10.1063/1.1643727

Google Scholar

Crossref

PubMed

77.

M.-H.

Ho

,

M. L.

Klein

, and

I.-F.

Kuo

, “

Bulk and interfacial aqueous fluoride: An investigation via first principles molecular dynamics

,”

J. Phys. Chem. A

113

,

2070

(

2009

).

https://doi.org/10.1021/jp808735x

Google Scholar

Crossref

PubMed

78.

G. L.

Warren

and

S.

Patel

, “

Electrostatic properties of aqueous salt solution interfaces: A comparison of polarizable and nonpolarizable ion models

,”

J. Phys. Chem. B

112

,

11679

(

2008

).

https://doi.org/10.1021/jp8038835

Google Scholar

Crossref

PubMed

79.

F.

Uhlig

,

O.

Marsalek

, and

P.

Jungwirth

, “

Electron at the surface of water: Dehydrated or not?

,”

J. Phys. Chem. Lett.

4

,

338

(

2013

).

https://doi.org/10.1021/jz3020953

Google Scholar

Crossref

PubMed

80.

G.

Scalmani

,

V.

Barone

,

K. N.

Kudin

,

C. S.

Pomelli

,

G. E.

Scuseria

, and

M. J.

Frisch

, “

Achieving linear-scaling computational cost for the polarizable continuum model of solvation

,”

Theor. Chem. Acc.

111

,

90

(

2004

).

https://doi.org/10.1007/s00214-003-0527-2

Google Scholar

Crossref

81.

F.

Lipparini

,

B.

Stamm

,

E.

Cancès

,

Y.

Maday

, and

B.

Mennucci

, “

Fast domain decomposition algorithm for continuum solvation models: Energy and first derivatives

,”

J. Chem. Theory Comput.

9

,

3637

(

2013

).

https://doi.org/10.1021/ct400280b

Google Scholar

Crossref

PubMed

82.

F.

Lipparini

,

L.

Lagardère

,

G.

Scalmani

,

B.

Stamm

,

E.

Cancès

,

Y.

Maday

,

J.-P.

Piquemal

,

M. J.

Frisch

, and

B.

Mennucci

, “

Quantum calculations in solution for large to very large molecules: A new linear scaling QM/continuum approach

,”

J. Phys. Chem. Lett.

5

,

953

(

2014

).

https://doi.org/10.1021/jz5002506

Google Scholar

Crossref

PubMed

83.

C. H.

Venner

and

A. A.

Lubrecht

,

Multi-Level Methods in Lubrication

, Volume 37 of Tribology, 1st ed. (

Elsevier

,

2000

).

Google Scholar

84.

L. D.

Jacobson

and

J. M.

Herbert

, “

Polarization-bound quasi-continuum states are responsible for the ‘blue tail’ in the optical absorption spectrum of the aqueous electron

,”

J. Am. Chem. Soc.

132

,

10000

(

2010

).

https://doi.org/10.1021/ja1042484

Google Scholar

Crossref

PubMed

85.

L. D.

Jacobson

and

J. M.

Herbert

, “

A one-electron model for the aqueous electron that includes many-body electron-water polarization: Bulk equilibrium structure, vertical electron binding energy, and optical absorption spectrum

,”

J. Chem. Phys.

133

,

154506

(

2010

).

https://doi.org/10.1063/1.3490479

Google Scholar

Crossref

PubMed

86.

J. M.

Herbert

and

L. D.

Jacobson

, “

Structure of the aqueous electron: Assessment of one-electron pseudopotential models in comparison to experimental data and time-dependent density functional theory

,”

J. Phys. Chem. A

115

,

14470

(

2011

).

https://doi.org/10.1021/jp206391d

Google Scholar

Crossref

PubMed

87.

F.

Uhlig

,

J. M.

Herbert

,

M. P.

Coons

, and

P.

Jungwirth

, “

Optical spectroscopy of the bulk and interfacial hydrated electron from ab initio calculations

,”

J. Phys. Chem. A

118

,

7507

(

2014

).

https://doi.org/10.1021/jp5004243

Google Scholar

Crossref

PubMed

88.

P.

Ren

and

J. W.

Ponder

, “

Polarizable atomic multipole water model for molecular mechanics simulation

,”

J. Phys. Chem. B

107

,

5933

(

2003

).

https://doi.org/10.1021/jp027815+

Google Scholar

Crossref

89.

A.

Grossfield

,

P.

Ren

, and

J. W.

Ponder

, “

Ion solvation thermodynamics from simulation with a polarizable force field

,”

J. Am. Chem. Soc.

125

,

15671

(

2003

).

https://doi.org/10.1021/ja037005r

Google Scholar

Crossref

PubMed

90.

M.

Harger

,

D.

Li

,

Z.

Wang

,

K.

Dalby

,

L.

Lagardère

,

J.-P.

Piquemal

,

J.

Ponder

, and

P.

Ren

, “

Tinker-OpenMM: Absolute and relative alchemical free energies using AMOEBA on GPUs

,”

J. Comput. Chem.

38

,

2047

(

2017

).

https://doi.org/10.1002/jcc.24853

Google Scholar

Crossref

PubMed

91.

A. W.

Lange

and

J. M.

Herbert

, “

Symmetric versus asymmetric discretization of the integral equations in polarizable continuum solvation models

,”

Chem. Phys. Lett.

509

,

77

(

2011

).

https://doi.org/10.1016/j.cplett.2011.04.092

Google Scholar

Crossref

92.

E.

Cancès

and

B.

Mennucci

, “

Comment on ‘Reaction field treatment of charge penetration’ [J. Chem. Phys. 112, 5558 (2000)]

,”

J. Chem. Phys.

114

,

4744

(

2001

).

https://doi.org/10.1063/1.1349091

Google Scholar

Crossref

93.

D. M.

Chipman

, “

Comparison of solvent reaction field representations

,”

Theor. Chem. Acc.

107

,

80

(

2002

).

https://doi.org/10.1007/s00214-001-0302-1

Google Scholar

Crossref

94.

A. W.

Lange

and

J. M.

Herbert

, “

A smooth, non-singular, and faithful discretization scheme for polarizable continuum models: The switching/Gaussian approach

,”

J. Chem. Phys.

133

,

244111

(

2010

).

https://doi.org/10.1063/1.3511297

Google Scholar

Crossref

PubMed

95.

C.-G.

Zhan

,

J.

Bentley

, and

D. M.

Chipman

, “

Volume polarization in reaction field theory

,”

J. Chem. Phys.

108

,

177

(

1998

).

https://doi.org/10.1063/1.475371

Google Scholar

Crossref

96.

D. M.

Chipman

, “

Simulation of volume polarization in reaction field theory

,”

J. Chem. Phys.

110

,

8012

(

1999

).

https://doi.org/10.1063/1.478729

Google Scholar

Crossref

97.

D. M.

Chipman

, “

Reaction field treatment of charge penetration

,”

J. Chem. Phys.

112

,

5558

(

2000

).

https://doi.org/10.1063/1.481133

Google Scholar

Crossref

98.

M.

Cossi

,

N.

Rega

,

G.

Scalmani

, and

V.

Barone

, “

Polarizable dielectric continuum model of solvation with inclusion of charge penetration effects

,”

J. Chem. Phys.

114

,

5691

(

2001

).

https://doi.org/10.1063/1.1354187

Google Scholar

Crossref

99.

Y.

Shao

,

Z.

Gan

,

E.

Epifanovsky

,

A. T. B.

Gilbert

,

M.

Wormit

,

J.

Kussmann

,

A. W.

Lange

,

A.

Behn

,

J.

Deng

,

X.

Feng

,

D.

Ghosh

,

M.

Goldey

,

P. R.

Horn

,

L. D.

Jacobson

,

I.

Kaliman

,

R. Z.

Khaliullin

,

T.

Kús

,

A.

Landau

,

J.

Liu

,

E. I.

Proynov

,

Y. M.

Rhee

,

R. M.

Richard

,

M. A.

Rohrdanz

,

R. P.

Steele

,

E. J.

Sundstrom

,

H. L.

Woodcock

III,

P. M.

Zimmerman

,

D.

Zuev

,

B.

Albrecht

,

E.

Alguire

,

B.

Austin

,

G. J. O.

Beran

,

Y. A.

Bernard

,

E.

Berquist

,

K.

Brandhorst

,

K. B.

Bravaya

,

S. T.

Brown

,

D.

Casanova

,

C.-M.

Chang

,

Y.

Chen

,

S. H.

Chien

,

K. D.

Closser

,

D. L.

Crittenden

,

M.

Diedenhofen

,

R. A.

DiStasi

, Jr.,

H.

Do

,

A. D.

Dutoi

,

R. G.

Edgar

,

S.

Fatehi

,

L.

Fusti-Molnar

,

A.

Ghysels

,

A.

Golubeva-Zadorozhnaya

,

J.

Gomes

,

M. W. D.

Hanson-Heine

,

P. H. P.

Harbach

,

A. W.

Hauser

,

E. G.

Hohenstein

,

Z. C.

Holden

,

T.-C.

Jagau

,

H.

Ji

,

B.

Kaduk

,

K.

Khistyaev

,

J.

Kim

,

J.

Kim

,

R. A.

King

,

P.

Klunzinger

,

D.

Kosenkov

,

T.

Kowalczyk

,

C. M.

Krauter

,

K. U.

Lao

,

A.

Laurent

,

K. V.

Lawler

,

S. V.

Levchenko

,

C. Y.

Lin

,

F.

Liu

,

E.

Livshits

,

R. C.

Lochan

,

A.

Luenser

,

P.

Manohar

,

S. F.

Manzer

,

S.-P.

Mao

,

N.

Mardirossian

,

A. V.

Marenich

,

S. A.

Maurer

,

N. J.

Mayhall

,

C. M.

Oana

,

R.

Olivares-Amaya

,

D. P.

O’Neill

,

J. A.

Parkhill

,

T. M.

Perrine

,

R.

Peverati

,

P. A.

Pieniazek

,

A.

Prociuk

,

D. R.

Rehn

,

E.

Rosta

,

N. J.

Russ

,

N.

Sergueev

,

S. M.

Sharada

,

S.

Sharma

,

D. W.

Small

,

A.

Sodt

,

T.

Stein

,

D.

Stück

,

Y.-C.

Su

,

A. J. W.

Thom

,

T.

Tsuchimochi

,

L.

Vogt

,

O.

Vydrov

,

T.

Wang

,

M. A.

Watson

,

J.

Wenzel

,

A.

White

,

C. F.

Williams

,

V.

Vanovschi

,

S.

Yeganeh

,

S. R.

Yost

,

Z.-Q.

You

,

I. Y.

Zhang

,

X.

Zhang

,

Y.

Zhao

,

B. R.

Brooks

,

G. K. L.

Chan

,

D. M.

Chipman

,

C. J.

Cramer

,

W. A.

Goddard

III,

M. S.

Gordon

,

W. J.

Hehre

,

A.

Klamt

,

H. F.

Schaefer

III,

M. W.

Schmidt

,

C. D.

Sherrill

,

D. G.

Truhlar

,

A.

Warshel

,

X.

Xu

,

A.

Aspuru-Guzik

,

R.

Baer

,

A. T.

Bell

,

N. A.

Besley

,

J.-D.

Chai

,

A.

Dreuw

,

B. D.

Dunietz

,

T. R.

Furlani

,

S. R.

Gwaltney

,

C.-P.

Hsu

,

Y.

Jung

,

J.

Kong

,

D. S.

Lambrecht

,

W.

Liang

,

C.

Ochsenfeld

,

V. A.

Rassolov

,

L. V.

Slipchenko

,

J. E.

Subotnik

,

T.

Van Voorhis

,

J. M.

Herbert

,

A. I.

Krylov

,

P. M. W.

Gill

, and

M.

Head-Gordon

, “

Advances in molecular quantum chemistry contained in the Q-Chem 4 program package

,”

Mol. Phys.

113

,

184

(

2015

).

https://doi.org/10.1080/00268976.2014.952696

Google Scholar

Crossref

100.

R. A.

Kendall

and

H. A.

Früchtl

, “

The impact of the resolution of the identity approximate integral method on modern ab initio algorithm development

,”

Theor. Chem. Acc.

97

,

158

(

1997

).

https://doi.org/10.1007/s002140050249

Google Scholar

Crossref

101.

F.

Weigend

,

A.

Köhn

, and

C.

Hättig

, “

Efficient use of the correlation consistent basis sets in resolution of the identity MP2 calculations

,”

J. Chem. Phys.

116

,

3175

(

2002

).

https://doi.org/10.1063/1.1445115

Google Scholar

Crossref

102.

N.

Kurahashi

,

S.

Karashima

,

Y.

Tang

,

T.

Horio

,

B.

Abulimiti

,

Y.-I.

Suzuki

,

Y.

Ogi

,

M.

Oura

, and

T.

Suzuki

, “

Photoelectron spectroscopy of aqueous solutions: Streaming potentials of NaX (X = Cl, Br, and I) solutions and electron binding energies of liquid water and X⁻

,”

J. Chem. Phys.

140

,

174506

(

2014

).

https://doi.org/10.1063/1.4871877

Google Scholar

Crossref

PubMed

103.

M.

Faubel

,

B.

Steiner

, and

J. P.

Toennies

, “

Photoelectron spectroscopy of liquid water, some alcohols, and pure nonane in free micro jets

,”

J. Chem. Phys.

106

,

9013

(

1997

).

https://doi.org/10.1063/1.474034

Google Scholar

Crossref

104.

B.

Winter

,

R.

Weber

,

W.

Widdra

,

M.

Dittman

,

M.

Faubel

, and

I. V.

Hertel

, “

Full valence band photoemission from liquid water using EUV synchrotron radiation

,”

J. Phys. Chem. A

108

,

2625

(

2004

).

https://doi.org/10.1021/jp030263q

Google Scholar

Crossref

105.

K.

Nishizawa

,

N.

Kurahashi

,

K.

Sekiguchi

,

T.

Mizuno

, and

Y.

Ogi

, “

High-resolution soft x-ray photoelectron spectroscopy of liquid water

,”

Phys. Chem. Chem. Phys.

13

,

413

(

2011

).

https://doi.org/10.1039/c0cp01636e

Google Scholar

Crossref

PubMed

106.

A. T. B.

Gilbert

,

N. A.

Besley

, and

P. M. W.

Gill

, “

Self-consistent field calculations of excited states using the maximum overlap method (MOM)

,”

J. Phys. Chem. A

112

,

13164

(

2008

).

https://doi.org/10.1021/jp801738f

Google Scholar

Crossref

PubMed

107.

N. A.

Besley

,

A. T. B.

Gilbert

, and

P. M. W.

Gill

, “

Self-consistent-field calculations of core excited states

,”

J. Chem. Phys.

130

,

124308

(

2009

).

https://doi.org/10.1063/1.3092928

Google Scholar

Crossref

PubMed

108.

The “FragMO” SCF guess discussed herein should not be confused with the fragment molecular orbital method, as described, for example, in

D. G.

Fedorov

, “

The fragment molecular orbital method: Theoretical development, implementation in GAMESS, and applications

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

7

,

e1322

(

2017

).

https://doi.org/10.1002/wcms.1322

Google Scholar

Crossref

109.

D.

Nordlund

,

M.

Odelius

,

H.

Bluhm

,

H.

Ogasawara

,

L. G. M.

Pettersson

, and

A.

Nilsson

, “

Electronic structure effects in liquid water studied by photoelectron spectroscopy and density functional theory

,”

Chem. Phys. Lett.

460

,

86

(

2008

).

https://doi.org/10.1016/j.cplett.2008.04.096

Google Scholar

Crossref

110.

Reference 72 considers a four-coordinate, tetrahedral ${(H_{2} O)}_{4}^{-}$ 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.¹⁵ 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.

111.

R.

Casasnovas

,

J.

Ortega-Castro

,

J.

Frau

,

J.

Donoso

, and

F.

Muñoz

, “

Theoretical pK_a calculations with continuum solvents, alternative protocols to thermodynamic cycles

,”

Int. J. Quantum Chem.

114

,

1350

(

2014

).

https://doi.org/10.1002/qua.24699

Google Scholar

Crossref

112.

J. V.

Coe

,

S. M.

Williams

, and

K. H.

Bowen

, “

Photoelectron spectra of hydrated electron clusters vs. cluster size

,”

Int. Rev. Phys. Chem.

27

,

27

(

2008

).

https://doi.org/10.1080/01442350701783543

Google Scholar

Crossref

113.

Y.

Tang

,

H.

Shen

,

K.

Sekiguchi

,

N.

Kurahashi

,

T.

Mizuno

,

Y. I.

Suzuki

, and

T.

Suzuki

, “

Direct measurement of vertical binding energy of a hydrated electron

,”

Phys. Chem. Chem. Phys.

12

,

3653

(

2010

).

https://doi.org/10.1039/b925741a

Google Scholar

Crossref

PubMed

114.

K. R.

Siefermann

,

Y.

Liu

,

E.

Lugovoy

,

O.

Link

,

M.

Faubel

,

U.

Buck

,

B.

Winter

, and

B.

Abel

, “

Binding energies, lifetimes and implications of bulk and interface solvated electrons in water

,”

Nat. Phys.

2

,

274

(

2010

).

https://doi.org/10.1038/nchem.580

Google Scholar

Crossref

115.

T.

Horio

,

H.

Shen

,

S.

Adachi

, and

T.

Suzuki

, “

Photoelectron spectra of solvated electrons in bulk water, methanol, and ethanol

,”

Chem. Phys. Lett.

535

,

12

(

2012

).

https://doi.org/10.1016/j.cplett.2012.03.051

Google Scholar

Crossref

116.

F.

Buchner

,

T.

Schultz

, and

A.

Lübcke

, “

Solvated electrons at the water-air interface: Surface versus bulk signal in low kinetic energy photoelectron spectroscopy

,”

Phys. Chem. Chem. Phys.

14

,

5837

(

2012

).

https://doi.org/10.1039/c2cp23305c

Google Scholar

Crossref

PubMed

117.

R. E.

Larsen

,

W. J.

Glover

, and

B. J.

Schwartz

, “

Does the hydrated electron occupy a cavity?

,”

Science

329

,

65

(

2010

).

https://doi.org/10.1126/science.1189588

Google Scholar

Crossref

PubMed

118.

L. D.

Jacobson

and

J. M.

Herbert

, “

Comment on ‘Does the hydrated electron occupy a cavity?

,’”

Science

331

,

1387

(

2011

).

https://doi.org/10.1126/science.1198191

Google Scholar

Crossref

PubMed

119.

L.

Turi

and

A.

Madarász

, “

Comment on ‘Does the hydrated electron occupy a cavity?

,’”

Science

331

,

1387

(

2011

).

https://doi.org/10.1126/science.1197559

Google Scholar

Crossref

PubMed

120.

R. E.

Larsen

,

W. J.

Glover

, and

B. J.

Schwartz

, “

Response to comment on ‘Does the hydrated electron occupy a cavity?

,’”

Science

331

,

1387

(

2011

).

https://doi.org/10.1126/science.1197884

Google Scholar

Crossref

121.

J. R.

Casey

,

R. E.

Larsen

, and

B. J.

Schwartz

, “

Resonance Raman and temperature-dependent electronic absorption spectra of cavity and noncavity models of the hydrated electron

,”

Proc. Natl. Acad. Sci. U. S. A.

110

,

2712

(

2013

).

https://doi.org/10.1073/pnas.1219438110

Google Scholar

Crossref

PubMed

122.

J. R.

Casey

,

A.

Kahros

, and

B. J.

Schwartz

, “

To be or not to be in a cavity: The hydrated electron dilemma

,”

J. Phys. Chem. B

117

,

14173

(

2013

).

https://doi.org/10.1021/jp407912k

Google Scholar

Crossref

PubMed

123.

C.-C.

Zho

and

B. J.

Schwartz

, “

Time-resolved photoelectron spectroscopy of the hydrated electron: Comparing cavity and noncavity models to experiment

,”

J. Phys. Chem. B

120

,

12604

(

2016

).

https://doi.org/10.1021/acs.jpcb.6b07852

Google Scholar

Crossref

PubMed

124.

C.-C.

Zho

,

E. P.

Farr

,

W. J.

Glover

, and

B. J.

Schwartz

, “

Temperature dependence of the hydrated electron’s excited-state relaxation. I. Simulation predictions of resonance Raman and pump-probe transient absorption spectra of cavity and non-cavity models

,”

J. Chem. Phys.

147

,

074503

(

2017

).

https://doi.org/10.1063/1.4985905

Google Scholar

Crossref

PubMed

125.

L. S.

Cederbaum

,

W.

Domcke

,

J.

Schirmer

, and

W.

von Niessen

, “

Correlation effects in the ionization of molecules: Breakdown of the molecular orbital picture

,”

Adv. Chem. Phys.

65

,

115

(

1986

).

https://doi.org/10.1002/9780470142899.ch3

Google Scholar

Crossref

126.

D.

Danovich

, “

Green’s function methods for calculating ionization potentials, electron affinities, and excitation energies

,”

Wiley Interdiscip. Rev.: Comput. Mol. Sci.

1

,

377

(

2011

).

https://doi.org/10.1002/wcms.38

Google Scholar

Crossref

127.

A. P.

Gaiduk

,

M.

Govoni

,

R.

Seidel

,

J. H.

Skone

,

B.

Winter

, and

G.

Galli

, “

Photoelectron spectra of aqueous solutions from first principles

,”

J. Am. Chem. Soc.

138

,

6912

(

2016

).

https://doi.org/10.1021/jacs.6b00225

Google Scholar

Crossref

PubMed

128.

J. W.

Knight

,

X.

Wang

,

L.

Gallandi

,

O.

Dolgounitcheva

,

X.

Ren

,

J. V.

Ortiz

,

P.

Rinke

,

T.

Körzdöfer

, and

N.

Marom

, “

Accurate ionization potentials and electron affinities of acceptor molecules. III. A benchmark of GW methods

,”

J. Chem. Theory Comput.

12

,

615

(

2016

).

https://doi.org/10.1021/acs.jctc.5b00871

Google Scholar

Crossref

PubMed

129.

O.

Dolgounitcheva

,

M.

Díaz-Tinoco

,

V. G.

Zakrzewski

,

R. M.

Richard

,

N.

Marom

,

C. D.

Sherrill

, and

J. V.

Ortiz

, “

Accurate ionization potentials and electron affinities of acceptor molecules. IV. Electron-propagator methods

,”

J. Chem. Theory Comput.

12

,

627

(

2016

);

https://doi.org/10.1021/acs.jctc.5b00872

Crossref

[PubMed]

Google Scholar

Erratum, 13,

389

–

391

(

2017

).

https://doi.org/10.1021/acs.jctc.6b01180

130.

J. V.

Ortiz

, “

Interpreting bonding and spectra with correlated, one-electron concepts from electron propagator theory

,”

Annu. Rep. Comput. Chem.

13

,

139

(

2017

).

https://doi.org/10.1016/bs.arcc.2017.06.002

Google Scholar

Crossref

131.

A.

Thierbach

,

C.

Neiss

,

L.

Gallandi

,

N.

Marom

,

T.

Körzdörfer

, and

A.

Goerling

, “

Accurate valence ionization energies from Kohn-Sham eigenvalues with the help of potential adjustors

,”

J. Chem. Theory Comput.

13

,

4726

(

2017

).

https://doi.org/10.1021/acs.jctc.7b00490

Google Scholar

Crossref

PubMed

132.

H.

Tjong

and

H.-X.

Zhou

, “

On the dielectric boundary in Poisson–Boltzmann calculations

,”

J. Chem. Theory Comput.

4

,

507

(

2008

).

https://doi.org/10.1021/ct700319x

Google Scholar

Crossref

PubMed

133.

H.-X.

Zhou

,

S.

Qin

, and

H.

Tjong

, “

Modeling protein–protein and protein–nucleic acid interactions: Structure, thermodynamics, and kinetics

,”

Annu. Rep. Comput. Chem.

4

,

67

(

2008

).

https://doi.org/10.1016/s1574-1400(08)00004-2

Google Scholar

Crossref

134.

X.

Pang

and

H.-X.

Zhou

, “

Poisson-Boltzmann calculations: van der Waals or molecular surface?

,”

Commun. Comput. Phys.

13

,

1

(

2013

).

https://doi.org/10.4208/cicp.270711.140911s

Google Scholar

Crossref

PubMed

135.

S.

Decherchi

,

J.

Colmenares

,

C. E.

Catalano

,

M.

Spagnuolo

,

E.

Alexov

, and

W.

Rocchia

, “

Between algorithm and model: Different molecular surface definitions for the Poisson-Boltzmann based electrostatic characterization of biomolecules in solution

,”

Commun. Comput. Phys.

13

,

61

(

2013

).

https://doi.org/10.4208/cicp.050711.111111s

Google Scholar

Crossref

PubMed

136.

A. V.

Onufriev

and

B.

Aguilar

, “

Accuracy of continuum electrostatic calculations based on three common dielectric boundary definitions

,”

J. Theor. Comput. Chem.

13

,

1440006

(

2014

).

https://doi.org/10.1142/s0219633614400069

Google Scholar

Crossref

PubMed

137.

A. P.

Gaiduk

,

T. A.

Pham

,

M.

Govoni

,

F.

Paesani

, and

G.

Galli

, “

Electron affinity of liquid water

,”

Nat. Commun.

9

,

247

(

2018

).

https://doi.org/10.1038/s41467-017-02673-z

Google Scholar

Crossref

PubMed

138.

D. A.

Scherlis

,

J.-L.

Fattebert

,

F.

Gygi

,

M.

Cococcioni

, and

N.

Marzari

, “

A unified electrostatic and cavitation model for first-principles molecular dynamics in solution

,”

J. Chem. Phys.

124

,

074103

(

2006

).

https://doi.org/10.1063/1.2168456

Google Scholar

Crossref

139.

See http://osc.edu/ark:/19495/f5s1ph73 for Ohio Supercomputer Center.

2018

Author(s)

Quantum chemistry in arbitrary dielectric environments: Theory and implementation of nonequilibrium Poisson boundary conditions and application to compute vertical ionization energies at the air/water interface

I. INTRODUCTION

II. NONEQUILIBRIUM POISSON FORMALISM

A. Self-consistent equilibrium solvation

B. State-specific nonequilibrium solvation

C. Dielectric function

1. Bulk environment

2. Interfacial environment

III. NUMERICAL SOLUTION OF POISSON’S EQUATION

A. Finite-difference scheme

B. Multigrid approach

IV. COMPUTATIONAL METHODS

A. Molecular dynamics simulations

B. Continuum solvation models

C. Electronic structure calculations

V. RESULTS

A. VIEs in bulk water

B. VIEs at the liquid/vapor interface

VI. DISCUSSION

VII. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Quantum chemistry in arbitrary dielectric environments: Theory and implementation of nonequilibrium Poisson boundary conditions and application to compute vertical ionization energies at the air/water interface

I. INTRODUCTION

II. NONEQUILIBRIUM POISSON FORMALISM

A. Self-consistent equilibrium solvation

B. State-specific nonequilibrium solvation

C. Dielectric function

1. Bulk environment

2. Interfacial environment

III. NUMERICAL SOLUTION OF POISSON’S EQUATION

A. Finite-difference scheme

B. Multigrid approach

IV. COMPUTATIONAL METHODS

A. Molecular dynamics simulations

B. Continuum solvation models

C. Electronic structure calculations

V. RESULTS

A. VIEs in bulk water

B. VIEs at the liquid/vapor interface

VI. DISCUSSION

VII. CONCLUSIONS

ACKNOWLEDGMENTS

REFERENCES

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only