Impact of solvation on the GW quasiparticle spectra of molecules

,


I. INTRODUCTION
Electrochemical reactions will underpin the decarbonization of the global economy by enabling the widespread use of fuel cells, sustainable production of hydrogen and ammonia, and reduction of CO 2 .][3] Improving the accuracy of first-principles electrochemical predictions requires better treatment of both electronic structure and solvation by the electrochemical environment. 4][7] This shortcoming of DFT is particularly problematic in studies of (photo)electrochemical reactions because accurate prediction of frontier state alignment for (photo)electrocatalysts relative to those of reactants/products is fundamental in the discovery of new catalysts and reaction mechanisms. 2 The incorporation of exact exchange in hybrid functionals partially corrects the SIE 8 but can add empirical parameters that decrease the transferability of such approaches.
0][11][12][13][14][15][16] This is primarily due to the complexity in sampling solvent molecules configurations in explicit solvation approaches, and the lack of a complete theoretical formulation and computational tools for GW calculations with implicit solvation approaches.
Explicit solvation accounts for solvent chemical interactions with the electrode and adsorbed species, but including more than the first solvation shell typically requires expensive molecular dynamics simulations to sample the large thermodynamic phase space of solvent configurations. 17Ab initio molecular dynamics of electrolytes are particularly challenging due to the need for long timescale and large simulation cells to capture statistically significant numbers of ions. 4 These challenges are substantially magnified when combined with GW calculations due to its significantly higher computational expense as compared to semi-local generalizedgradient approximation (GGA) calculations.Nevertheless, this approach has been used to predict the G 0 W 0 @(GGA/hybrid) photoelectron spectra of NaCl + 52 H 2 O, 17 G 0 W 0 @(PBE, r 2 SCAN, PBE0) IPs of solvated transition metal ions, 18 and stochastic G 0 W 0 @PBE photoemission spectra of organic molecules. 19,20mplicit solvation approaches, which replace explicit solvent configurations with polarizable continua or liquid density distributions, are more computationally expedient and widely applied in conjunction with DFT in first-principles electrochemistry.These models bypass the phase-space sampling issue and capture the impact of the liquid in a single electronic structure calculation, which is particularly attractive for combining with high-fidelity electronic structure methods such as GW.Some of the ingredients for solvated GW emerged in explorations of electrochemical systems with GW using solvated-DFT starting points, 21 and in dielectric screening models to account for substrate effects in GW calculations of low-dimensional materials. 226][27][28][29] However, the impact of dynamic screening by the solvent in GW calculations has not yet been explored, and solvated GW is not yet available in the widely-applied first-principles electrochemistry tools.
Here, we fully combine the GW and implicit solvation approaches by incorporating a full-frequency description of the solvent polarizability into the screened coulomb potential.This approach is generalizable to any solvation model, and applicable to systems of any dimensionality and charge state described with the joint DFT approach. 30We demonstrate this approach using the nonlocal solvation model, the spherically-averaged liquid susceptibility ansatz (SaLSA), 31 which includes no empirical parameters in the screening response, accurately describes neutral and charged surfaces in solution, and as we show here, facilitates a convenient representation of the frequency-dependent solvent polarizability for GW calculations.We show using GW calculations that solvent screening strongly affects the frontier orbital energies in small molecules due to the complicated nature of their self-energies, 32 and thereby modifies the HOMO-LUMO gaps of solvated molecules by several eV.Finally, we implement these methods by interfacing the widelyused open-source software packages for GW and solvated DFT, BerkeleyGW 33 and JDFTx 34 respectively, in order to facilitate widespread application of solvated beyond-DFT methods for electrochemistry.

