Cooperation and Environment Characterize the Low-Lying Optical

molecule, the joint density of states of the liquid is much broader than that of the ..... a velocity rescaling thermostat) in order to effectively ac...
0 downloads 0 Views 743KB Size
Download PDF
Letter pubs.acs.org/JPCL

Cooperation and Environment Characterize the Low-Lying Optical Spectrum of Liquid Water Sudheer Kumar P., Alessandro Genova, and Michele Pavanello* Department of Chemistry, Rutgers University, Newark, New Jersey 07102, United States S Supporting Information *

ABSTRACT: The optical spectrum of liquid water is analyzed by subsystem timedependent density functional theory. We provide simple explanations for several important (and so far elusive) features. Due to the disordered environment surrounding each water molecule, the joint density of states of the liquid is much broader than that of the vapor, thus explaining the red-shifted Urbach tail of the liquid compared to the gas phase. Confinement effects provided by the first solvation shell are responsible for the blue shift of the first absorption peak compared to the vapor. In addition, we also characterize many-body excitonic effects. These dramatically affect the spectral weights at low frequencies, contributing to the refractive index by a small but significant amount.

W

the latter, many-body refers to electron−hole interactions. We set out to investigate how many-body interactions affect the optical spectrum and other related quantities (such as the refractive index). These can be cooperative or anticooperative in nature. Many-body effects have been discussed before in terms of the screening properties of the bulk in the computation of selfenergies for GW calculations.5,6 It was found that screening is independent of the particular configuration of water considered. Thus, it can be inferred that it is not affected by the structure of the environment surrounding the water molecules. It has also been shown that for ice cooperative many-body effects in the form of excitonic couplings increase the oscillator strength of low-lying excitations and are responsible for an increase of the index of refraction with increasing pressure.7 We also set out to investigate coupling between the first absorption band and the structure of the liquid. This helps us understand what influences the peak position and shape (broadening). The underlying reasons for a blue shift of the first absorption band of water from 7.4 eV in the gas phase8 to 8.3−8.6 eV in the liquid phase9−11 are still a matter of debate. G0W0 calculations based on PBE orbitals predict a shift of about 0.7 eV and attribute it to excitonic effects6,12 (i.e., the isolated water molecule features a more localized exciton than that of the liquid phase). Calculations based on clusters, instead, define the blue shift as a result of electrostatic embedding13 (a sort of confinement effect). This view is corroborated by excitonically coupled coupled cluster (CCSD) and semilocal time-dependent density functional theory

ater is the most important liquid on Earth. It plays roles in most biological transformations and in the properties of materials. For example, material’s interfaces with water are ubiquitous. Thus, the ability to predict the response of aqueous solutions to external perturbations (such as electromagnetic fields) constitutes a first important step in the predictive modeling of the majority of life and energy processes. Particularly important is the ability to determine the local response of water for both the electrons1,2 and the vibrations.3,4 This is because the so-called “specific interactions” of water molecules with other entities often determine their optical properties. In this respect, formally decomposing the response into local contributions can be a powerful tool for describing how water behaves when it interacts with localized external perturbations. Although there has been tremendous progress in the ab initio modeling of optical properties of water, difficulties arise because the disordered nature of liquid water requires large simulation cells. In turn, this forces the use of either approximate methods on large cells or accurate methods on cells that are too small. As a result of this limitation, there are several open questions and interesting features of the optical absorption spectrum of water that are yet to be fully explained. In this work, we explore two themes: (1) many-body excitonic interactions between the water molecules and (2) coupling of the first absorption band to the nuclear degrees of freedom describing the liquid structure. By many-body effects, we mean the effects that arise when single water molecules (single bodies) in the liquid interact with each other (other bodies) both in the ground state as well as in their excited states. Although related, this definition is different in spirit from the kinds of many-body effects that a Bethe−Salpeter equation (BSE) treatment would recover. In © XXXX American Chemical Society

Received: August 21, 2017 Accepted: October 2, 2017 Published: October 2, 2017 5077

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083

Letter

