Impact of solvation on the GW quasiparticle spectra of molecules

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$⁠. First-principles density functional theory (DFT) modeling of electrochemical reactions has proven invaluable in the discovery of new catalysts for electrochemical and photochemical energy conversion, as well as battery electrode materials.^1–3 

Improving the accuracy of first-principles electrochemical predictions requires better treatment of both electronic structure and solvation by the electrochemical environment.⁴ On the electronic structure side, semi-local DFT methods have notorious limitations in correctly predicting the localization of charge, band/HOMO–LUMO gaps, and alignment of frontier states due to self-interaction error (SIE).^5–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.² The incorporation of exact exchange in hybrid functionals partially corrects the SIE⁸ but can add empirical parameters that decrease the transferability of such approaches.

Many-body perturbation theory calculations using the GW approximation starting from DFT wavefunctions (GW@DFT) are thus appealing due to their fully ab initio approach and well-known accuracy in predicting one-particle properties, including ionization potential (IP), electron affinity (EA), and energy gaps in a wide variety of systems.^5,9–16 Calculations using the GW approximation for periodic materials have typically been applied toward the study of bulk solids or molecules/surfaces in vacuum calculations.^5,9–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.^17, Ab 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.⁴ These challenges are substantially magnified when combined with GW calculations due to its significantly higher computational expense as compared to semi-local generalized-gradient 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,¹⁷ G $0$W $0$@(PBE, r $2$SCAN, PBE0) IPs of solvated transition metal ions,¹⁸ and stochastic G $0$W $0$@PBE photoemission spectra of organic molecules.^19,20

Implicit 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²¹ and in dielectric screening models to account for substrate effects in GW calculations of low-dimensional materials.²² Recent studies have integrated solvation models in the DFT starting point and as a dielectric profile modifying either the static dielectric matrix,^23,24 Coulomb potential, or static COHSEX self-energy.^25–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.³⁰ We demonstrate this approach using the nonlocal solvation model, the spherically averaged liquid susceptibility ansatz (SaLSA),³¹ 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³² and thereby modifies the HOMO–LUMO gaps of solvated molecules by several eV. Finally, we implement these methods by interfacing the widely used open-source software packages for GW and solvated DFT, BerkeleyGW,³³ and JDFTx,³⁴ respectively, in order to facilitate widespread application of solvated beyond-DFT methods for electrochemistry.

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. When inverting the dielectric matrix with $χ solv$⁠, we are including the solute and solvent degrees of freedom in W. 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$ can be directly combined with $χ 0$ to account for solvent screening. However, for these local models, $χ lq G G ′∝ G⋅ G ′$ in a plane-wave 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.

See Ref. 31 for the derivation and complete parameterization of the SaLSA solvation model. Briefly, SaLSA is derived from the linear response limit of a polarizable classical density functional approximation of the nonlocal and nonlinear response of the solvent in general.⁴⁷ The density of the molecule centers, $N 0( r →)=N bulks( r →)$⁠, is set based on the bulk number density of molecules in the fluid, $N bulk$⁠, and a cavity shape function $s( r →)$ that switches from 0 in the solute region to 1 in the solvent region, with the orientation probability assumed to be uniform in the absence of electrostatic interactions with the solute. SaLSA is derived as the first-order perturbation in the rotational and polarization response with respect to this zeroth-order starting point,³¹ which leads to an expansion in multipole components of the solvent charge density (⁠ $ρ mol l m$⁠) for the rotational component and a dipole contribution alone for the polarization component.

Most importantly, the nonlocality of the SaLSA model leads to $χ lq G G ′→0$ as $G→∞$ because the nonlocal weight functions $ρ ~ mol l m( G)$ and $w ~ α( G)→0$ as $G→∞$ [Fig. S2(a)]. This allows the liquid susceptibility to be converged with a small plane wave cutoff, in comparison to the $χ lq G G ′∝ G⋅ G ′$ of local models [Fig. S2(b)], that require additional numerical care in constructing and inverting $ϵ G G ′$⁠. We implemented calculation of $χ lq G G ′( q;ω)$ in the plane-wave basis in the open-source JDFTx software,³⁴ 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 $ϵ G G ′( q;ω)$ already implemented within the open-source BerkeleyGW software using the static subspace approximation.^33,50 Although 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.

All DFT calculations were performed using the JDFTx software package³⁴ using the Perdew–Burke–Ernzerhof (PBE) GGA functional,⁵¹ 60 Ryd planewave cutoff, and the SG15 ONCV norm-conserving pseudopotentials.⁵² Implicit solvation effects were included using the SaLSA model.³¹ The Coulomb interaction between periodic cells was numerically truncated at both the DFT and GW levels.⁵³ All simulation boxes measured $10×10×10 Å 3$ and the molecules were relaxed using the PBE functional. The full-frequency GW calculations were run in BerkeleyGW 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,33 The polarizability used for GW calculations was computed 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.^54,55 Increasing 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 14 000 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.^56,57 The 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 was $−$16 and 37 meV, respectively.