II. THEORY
The GW approach solves the Dyson equation for the quasiparticle (QP) energies and wavefunctions (using Hartree atomic units), where ∇ is the kinetic energy operator, V ion is the ionic potential, V H is the Hartree potential, Σ is the nonlocal, energy-dependent self energy within the GW approximation, and E QP nk and ψ QP nk are the quasiparticle energies and wavefunctions, respectively.We use the singleparticle wavefunctions from DFT as our quasiparticle wavefunctions (the diagonal approximation) and iterate to self-consistency in the quasiparticle energies that enter into the Green's function G in the GW approximation.GW calculations carried out using the DFT vacuum or SaLSA wavefunctions are labeled as GW@GGA or GW@GGA+SaLSA, respectively.
The self-energy operator Σ depends on the screened Coulomb interaction, where v is the bare Coulomb interaction, ε is the fullfrequency dielectric matrix.Within the random phase approximation [35][36][37] ε is given by Above, χ 0 is the independent-electron polarizability, written in the plane-wave basis as ) with the matrix elements M defined by Here, q is a vector in the first Brillouin zone, G and G ′ are reciprocal-lattice vectors, ω is frequency, δ is an infinitesimal that locates the poles of χ 0 relative to the real frequency axis, ψ nk and ϵ nk are the Kohn-Sham orbitals and energies from the DFT starting point, and n and n ′ count occupied and empty bands, respectively.See Refs.[38-46] for a detailed review of the GW method.We focus below on the modifications for solvation.
To account for solvent screening, we replace χ 0 in Eq. 3 with χ solv = χ 0 + χ lq , where χ lq is the independent-particle susceptibility of the liquid using an implicit solvation model.This approach is labeled GW@GGA+SaLSA+χ lq below.Note that χ lq is also a response to the total Coulomb potential, analogous to χ 0 for the electrons, which is necessary for this additivity.Most implicit solvation models incorporate a local dielectric response with a spatially-dependent relative permittivity ϵ lq (r), leading to χlq being a differential operator defined by χlq ϕ(r) = 1 4π ∇ • (ϵ lq (r) − 1)∇ϕ(r), where ϕ(r) is the electrostatic potential.This χ lq which can be directly combined with χ 0 to account for solvent screening.However, for these local models χ lq GG ′ ∝ G•G ′ in a planewave basis, which then requires a large plane-wave cutoff to properly account for the solvation effects in the dielectric matrices and thus complicate convergence of the GW calculation.
Here, we instead use the nonlocal susceptibility from the SaLSA solvation model, 31 The first term captures the rotational response of solvent molecules with density N 0 (r) ( Ñ0 is its Fourier transform) at temperature T , with ρ lm mol being the charge density distribution of a solvent molecule projected to each spherical harmonic component indexed by l, m.The second term captures the electronic polarization response of solvent molecules from atomic sites with densities N 0 α (⃗ r), each with susceptibilities of strength χ α and non-local weight function w α (r).All the above parameters are extracted from DFT calculations of solvent molecules, while the correlation strength functions C rot (ω) and C pol (ω) are constrained by the experimental rotational and electronic dielectric response of the solvent.See Ref. 31  for the derivation and complete parameterization of the SaLSA solvation model.
Most importantly, the nonlocality of the SaLSA model leads to χ lq GG ′ → 0 as G → ∞ because the nonlocal weight functions ρlm mol (G) and wα (G) → 0 as G → ∞ (Figure S1a).This allows the liquid susceptibility to be converged with a small plane wave cutoff, in comparison to the χ lq GG ′ ∝ G • G ′ of local models (Figure S1b), that require additional numerical care in constructing and inverting ε GG ′ .We implemented calculation of χ lq GG ′ (q; ω) in the plane-wave basis in the open-source JDFTx software, 34 exported directly as a matrix in the plane-wave basis for each q and ω, which can then be included in the construction and inversion of ε GG ′ (q; ω) already implemented within the open-source BerkeleyGW software using the static subspace approximation. 33,47lthough the calculations in this work only used a single k/q-point, we emphasize that this formalism is general and can be used to predict electronic structure properties of solvated surfaces as well.
Finally, note that the rotational response contributes only at frequencies that are orders of magnitude lower than the electronic frequencies, and can therefore be dropped from Eq. 6 in calculating the screening for GW.However, the rotational response dominates in the static limit (ω → 0) and must be retained in the solvated DFT calculation used as the starting point.

III. COMPUTATIONAL METHODS
All DFT calculations were performed using the JDFTx software package 34 using the Perdew-Burke-Ernzerhof (PBE) GGA functional 48 , 60 Ryd planewave cutoff, and the SG15 ONCV norm-conserving pseudopotentials. 49mplicit solvation effects were included using the SaLSA model. 31The Coulomb interaction between periodic cells was numerically truncated at both the DFT and GW levels. 50All simulation boxes measured 10 Å×10 Å×10 Å and the molecules were relaxed using the PBE functional.The full-frequency GW calculations were run in Berke-leyGW using the static subspace approximation and the static remainder approach to correct for the slow convergence of the self-energy with number of empty states. 9,33he polarizability used for GW calculations was com-puted with a cutoff of 25 Ryd, and 81 real frequencies spanning 20 eV and 15 imaginary frequencies.Real and imaginary frequencies are required for evaluation of the self energy using contour deformation. 51,52Increasing the number of imaginary frequencies to 25 lowered the CO GW 0 HOMO-LUMO gap by 8 meV.The total number of unoccupied bands in the calculation was approximately 14000 and reduced to approximately 500-800 bands using a stochastic pseudobands approach for bands at least 3 Ryd above the LUMO using a 2% accumulation window. 53,54The quasiparticle eigenvalues were converged to within 10 meV by performing two additional iterations in G following the initial G 0 W 0 calculation.We found that increasing the planewave and screened cutoffs to 100 and 30 Ryd, respectively, increased the predicted CO GW 0 HOMO-LUMO gap by 24 and 18 meV, respectively.The change in CO and H 2 O GW 0 HOMO-LUMO gap calculated using experimental structure parameters from those calculated using their PBE-relaxed structures were -16 and 37 meV, respectively.