The Journal of Physical Chemistry Letters (TDDFT) calculations showing that the first absorption band is composed of localized states.14,15 A recent analysis of the nonself-consistency of the previously mentioned G0W0 results uncovered some possible biases that are carried over from using DFT orbitals in the GW+BSE part of the calculation,16 putting into question the previously predicted exciton sizes of liquid and vapor. We set out to attack these open questions with simulations based on subsystem DFT.17−19 The electron density is decomposed into subsystem contributions, ρ(r) = ∑NI S ρI(r), where Ns is the number of subsystems. The subsystem densities are recovered variationally solving subsystem-coupled Kohn−

(

1

propagator, totaling 50 fs of electron dynamics carried out for each of the three Cartesian directions and for each subsystem. We computed the optical spectrum of 10 snapshots sampled from subsystem DFT-based Born−Oppenheimer (BO) dynamics of bulk liquid water modeled by 64 independent water molecules in a cubic box (a = 12.4278 Å, presented in previous work23 and chosen to reproduce the experimental density of liquid water at ambient conditions). The snapshots were chosen making sure that (1) the structures of the snapshots are not correlated with each other and (2) the average and variance of the number of hydrogen bonds per water molecule are similar to the same values computed along the entire BO trajectory. We provide further details about the snapshot selection in the Supporting Information. To summarize ref 23, the dynamics featured 81000 steps with a time step of 26 au. The first 12000 steps were discarded, and the remaining 69000 steps were used for snapshot selection. The simulation temperature was chosen to be slightly elevated (340 K, imposed by a velocity rescaling thermostat) in order to effectively account for nuclear quantum effects. A more accurate assessment of nuclear quantum effects on the optical properties of water have been recently discussed by others.31 The equations of motion were propagated by the Verlet algorithm. Although it has been shown that dispersion corrections are important to quantitatively reproduce the structure of liquid water,27 in ref 23 we indicated that even without including dispersion corrections the structure of water as computed by eQE is in satisfactory agreement with current experimental characterizations of the liquid. In essence, error cancellation between the xc and nonadditive Ts functionals tilt the interactions toward a good description of hydrogen bonds and, overall, a good structure and dynamics of the liquid. We employed the LC9432 nonadditive Ts functional, which was shown to satisfactorily reproduce energy surfaces of CCSD(T) benchmarks for the water dimer as well as the structural parameters of the liquid (e.g., radial and angular distribution functions).23 The approach implemented in eQE allows us to use a subsystem-specific auxiliary periodic cell exclusively to represent the subsystem Kohn−Sham Hamiltonian and eigenfunctions in a reduced plane wave basis. The Hartree potential is still computed in the large (physical) cell, and thus, the intrinsic periodicity of the system is well represented. The subsystem auxiliary cells were chosen to be 60% of the length of the physical cell. This approach led us to use 1.9 million plane waves less in the representation of the Hamiltonian and waves. Note: the number of plane waves generated in a cell with a given cutoff is proportional to the cube of the lattice constant. Thus, this procedure generates massive computational savings.24 The real-time subsystem dipole change, δμI(t) = ∫ r (ρ(r, t) − ρ(r,t = 0)) dr, is Fourier transformed to frequency space to yield the isotropic dipole polarizability, αI(ω), which is related to the frequency-dependent dielectric constant,11,33 Im[ε(ω)−1] ∝ Im[⟨α(ω)⟩], where the average is carried out over all water molecules for all snapshots,

)

I (r) ϕi I(r) = εiIϕi I(r) , Sham equations, − 2 ∇2 + νsI(r) + νemb