We now discuss how molecular HOMO–LUMO gaps and frontier orbital energies change with computational methodology. Summaries of all of the gap and frontier orbital energies discussed in this section are provided in Figs. 1, 2, S4, and S5 and Tables SI, SII, and SIII. We note that depending on method, CO, F $2$⁠, N $2$⁠, benzene, phenol, thiophene, and 1,2,5-thiadiazole can have an orbital predicted to be the LUMO with GW that was not the LUMO in the initial GGA calculation. Because we focus here the impact of this formalism on GW-predicted frontier states, all GGA and GGA+SaLSA results shown in Figs. 1 and 2 used the same orbital indices predicted to be the HOMO/LUMO by the GW $0$@GGA and GW $0$@GGA+SaLSA+ $χ lq$ methods, respectively. For similar reasons, all summary statistics were calculated using this data set. Alternative plots using the HOMO/LUMO levels predicted from each method individually are shown in Figs. S4 and S5. Additional eigenvalues outside of the frontier orbitals are provided in Table SIV.

As expected, the GW $0$@GGA HOMO–LUMO gaps of the molecules in vacuum significantly increase (median increase of 6.03 eV) from the GGA-predicted gaps (Fig. 1). For all molecules except 1,2,5-thiadiazole and F $2$⁠, 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.30 eV) (Fig. 2). Our predicted vacuum GW₀@GGA gaps are in good agreement with experiment (Table SI) and previous studies using similar methodologies.^5,10

The inclusion of implicit solvation does not significantly affect the HOMO–LUMO gaps nor frontier orbital placements of the molecules at the DFT level (Fig. S5). 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. We note that because SaLSA is an implicit solvation model that primarily approximates the electrostatic interaction between the solute and the solvent, it does not typically introduce state reordering or affect the state degeneracies. 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 wavefunctions however, the GW $0$@GGA+SaLSA+ $χ lq$ method LUMO energies are typically lowered from their GGA energies (median decrease of 2.30 eV) rather than being increased. Consequently, the GW $0$@GGA+SaLSA+ $χ lq$ LUMO energy for each molecule is negative instead of typically being 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 (Fig. 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,32 The impact of using the frequency-dependent solvent polarizability on the $G= G ′=0$ component of the real inverse dielectric for H₂O and benzene is shown in Fig. 3. We observed that the low frequency region of Re(⁠ $ϵ 00 − 1$(⁠ $ω$⁠)) was not strongly affected by $χ$^lq. Instead, the increased screening due to the solvent typically adjusted the higher frequency components of Re(⁠ $ϵ 00 − 1$(⁠ $ω$⁠)) to be closer to 1.

Although the primary goal of this work is the development 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. Measurements of experimental solvated molecular IPs are non-trivial, however, and the development of advanced techniques that increase accuracy is still an ongoing area of research.^58–64 Additionally, other considerations such as PES peak broadening due to solvent and the extrapolation procedure used to obtain PES signal onset energies can further obscure direct experimental and computational comparison.^61,65 To our knowledge, 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₀@GGA+SaLSA+ $χ$^lq predicted IPs to the available PES spectra peak energies rather than the linear extrapolations of the threshold values.

Table 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 bandgap of liquid water has been the subject of several investigations, no clear consensus exists on the exact gap. Briefly, a commonly cited experimental gap estimate is $8.7±0.5$ eV with a more recent estimate of $9.0±0.2$ eV, while the computational gap has ranged from 7.8 to 10.5 eV.^65–71 Our predicted H₂O GW₀@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,72 We also note that because our approach considers off-diagonal and frequency-dependent elements in the solvent polarizability (Fig. S3), it is consistent with the well-known finding that the GW approximation is most accurate when both local field and dynamical effects are included.⁹ 

Despite the accuracy of the predicted IPs, it is important to note that the predicted GW₀@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 to $−$0.1 eV by beyond-DFT computational studies using fully explicit water models.^65,68,69 This 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, (3) the more diffuse LUMO states are more affected by the environment screening than occupied states, which are expected to be more affected by screening from the molecule itself, and/or (4) the SaLSA implicit solvation model is known to poorly describe the asymmetry in solvation of cations and anions⁷³ and thus could be expected to less accurately predict an N+1 particle property such as the EA (⁠ $=− ε L U M O$⁠). Future work elucidating the impact of nearby explicit water molecules, functional starting point on the bound charge and wavefunction localization, and improvements in implicit solvation models is 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.

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 LUMOs. 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.

See the supplementary material for figures showing the frequency-dependent liquid polarizability parameterization, data contained within the liquid polarizability matrix, alternative HOMO–LUMO gap and edge plots, three summary tables of the vacuum and solvated (with and without χ lq) results, and a full table of DFT and QP eigenvalues calculated with the different methods.

This work was supported by the U.S. Department of Energy, Office of Science, Basic Energy Sciences, under Award No. DE-SC0022247. We used computational resources sponsored by the Department of Energy’s Office of Energy Efficiency and Renewable Energy, located at the National Renewable Energy Laboratory, and at the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility located at Lawrence Berkeley National Laboratory, operated under Contract No. DE-AC02-05CH11231 using NERSC Award No. ERCAP0020105.

The authors have no conflicts to disclose.

Jacob M. Clary: Conceptualization (supporting); Data curation (lead); Formal analysis (lead); Funding acquisition (supporting); Investigation (lead); Software (supporting); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (equal). Mauro Del Ben: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal). Ravishankar Sundararaman: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (equal); Visualization (equal); Writing – review & editing (equal). Derek Vigil-Fowler: Conceptualization (equal); Funding acquisition (equal); Methodology (equal); Project administration (equal); Software (equal); Supervision (lead); Visualization (equal); Writing – review & editing (equal).

The data that support the findings of this study are available from the corresponding authors upon reasonable request.