IV. RESULTS AND DISCUSSION
We now discuss how molecular HOMO-LUMO gaps and frontier orbital energies change with computational methodology.Summaries of all of the frontier orbital and gap energies discussed in this section are provided in Tables SI, SII, and SIII.Additional eigenvalues are provided in Table SIV.As expected, the HOMO-LUMO gaps of the molecules in vacuum significantly increases (median increase of 6.36 eV) from the GGA-predicted gaps (Figure 1).For all molecules except F 2 and H 3 O + , the majority of the gap opening occurs due to a much larger decrease in the HOMO energy (median decrease of 4.96 eV) than the increase in the LUMO energy (median increase of 1.74 eV) (Figure 2).Our predicted vacuum GW 0 @GGA gaps are in good agreement with experiment (Table SI) and previous studies using similar methodologies. 5,10he inclusion of implicit solvation does not significantly affect the HOMO-LUMO gaps nor frontier orbital placements of the molecules at the DFT level (Figure 2).This is expected because the contribution to the electrostatic potential from the solvation model is mostly constant / varies slowly in the spatial region of the solute, and therefore contributes an almost rigid shift in the DFT orbital energies of molecules.Similarly, implicit solvation incorporated at the GW level only via modification of the DFT wavefunctions (GW@GGA+SaLSA), produces HOMO/LUMO gaps and energies that are similar to the vacuum GW 0 results.
In contrast, however, augmenting χ 0 with the frequency-dependent solvent polarizability significantly changes the predicted G 0 W 0 and GW 0 eigenvalues.The GW 0 @GGA+SaLSA+χ lq method significantly lowers the GGA+SaLSA HOMO levels (median decrease of 4.30 eV).Unlike in vacuum or when using only the SaLSA FIG. 1. Summary of predicted HOMO-LUMO gaps of small molecules using GGA in vacuum (black), GGA plus SaLSA implicit solvation (red), GW0 with GGA wavefunctions in vacuum (gold), GW0 with GGA plus SaLSA implicit solvation wavefunctions (blue), and GW0 with GGA plus SaLSA implicit solvation wavefunctions and liquid polarizability (purple).The dotted lines are included as guides to the eye.
wavefunctions however, the GW 0 @GGA+SaLSA+χ lq method LUMO energies are typically lowered from their GGA energies (median decrease of 1.36 eV) rather than being increased.Consequently, the GW 0 @GGA+SaLSA LUMO energy for all molecules is negative instead of positive like for the GW 0 @GGA and GW 0 @GGA+SaLSA methods.The combined effects of χ lq on the molecules' HOMO and LUMO energies only partially offset, leading to significantly reduced GW 0 @GGA+SaLSA+χ lq HOMO-LUMO gaps as compared to the other GW 0 methods (Figure 1).Overall, when compared to the GW 0 @GGA+SaLSA method, the median effect of the GW 0 @GGA+SaLSA+χ lq method increased the HOMO energy by 0.71 eV, decreased the LUMO energy by 3.04 eV, and decreased the HOMO-LUMO gap by 3.70 eV.
Solvated chemical systems with complicated inverse dielectric functions and self-energies across the frequency axis are expected to most benefit from the approach in this work. 10,32The impact of using the frequencydependent polarizability on the G = G ′ = 0 component of the real inverse dielectric for H 2 O and benzene is shown in Figure 3.We observed that the low frequency region of Re(ε −1 00 (ω)) was not strongly affected by χ lq .Instead, the increased screening due to the solvent typically adjusted the higher frequency components of Re(ε −1 00 (ω)) to be closer to 1.
Although the primary goal of this work is the devel- opment of a general formalism needed that integrates quasiparticle calculations with implicit solvation model screening, comparisons to experimental molecular IPs using photoelectron spectroscopy (PES) can be made to validate the accuracy of this approach for occupied states.6][57][58][59][60][61] Additionally, other considerations such as PES peak broadening due to solvent and  I summarizes the available experimental IPs and the corresponding values predicted by the current approach.The results show that this methodology predicts reasonably accurate solvated IPs given existing measurement uncertainty.This suggests the SaLSA implicit solvation model may capture the ensemble solvent interaction with solute HOMOs, potentially reducing the need for full explicit solvation and requiring a less extensive degree of statistical sampling.Similarly, although the band gap of liquid water has been the subject of several investigations, no clear consensus exists on the exact gap.3][64][65][66][67][68][69] Our predicted H 2 O GW 0 @GGA+SaLSA+χ lq HOMO-LUMO gap of 8.669 eV is similar to these values.The accuracy of our predicted solvated IPs is consistent with the observation that the GW approximation performs surprisingly well in predicting accurate molecular IPs in vacuum given the relative simplicity of the method. 10,70We also note that because our approach considers off-diagonal and frequency-dependent elements in the solvent polarizability (Figure S2), it is consistent with the well-known finding that the GW approximation is most accurate when both local field and dynamical effects are included. 9espite the accuracy of the predicted IPs, it is important to note that the predicted GW 0 @GGA+SaLSA+χ lq LUMO values all lie between approximately -2 and -3 eV, when the solvated water LUMO has been previously predicted to be -0.8 -0.1 eV by beyond-DFT computational studies using fully explicit water models. 62,66,67This discrepancy could originate from one or more sources, including 1) the PBE starting point for the unoccupied states may be particularly inaccurate, 2) unoccupied states are significantly influenced by interactions with the explicit water molecules included in previous studies, and/or 3) the SaLSA implicit solvation model is known to poorly describe the asymmetry in solvation of cations and anions 71 and thus could be expected to less accurately predict an N+1 particle property such as the EA (=ϵ LU M O ).Future work elucidating the impact of nearby explicit water molecules, functional starting point, and improvements in implicit solvation models are currently ongoing.Nevertheless, we believe that the approach described here will help facilitate beyond-DFT modeling of solvated molecular and electrocatalytic systems within the joint-DFT framework and only continue to improve as advancements in solvation models and DFT functionals are made.