where νIs is the Kohn−Sham potential of the isolated subsystem evaluated at ρI(r). The embedding potential, νIemb, contains exact Coulomb interactions with surrounding subsystems, as well as functional derivatives of nonadditive exchange− correlation (xc) and a nonadditive noninteracting kinetic energy functional (Ts).17−19 To access information about excited states, subsystem TDDFT can be formulated either in the frequency20 or time domain.19 Subsystem TDDFT allows us to approach larger supercells than before and naturally gives us access to many-body effects without compromising the accuracy of the model. Many-body effects arise as intersubsystem static and time-dependent interactions.19,21,22 These favorable qualities are further complemented by the ability of subsystem DFT to naturally localize the electronic structure of the subsystems while still allowing electron density overlap. We have shown, for example, that it can quantitatively reproduce the dynamics and structure of liquid water because the too strong hydrogen bond resulting from self-interaction in semilocal xc functionals is offset by the subsystem-local electronic structure.23 While we do not expect that subsystem electron density localization be a generally beneficial feature of the model, it is beneficial in the limit of simulating systems composed of noncovalently bound subunits. Such as molecular liquids, crystals, and layered periodic systems.24 Simulations were carried out with the embedded QuantumESPRESSO (eQE) package,24 employing the PBE xc functional and Ultrasoft pseudopotentials (O.pbe-rrkjus.UPF/H.pberrkjus.UPF from the Quantum-ESPRESSO PP Library25). We employed plane wave kinetic energy cutoffs of 50 and 500 Ry for the waves and charge densities, respectively. While a 50 Ry cutoff for the expansion in plane waves of the pseudowaves may seem low, we stress that we are employing Ultrasoft pseudopotentials that, at the expense of relaxing norm conservation (recovered using an appropriate generalized orthonormality condition), allow the use of a much reduced plane wave cutoff while maintaining a quantitative description of the electronic structure of the system.26 This is in contrast to norm-conserving pseudopotentials, which are known to require a cutoff much higher than 50 Ry for water,27 especially when wishing to approach the basis set completeness limit.28,29 In the Supporting Information, we further dwell on our choice of energy cutoff and present additional figures featuring additive and interaction energy vs cutoff. Real-time subsystem TDDFT was implemented in the linear response regime, applying a δ electric field perturbation30 of 0.02 Ry/Å, a time step of 1 as, and assuming the adiabatic approximation for all of the density-dependent potentials (xc as well as nonadditive Ts). The time-dependent KS orbitals were propagated for 50 000 steps with a Crank−Nicholson

⟨α(ω)⟩ =

1 Ns

N

∑I s αI(ω)

. To determine the proporSnapshots

tionality constant between the imaginary inverse dielectric constant and the imaginary frequency-dependent polarizability, the sum rule for the polarizability was normalized to the experimental value in the range of 0−25 eV. Each frequency5078

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083

Letter

The Journal of Physical Chemistry Letters

Subsystem TDDFT allows us to inspect the dynamical interactions between subsystems and allows us a glimpse into the cooperative interactions that arise when the liquid responds to an external perturbation. Formally exact decomposition of the interacting electronic response function of a system into subsystem contributions20,35,36 can be implemented. Namely

dependent polarizability is taken to be isotropic, that is, αIxx(ω) + αIyy(ω) + αIzz(ω) . 3

αI(ω) = The Fourier transformation was carried out including an artificial peak broadening (see the Supporting Information). In Figure 1, we show the comparison between our calculation and the most recent experiment for the real (ε1) and imaginary

Ns

χ=

∑ χIc

(1)

I

with each interacting subsystem response function given by a local interacting one-body term (χuI ) and a nonlocal many-body term χIc =

Ns

χIu

+

one⏟ body

∑ χIu KIJχJc J≠I  

many body

(2)

where hereafter superscripts c and u stand for “coupled” (i.e., employing the full many-body response) and “uncoupled” (i.e., employing only the one-body response). The subsystem TDDFT kernel for I ≠ J is given by KIJ(r1,r2,t−t′) =

δ(t − t ′) |r1 − r2|

+

δ 2Exc δρ(r1, t )δρ(r2, t ′)

+

δ 2Ts .20,21 δρ(r1, t )δρ(r2, t ′)

The

dipole polarizability is related to the response function by αxy I (ω) = −∫ x′χI(r′,r,ω)y dr dr′. The uncoupled subsystemlocal response function is computed with response equations that include occupied-virtual orbital excitations of only one subsystem. Alternatively, to determine χuI in the real-time approach, the embedding potential is computed at each time step, only updating the subsystem time-dependent density and leaving the density of the surrounding subsystems frozen at t = 0. To better understand the nature of the many-body interactions invoked in this work, we recall that the manybody terms include couplings between subsystem excitations at the level of two and higher bodies37

Figure 1. Real (ε1) and imaginary (ε2) parts of the frequencydependent dielectric constant of liquid water. IXS stands for inelastic X-ray scattering data from ref 11. In red are our subsystem TDDFT results, which were computed with a large broadening parameter (η = 0.5 eV) to better compare to the experiment.

(ε2) parts of the frequency-dependent dielectric constant. Although differences are evident, our subsystem TDDFT calculation recovers the overall trend and also reproduces some interesting features. Similarly to previous simulations based on supramolecular TDDFT of water-32,34 the computed ε2 is consistently red-shifted; however, the overall multimodal shape is reproduced. For ε1, our subsystem TDDFT simulations satisfactorily reproduce the ω = 0 limit, as well as the peak maximum.

χIc − χIu = χIu KIJ χJu + χIu KIJ χJu KJK χKu ...

(3)

where we have employed Einstein’s summation rule.

Figure 2. One-body (i.e., employing only the uncoupled response functions) and full many-body (i.e., employing eqs 1 and 2)) real (ε1) and imaginary (ε2) parts of the frequency-dependent dielectric constant of liquid water. Filled curves highlight the difference between the one-body and the full many-body result. 5079

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083

Letter

The Journal of Physical Chemistry Letters

Figure 3. (a) JDOS of the first excitation energy of the gas phase (dotted line) and first absorption band of liquid water (blue solid line). (b) First excitation energy of isolated water molecules against the symmetric stretch mode. (c) First absorption band of liquid water molecules against the symmetric stretch mode. (d) First absorption band of liquid water molecules against the environmental order parameter. Linear fits of the correlation scatter plots (red lines) are shown. Insets (b−d) are represented by hexagonal bins. Darker hexagons correspond to higher count value.

dielectric constant but do not change the position of the first peak in the optical spectrum. In Figure 3, we summarize our results pertaining to the first absorption band. Inset (a) of the figure displays smooth histograms of the first excitation energy on the x-axis and the subsystem count on the y-axis for the liquid (in blue) and for the vapor (dotted line). These histograms can be interpreted as the vibrationally modulated joint density of states (JDOS). The gas-phase data is generated with water molecules having the same geometry as the liquid, but their electronic structure is computed in the absence of the environment with the ADF39 software employing the PBE xc functional and a quadruple-ζ Slater-type orbital basis set. From the JDOS, we evince that the gas-phase first excitation energy is sharply centered around the average (there are 640 data points in the histogram) at ωiso 1 = 6.0 eV (with a resulting standard deviation of σiso 1 = 0.07 eV and absorption onset at 5.8 eV), while the subsystem TDDFT places the band maximum at ωu1 = 6.4 eV, featuring a very broad range of excitation energies (σu1 = 0.4 eV, onset at 5.4 eV). Such a broadening (σu1 ≫ σiso 1 ) is entirely due to interactions of the water molecules with their environment. Because the environment is heterogeneous and dynamic in the liquid, a broad distribution of excitation energies emerges. Figure 3a clearly shows that the onset of the absorption band (Urbach tail) is red-shifted compared to the isolated water molecule case by 0.4 eV, again due to vibronic coupling to environmental degrees of freedom. To justify the above claim, we provide in Figure 3b,c the scatter plot of the computed excitation energy for each of the isolated and embedded (liquid-phase) water molecules against the symmetric stretching degree of freedom. We find that the excitation energy for the isolated water molecules almost perfectly anticorrelates with the stretching mode. We also find a similar anticorrelation with the HOMO−LUMO energy difference. A sanity check shows that the HOMO−LUMO energy differences for the isolated (gas-phase) molecules