V. CONCLUSIONS
In this work, we demonstrated that the electronic polarizability in full-frequency GW calculations can be augmented by the frequency-dependent fluid polarizability from the SaLSA implicit solvation model to capture the effects of solvent screening in beyond-DFT electronic structure calculations without the need for sampling of explicit water structures.This approach is generalizable to any solvation model and dimensionality of chemical system, and has been tested extensively here for molecular systems that are expected to be most susceptible to the liquid screening effects.The solvent polarizability significantly alters the frontier orbital positions and in particular lowers the energy of all LUMO energies.The change in frontier orbital energies leads to a net decrease in the molecular HOMO-LUMO gaps by 3 -5 eV on average.We showed that solvent screening is central to these electronic structure changes, and that wavefunctions generated from solvated DFT calculations alone do not result in GW frontier orbital energies much different than DFT wavefunctions generated in vacuum.The technique demonstrated here can be applied to quasiparticle energies of electrochemical interfaces in future work, laying the groundwork for beyond-DFT first-principles electrochemistry.

FIG. 2 .FIG. 3 .
FIG.2.Summary of predicted HOMO (filled points) and LUMO (unfilled points) levels of small molecules using GGA in vacuum (black), GGA plus SaLSA implicit solvation (red), GW0 with GGA wavefunctions in vacuum (gold), GW0 with GGA plus SaLSA implicit solvation wavefunctions (blue), and GW0 with GGA plus SaLSA implicit solvation wavefunctions and liquid polarizability (purple).The dotted lines are included as guides to the eye.

TABLE I .
[56][57][58]f GW0@GGA+SaLSA+χ lq IPs with experimentally measured solvated IPs using photoelectron spectroscopy.The experimental IPs are measured using the peak maximums in the PES spectra.The 2 nd IPs of H2O are the average of the 3a1H and 3a1L values.[56][57][58]Toourknowledge, the only molecules in the current test set with existing experimental solvated IP measurements are H 2 O, H 3 O + , and phenol.Finally, because the systems studied in this work consist of only a single explicit solute molecule, we compared the GW 0 @GGA+SaLSA+χ lq predicted IPs to the available PES spectra peak energies rather than the linear extrapolations of the threshold values.Table