It is evident from Figure 2 that the many-body contributions are cooperative at low frequencies and anticooperative at high frequencies. Although it is difficult to pinpoint the specific reasons for this behavior, we note that increased spectral weights at low frequencies are typical of excitonically coupled systems.38 We should remark that the many-body term in eq 2 must be associated with an overall zero sum rule.37 Thus, if there is a many-body enhancement in one region of the spectrum, then there must be a many-body depletion in another region of the spectrum. In addition, the predicted many-body enhancement at low frequencies is consistent with the finding that the refractive index in water increases with increasing pressure, ascribed to an increase of the oscillator strengths with pressure.7 From our results, we estimate the refractive index as n = ε1(ω = 0) and obtain n = 1.25 and 1.30 for the one-body and the full many-body results. Thus, many-body effects increase n by about 4%. After having characterized the role of many-body dynamical interactions between water molecules in the liquid, we now analyze the first absorption band. Subsystem TDDFT places the average peak position of the first absorption band at ωu1 = 6.4 eV. Compared to the experimental 8.3−8.6 eV, it underestimates it by ∼2 eV. This is expected5 as we employ the PBE semilocal xc functional and the adiabatic approximation. Due to the more localized electronic structure compared to semilocal TDDFT of the supersystem, it is also expected that subsystem TDDFT yields a larger gap (semilocal TDDFT of the supersystem finds a gap of ∼5 eV34). For the same reasons, the subsystem TDDFT value for ε1(ω = 0) = 1.68 slightly underestimates the DFT-PBE value of 1.72.15 We find that the many-body terms have a negligible effect on the peak position of the first absorption band, which only shifts by about 0.02 eV when they are included. In other words, ωc1 − ωu1 = −0.02 eV. This points to a unique feature of liquid water, that is, the many-body terms affect the magnitude of the 5080

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083

Letter

The Journal of Physical Chemistry Letters

setting. For this, we also employed neural network regression (3 layers, 5 neurons per layer). The neural network regression analysis led to an average (statistical) correlation to the EOP in perfect agreement with the linear regression result of 0.65. We further explain the neural network analysis in the Supporting Information. This confirms that the nonlinear combinations of the initial order parameters (i.e., by slicing) describe fully the nonlinearity in the relation between the EOP and ωu1. Analysis of the major components of EOP highlights two important sets of parameters: first, the number of accepted hydrogen bonds and, second, the distances to the four closest oxygens in the environment. That is, the first solvation shell surrounding a water molecule determines the largest part of the first excitation energy peak position. In Figure 4, we report the

computed with ADF or eQE are essentially identical, with a small ∼0.1 eV offset due to the pseudopotentials involved in the eQE computation. Conversely, the first excitation of the embedded water molecules does not correlate at all with this internal degree of freedom. Inspection of Figure 3b,c unequivocally determines that the symmetric stretching, which was a determinant for the spectrum of gas phase-water, no longer plays a role in the spectrum of the liquid. The other local degrees of freedom (asymmetric stretch and bond angle) do not correlate with the absorption peak of either the gas or liquid phase. The same conclusions are reached when a similar analysis is carried out on the Kohn−Sham HOMO−LUMO gap of the embedded water molecules. Thus, we set out to find the degree of freedom describing the environment surrounding of a water molecule in the liquid that correlates the most to the first absorption peak of each subsystem. Among the ubiquitous order parameters (such as the number of donated/accepted hydrogen bonds, coordination number, and tetrahedral order parameter27), we found that only the number of accepted hydrogen bonds features a significant correlation to the first peak position. However, the computed correlation value (correlation =

σ(X , Y ) σ(X , X )σ(Y , Y )

, with

σ being the variance and X = excitation energy and Y = degree of freedom) is only 0.31, and thus, it is not large enough to confidently attribute the presence of a correlation. In Figure 3d, we show the scatter plot of the first excitation energy of each subsystem against a composite order parameter, which we call the “environment order parameter” (EOP), which includes the degrees of freedom of molecules in the environment. The EOP was constructed from 45 independent descriptors. They include O−O distances with the closest 9 waters and 12 respective angles with the hydrogen atoms (i.e., the O1−O2−H2 and O2−O1−H1 angles). In addition, they also include the tetrahedral order parameter,27 accepted and donated number of hydrogen bonds, coordination number, and 24 additional parameters estimating the lone pair positions. These initial descriptors were combined in nonlinear ways (for example, slicing the O−O distance in this way: O−O distance × (θ − θ0), with θ0 being a constant angle value) to generate 339 linearly dependent environmental degrees of freedom. These were reduced by singular value decomposition of their covariance matrix (i.e., a principal component analysis, or PCA) to only 25. Prior to performing PCA, the order parameters were filtered and ranked by correlation to ωu1. Overall, this procedure led us to a massive reduction in dimensionality, especially in view of the fact that the total number of possible binary nonlinear combinations of order parameters with 3 slices per set (i.e., 3 values of θ0) with 63 molecules in the environment equals 63 × 3 × 62 × 3 = 35154 parameters. Although we do not expect to have a complete and full description of the environment with only the given 25 principal components, these components are likely among the most important order parameters describing the environment. The 25 linearly independent order parameters resulting from the PCA were linearly combined to yield an EOP in such a way to maximize the correlation to the first excitation energy of each subsystem. The maximization was carried out by linear regression. EOP results in a correlation to the first absorption band peak position of 0.65. This is more than double the original best value from simple order parameters. We have also inspected whether linear regression is appropriate in this

Figure 4. Correlation coefficient values between ωu1 and the first nine oxygen−oxygen distances.

correlation coefficient of the first nine oxygen−oxygen distances to the first absorption band of liquid water. We clearly notice that the fourth and third distances correlate more strongly than the others. This further reinforces our assessment above, particularly because the distances are correlated with each other. That is, a closer fourth oxygen atom will, on average, be associated with a smaller first solvation shell, which is associated with shorter oxygen−oxygen distances. In sum, the first solvation shell in large part determines the absorption peak of water in the liquid phase. In Figure S4 of the Supporting Information, we show that 95% of the EOP is determined within the first solvation shell (as compared to the oxygen− oxygen radial pair distribution function). We should point out that the short-range nature of the EOP is nontrivial. In disordered molecular systems (like water), the vibrational spectrum is known to be short-ranged.4 However, the same general statement cannot be argued for electronic spectra. Electronic transitions interact at long-range and can have a fairly delocalized character even though the medium is disordered. Despite this, in this work, we show that at least for the first absorption band also the optical response of water is short-ranged (as determined by the EOP). We do not expect to recover a similar behavior for the higher-lying transitions in water (occurring at energies above 9 eV, which are now subject of analysis in our laboratory and will be presented in a future work). Our analysis explains the not large enough correlation to the total number of hydrogen bonds because the hydrogen bond definition is too simplistic and only combines a distance cutoff criterion in conjunction with an angle cutoff, therefore missing the complexity of the interaction. In the Supporting 5081

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083

Letter

The Journal of Physical Chemistry Letters

(6) Ziaei, V.; Bredow, T. Red and blue shift of liquid water’s excited states: A many body perturbation study. J. Chem. Phys. 2016, 145, 064508. (7) Pan, D.; Wan, Q.; Galli, G. The refractive index and electronic gap of water and ice increase with increasing pressure. Nat. Commun. 2014, DOI: 10.1038/ncomms4919. (8) Chan, W.; Cooper, G.; Brion, C. The electronic spectrum of water in the discrete and continuum regions. Absolute optical oscillator strengths for photoabsorption (6−200eV). Chem. Phys. 1993, 178, 387−400. (9) Elles, C. G.; Rivera, C. A.; Zhang, Y.; Pieniazek, P. A.; Bradforth, S. E. Electronic structure of liquid water from polarization-dependent two-photon absorption spectroscopy. J. Chem. Phys. 2009, 130, 084501. (10) Heller, J. M.; et al. Collective oscillation in liquid water. J. Chem. Phys. 1974, 60, 3483. (11) Hayashi, H.; Hiraoka, N. Accurate Measurements of Dielectric and Optical Functions of Liquid Water and Liquid Benzene in the VUV Region (1−100eV) Using Small-Angle Inelastic X-ray Scattering. J. Phys. Chem. B 2015, 119, 5609−5623. (12) Hahn, P. H.; Schmidt, W. G.; Seino, K.; Preuss, M.; Bechstedt, F.; Bernholc, J. Optical Absorption of Water: Coulomb Effects versus Hydrogen Bonding. Phys. Rev. Lett. 2005, 94, 037404. (13) Hermann, A.; Schmidt, W. G.; Schwerdtfeger, P. Resolving the Optical Spectrum of Water: Coordination and Electrostatic Effects. Phys. Rev. Lett. 2008, DOI: 10.1103/PhysRevLett.100.207403. (14) Mata, R. A.; Stoll, H.; Cabral, B. J. C. A Simple One-Body Approach to the Calculation of the First Electronic Absorption Band of Water. J. Chem. Theory Comput. 2009, 5, 1829−1837. (15) Lu, D.; Gygi, F.; Galli, G. Dielectric Properties of Ice and Liquid Water from First-Principles Calculations. Phys. Rev. Lett. 2008, DOI: 10.1103/PhysRevLett.100.147601. (16) Blase, X.; Boulanger, P.; Bruneval, F.; Fernandez-Serra, M.; Duchemin, I. GW and Bethe-Salpeter study of small water clusters. J. Chem. Phys. 2016, 144, 034109. (17) Jacob, C. R.; Neugebauer, J. Subsystem density-functional theory. WIREs: Comput. Mol. Sci. 2014, 4, 325−362. (18) Wesolowski, T. A.; Shedge, S.; Zhou, X. Frozen-Density Embedding Strategy for Multilevel Simulations of Electronic Structure. Chem. Rev. 2015, 115, 5891−5928. (19) Krishtal, A.; Sinha, D.; Genova, A.; Pavanello, M. Subsystem Density-Functional Theory as an Effective Tool for Modeling Ground and Excited States, their and Many-Body Interactions. J. Phys.: Condens. Matter 2015, 27, 183202. (20) Casida, M. E.; Wesolowski, T. A. Generalization of the KohnSham Equations with Constrained Electron Density Formalism and Its Time-Dependent Response Theory Formulation. Int. J. Quantum Chem. 2004, 96, 577−588. (21) Neugebauer, J. Chromophore-Specific Theoretical Spectroscopy: From Subsystem Density Functional Theory to Mode-Specific Vibrational Spectroscopy. Phys. Rep. 2010, 489, 1−87. (22) Krishtal, A.; Ceresoli, D.; Pavanello, M. Subsystem Real-Time Time Dependent Density Functional Theory. J. Chem. Phys. 2015, 142, 154116. (23) Genova, A.; Ceresoli, D.; Pavanello, M. Avoiding Fractional Electrons in Subsystem DFT Based Ab-Initio Molecular Dynamics Yields Accurate Models For Liquid Water and Solvated OH Radical. J. Chem. Phys. 2016, 144, 234105. (24) Genova, A.; Ceresoli, D.; Krishtal, A.; Andreussi, O.; DiStasio, R., Jr.; Pavanello, M. eQE: An Open-Source Densitiy Functional Embedding Theory Code for the Condensed Phase. Int. J. Quantum Chem. 2017, 117, e25401. (25) Rappe, A. M.; Rabe, K. M.; Kaxiras, E.; Joannopoulos, J. D. Optimized pseudopotentials. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 1227−1230. (26) Vanderbilt, D. Soft self-consistent pseudopotentials in a generalized eigenvalue formalism. Phys. Rev. B: Condens. Matter Mater. Phys. 1990, 41, 7892−7895.

Information, we provide further details about the procedure to obtain the EOP. In conclusion, we carried out real-time subsystem TDDFT simulations of liquid water and compared the results to gasphase water. We show that, although many-body excitonic effects have no effect on the position of the first absorption band peak, they dramatically enhance the spectral intensities in the low-frequency range and deplete the intensities at high frequencies. The simulations also predict that the liquid’s Urbach tail to be 0.4 eV red-shifted compared to the gas phase due to the strong coupling of the first absorption band with the environment within the first solvation shell. This reproduces the experimental observation and the original theoretical model of ref 12. Our analysis finds that interaction with the environment dramatically broadens the JDOS and blue shifts the band peak also by 0.4 eV. In sum, the results reproduce semiquantitatively the experimental red shift of the Urbach tail and the blue shift of the band peak.

■

ASSOCIATED CONTENT

* Supporting Information S

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpclett.7b02212. Computational details for the dielectric constant, environment order parameter, and comparison of linear regression and neural network regression (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Michele Pavanello: 0000-0001-8294-7481 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We thank Robert A. DiStasio, Jr. for helpful discussions, particularly for pointing out the need to slice the O−O distances to enhance correlations. This material is based upon work supported by the U.S. Department of Energy, Office of Basic Energy Sciences, under Award Number DE-SC0018343.

■

REFERENCES

(1) Cabral do Couto, P.; Chipman, D. M. Insights into the ultraviolet spectrum of liquid water from model calculations: The different roles of donor and acceptor hydrogen bonds in water pentamers. J. Chem. Phys. 2012, 137, 184301. (2) Cabral do Couto, P.; Chipman, D. M. Insights into the ultraviolet spectrum of liquid water from model calculations. J. Chem. Phys. 2010, 132, 244307. (3) Iftimie, R.; Tuckerman, M. E. Decomposing total IR spectra of aqueous systems into solute and solvent contributions: A computational approach using maximally localized Wannier orbitals. J. Chem. Phys. 2005, 122, 214508. (4) Medders, G. R.; Paesani, F. On the interplay of the potential energy and dipole moment surfaces in controlling the infrared activity of liquid water. J. Chem. Phys. 2015, 142, 212411. (5) Garbuio, V.; Cascella, M.; Reining, L.; Sole, R. D.; Pulci, O. Ab Initio Calculation of Optical Spectra of Liquids: Many-Body Effects in the Electronic Excitations of Water. Phys. Rev. Lett. 2006, DOI: 10.1103/PhysRevLett.97.137402. 5082

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083

Letter

The Journal of Physical Chemistry Letters (27) DiStasio, R. A.; Santra, B.; Li, Z.; Wu, X.; Car, R. The individual and collective effects of exact exchange and dispersion interactions on the ab initio structure of liquid water. J. Chem. Phys. 2014, 141, 084502. (28) Lee, H.-S.; Tuckerman, M. E. Structure of liquid water at ambient temperature from ab initio molecular dynamics performed in the complete basis set limit. J. Chem. Phys. 2006, 125, 154507. (29) Lee, H.-S.; Tuckerman, M. E. Dynamical properties of liquid water from ab initio molecular dynamics performed in the complete basis set limit. J. Chem. Phys. 2007, 126, 164501. (30) Yabana, K.; Bertsch, G. F. Time-dependent local-density approximation in real time. Phys. Rev. B: Condens. Matter Mater. Phys. 1996, 54, 4484−4487. (31) Chen, W.; Ambrosio, F.; Miceli, G.; Pasquarello, A. Ab initio Electronic Structure of Liquid Water. Phys. Rev. Lett. 2016, DOI: 10.1103/PhysRevLett.117.186401. (32) Lembarki, A.; Chermette, H. Obtaining a Gradient-Corrected Kinetic-Energy Functional from the Perdew-Wang Exchange Functional. Phys. Rev. A: At., Mol., Opt. Phys. 1994, 50, 5328. (33) Bertsch, G. F.; Iwata, J.-I.; Rubio, A.; Yabana, K. Real-space, realtime method for the dielectric function. Phys. Rev. B: Condens. Matter Mater. Phys. 2000, 62, 7998−8002. (34) Tavernelli, I. Electronic density response of liquid water using time-dependent density functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2006, DOI: 10.1103/PhysRevB.73.094204. (35) Pavanello, M. On the Subsystem Formulation of LinearResponse Time-Dependent DFT. J. Chem. Phys. 2013, 138, 204118. (36) Neugebauer, J. Couplings between electronic transitions in a subsystem formulation of time-dependent density functional theory. J. Chem. Phys. 2007, 126, 134116. (37) Krishtal, A.; Pavanello, M. Revealing Electronic Open Quantum Systems with Subsystem TDDFT. J. Chem. Phys. 2016, 144, 124118. (38) Onida, G.; Reining, L.; Rubio, A. Electronic excitations: densityfunctional versus many-body Green’s-function approaches. Rev. Mod. Phys. 2002, 74, 601−659. (39) te Velde, G.; Bickelhaupt, F. M.; Baerends, E. J.; van Gisbergen, S. J. A.; Fonseca Guerra, C.; Snijders, J. G.; Ziegler, T. Chemistry with ADF. J. Comput. Chem. 2001, 22, 931−967.

5083

DOI: 10.1021/acs.jpclett.7b02212 J. Phys. Chem. Lett. 2017, 8, 5077−5083