Modeling Liquid Photoemission Spectra: Path ... - ACS Publications

Sep 21, 2016 - Dynamics Combined with Tuned Range-Separated Hybrid ... tional price. The valence photoelectron spectrum of water is a convenient test ...
1 downloads 0 Views 843KB Size
Subscriber access provided by CORNELL UNIVERSITY LIBRARY

Article

Modeling Liquid Photoemission Spectra: Path Integral Molecular Dynamics Combined with Tuned Range-Separated Hybrid Functionals Daniel Hollas, Eva Muchova, and Petr Slavicek J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.6b00630 • Publication Date (Web): 21 Sep 2016 Downloaded from http://pubs.acs.org on September 27, 2016

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Modeling Liquid Photoemission Spectra: Path Integral Molecular Dynamics Combined with Tuned Range-Separated Hybrid Functionals Daniel Hollas, Eva Muchová and Petr Slavíček* University of Chemistry and Technology Prague, Department of Physical Chemistry, Technická 5, 16628 Prague 6, Czech Republic KEYWORDS path integral molecular dynamics; photoemission spectroscopy; range-separated hybrid functional; water; photoionization; ab initio methods

We present a computational protocol for modeling valence photoemission spectra of liquids. We use water as experimentally well-characterized model system and we represent its liquid state by larger finite-size droplets. The photoemission spectrum is evaluated for an ensemble of structures along molecular dynamics simulation. The nuclear quantum effects are accounted for by ab initio based path-integral molecular dynamics simulations greatly accelerated with the so called colored noise thermostat (PI+GLE) method. The ionization energies for the valence electrons are evaluated as orbital energies of optimally tuned range-separated hybrid functionals (OT-RSH). This approach provides Koopmans-type ionization energies including relaxation energy. We

ACS Paragon Plus Environment

1

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 2 of 34

show that the present protocol can describe quantitatively the valence photoemission spectrum of liquid water, i.e. the position, shapes and widths of the photoemission peaks. With the PI+GLE simulations, even the subtle isotope effects recently observed experimentally can be modeled. The electronic properties of finite-size droplets are shown to converge rapidly to that of liquids. We discuss the importance of proper tuning of the range-separation parameter in the OT-RSH as well as possible sources of errors in our simulations. The present approach seems to be a viable route to modelling photoemission spectra of liquids, especially in conjunction with efficient implementation of density functional methods on graphical processing units.

Introduction Photoemission spectroscopy is the most straightforward approach to study electronic structure of molecules and solid state materials.1 Practical application of photoemission spectroscopy for solutions and liquid interfaces is, however, relatively new.2 Technology developed for the solid state can be in principle used for non-volatile liquids (e.g. for ionic liquids3), but high vapor pressure of, e.g. water, makes the detection of outgoing electrons difficult. In the last decade, the field has witnessed a fast growth due to the introduction of new experimental setups (such as liquid microjet technique4 or near-ambient pressure photoelectron spectroscopy5) and wider availability of high energy photon sources. The experimental development also challenges the theory as the data have to be interpreted with ab initio based techniques. Modeling photoemission spectra in liquids is difficult from several perspectives. First, liquids are disordered systems and one cannot represent the liquid state by a single structure. Second, as charge of a solute is changed within the process of photoemission, even rather distant molecules contribute to energetics due to long-range nature of Coulombic forces. Finally, we need to describe a large number of excited final electronic states formed upon the electron ejection. In

ACS Paragon Plus Environment

2

Page 3 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

this work, we focus on modeling valence photoelectron spectra (PE) of liquid water using pathintegral based molecular dynamics (PIMD) simulations combined with optimally tuned rangeseparated hybrid functionals (OT-RSH), representing the liquid state by finite-size droplets. We show that this approach addresses most of the issues raised above and the method provides a quantitative estimate of photoemission spectra at affordable computational price. Valence photoelectron spectrum of water is a convenient test bed for developing new computational algorithms in liquid phase. Experimental data are available for the gas phase molecules,6 finite size clustersclusters7 as well as for liquid water;8 the data were measured even for isotopically substituted D2O.9 Theoretical studies on photoelectron spectra of liquid water are surprisingly scarce.10 Many studies are focusing on ab initio calculations of ionization energies of finite size water aggregates,7,11–16 typically without accounting for the inhomogeneous broadening of the spectra. However, important physical information is contained not only in the position of the photoelectron peaks, but also in peak widths and intensities. All these quantities reflect the structure of the liquid in the ground state as well as its electronic structure and shape of the potential energy surface (PES) in the ionized states. In our approach, we model the photoelectron spectra of liquid water using the so called reflection principle17 combined with path-integral molecular dynamics simulations (RP-PIMD). In the reflection principle, the ground state density is projected onto an ionized (or more generally excited) state. This method provides only low resolution spectra ignoring the vibrational resonances. On the other hand, it provides reliable spectral shapes at acceptable computational costs for extended molecular systems. The RP-PIMD approach was used before for calculating electronic spectra of molecules18–22 and molecular clusters;14,21,23–25 here, we apply the approach for liquid phase. The RP-PIMD scheme, though semi-classical because of the

ACS Paragon Plus Environment

3

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 34

reflection principle, is capable to capture the quantum nature of water which is now becoming widely appreciated26 and these effects should be seen also in the photoelectron spectra. Note at this point that path-integral simulations are computationally more demanding than the corresponding classical simulations. At the same time, recent development in the field of ab initio based PIMD suggests that the computational overhead can be largely eliminated.27–29 Here we use PIMD simulations accelerated with the so called colored noise thermostat.30 An obvious problem connected to the calculations of photoelectron spectra in the condensed phase is the necessity to include a large number of ionized states; each peak (e.g. the HOMO peak) consists of many ionized states which are distributed in energy. The arising question is if we can efficiently calculate the energy of these states. The ∆SCF or ∆KS approaches provide only the lowest ionization energy which is hardly satisfactory for liquids. There exists a whole variety of techniques to obtain even the higher states, e.g. the energies can be obtained from the excited states of ions, various Green’s function based methods can be used or highly accurate methods such as equation of motion coupled cluster theories can be applied.31–33 However, these approaches are computationally demanding when applied for the liquid phase. A thrifty route to ionization energies is represented by Koopmans’ theorem which provides the ionization energies as negative orbital energies of the neutral precursor of the ion. As liquid water is now routinely modeled at the ab initio level, it is natural to use the on-the-fly orbital energies for calculating the spectra. Yet there is a well-known drawback of the approach related to the missing relaxation effect in the Koopmans’ approach. Approximate density functional methods do not, in general, fulfill the Koopmans’ theorem. The Kohn–Sham orbital energies can still be used for modeling PE spectra, e.g. using G0W0 correction within the many body perturbation theory framework.10,34,35 A simple approach to

ACS Paragon Plus Environment

4

Page 5 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

improve the Koopmans’ theorem is based on the concept of tuning the range-separated hybrid functionals so that the orbital energy εi can be identified with the ionization energies of the i-th electron. The optimal tuning is achieved by setting the range parameter  so that the ∆KS and HOMO orbital energies are equal   =   −    The above condition leads to a value of HOMO binding energy which involves orbital relaxation. It has been shown with numerical and theoretical arguments that the orbital energies of the OT-RSH functionals represent faithfully the whole valence spectrum of ionization energies.36,37 Alternatively, one can include LUMO orbital energy and electron affinity to the optimized function to obtain superior electronic quantities such as band gaps and excitation energies.38–40 The combination of ab initio PIMD simulations with OT-RSH seems to be a convenient approach for modeling photoelectron spectra. As the orbital energies are provided on-the-fly during the simulation run, there is no extra cost with recalculating the ionization energies for each snapshot. At the same time, we obtain all of the ionization energies.18 The general questions addressed in this work can be summarized as follows: (i) Is the shape of liquid water PE spectrum correctly reproduced with the RP-PIMD Koopmans-based approach? (ii) What is the role of nuclear quantum effects on the spectrum? (iii) Can we also correctly calculate the positions of the peaks in liquid water? (iv) How does the PE spectrum converge with the number of water molecules? The application of the OT-RSH should not be accepted without caution. For example, can we use the optimal range parameter at a single snapshot for all the other geometries? Previous study indicates that the range parameters do indeed depend on a particular molecular geometry.18 Can we also transfer the range parameter optimized for an isolated

ACS Paragon Plus Environment

5

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 6 of 34

molecule or a small cluster and use it for liquids? We attempt to answer all the above questions by comparing the simulated data with experimentally measured spectra.

Computational Methods Photoelectron spectra were modeled within the reflection principle. In this semi-classical approach, the ground state nuclear density multiplied by the transition dipole moment is reflected onto the final state (ionized state in the present case).19,22,25 The ground state density translates into cross section for the photoemission via the relation    =

 &     !" #$", # ' − (, −  ) !" , 3ℏ  

where   !"  is the ground state nuclear density, E is the photon energy (typically in the X-ray range),  is the kinetic energy of the outgoing electron and (, is the binding energy of the i-th

electron in a molecule at geometry !" (semi-classical coordinate space representation)19 &

and #$", # is the transition dipole moment between the final (corresponding to an ionized molecule and outgoing electron with energy  ) and initial electronic states (denoted as a) at a

given geometry. We then display the resulting spectrum as a function of electron binding energies. In our calculations, we assume that the transition dipole moment is independent of geometry (Franck–Condon approximation) and is furthermore equal for all ionized states. The input for our approach is thus represented by two quantities: (i) Nuclear density   !" , which is generated via molecular dynamics techniques (ii) Ionization energies for each electron. We discuss technical details of the calculations of these quantities further in this section. Generating ground state densities for finite size models. We represent liquid water by finitesize molecular clusters, using however rather large aggregates up to 100 water molecules. Such

ACS Paragon Plus Environment

6

Page 7 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

clusters cannot completely mimic the structure of liquid water, but they are large enough to capture the most important features. Using an efficient GPU-based DFT implementation in TeraChem code,41 we were able to perform all MD simulations on a single GPU. Larger water clusters have extremely rich potential energy surface and their description is rather complicated. For the ground state sampling we used either classical MD (in this work, the term classical means that nuclei are treated classically but electrons are treated quantum mechanically within DFT) or path-integral based MD simulations. The MD simulations were performed for cluster sizes ranging from isolated water to (H2O)100. Forces in molecular dynamics were obtained on-the-fly with the LRC-ωPBE functional42 (further denoted as ωPBE) as implemented in the Terachem code,41 with 6-31g* basis set and without dispersion corrections. Such a choice of ab initio method is convenient as we use the same method to calculate ionization energies. Therefore, no further recalculation of ionization energies is needed and we get the spectra on-the-fly. The optimal range-separation parameter was evaluated for each cluster size from a series of snapshots from pilot MD simulations with empirical force field as described below. The temperature in all simulations was set to 300 K and the time step was set to 20 a.u (~0.5 fs). All dynamical simulations were performed using our in-house code ABIN.43 Path integral simulations can be quite demanding, especially in the ab initio context. Typically around 30 beads are needed for reasonably converged results for liquid water at room temperature for properties such as radial distribution functions and quantum kinetic and potential energies.26 However, it has been shown that the convergence can be significantly accelerated with the thermostat based on Generalized Langevin Equation (GLE) formalism.27,30 Using this so called PI+GLE27 approach, only 4 beads were needed to obtain converged densities in our simulations. The parameters for the GLE thermostat were taken from an online database.44 The

ACS Paragon Plus Environment

7

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 8 of 34

length of the PI+GLE simulations was typically 10 ps after the system was equilibrated for several picoseconds. Classical molecular dynamics simulations were performed with the velocity Verlet algorithm and Nosé–Hoover thermostat. The length of simulations was 15 ps. Calculations of ionization energies. The ionization energies were calculated at the ab initio level employing density functional theory with OT-RSH functionals. The orbital energies were collected for each geometry from the MD simulations using ωPBE42 functional with 6-31g* basis set. For PI+GLE simulations the total number was 80 000 geometries, for classical simulations, the total number of geometries was 30 000 geometries. Collected orbital energies were processed in the form of a histogram without Gaussian broadening of the lines. For a detailed comparison with the experimental spectra, the histograms were fitted to Gaussian functions as described in the results section. We optimized the value of the range-separation parameter for different structures and geometries. The tuning was performed for various cluster sizes, ranging from H2O to (H2O)100. For smaller water clusters (H2O)n, with n in the range 1–20, we tuned the parameter for a single optimized structure taken from Cambridge cluster database.45 These structures were re-optimized at the ωPBE(ω=0.45)/6-31g* level. For clusters (H2O)10 to (H2O)100, the tuning was performed for a set of 50 geometries obtained from PI+GLE simulations using q-TIP4P/F force field.46 This force field was specially designed for path-integral simulations so that the nuclear quantum effects are not double counted. The values of the range-separation parameter that were subsequently used in ab initio molecular dynamics were obtained as an average for the set of geometries.

Results and Discussion

ACS Paragon Plus Environment

8

Page 9 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Optimal Tuning of the RSH Functionals. The final goal of this work is to simulate photoemission spectrum of liquid water using the orbital energies from the range-separated DFT methods (ωPBE functional in the present case). The electronic structure method needs to be first calibrated and benchmarked. The first question to be addressed is what value of the parameter could be used for the simulations. We have previously tested tuning of the range parameter for small molecules in the context of X-ray absorption spectroscopy.18 We can conclude that (i) the parameter is different for different molecules (ii) the geometry variations of the optimal range parameter are relatively small. Molecular clusters with different numbers of monomeric units represent different molecules; the optimal range-separation parameter should thus be size dependent. Note at this point that with the range-separation parameter tuning procedure we effectively include orbital relaxation effect into the orbital energy. Using incorrect rangeseparation parameter means adding an undefined fraction of the orbital relaxation. Fig. 1 shows the dependence of optimal range-separation parameter on cluster size for the ωPBE functional. The dependence is very strong for small clusters as can be seen from the data for optimal geometries (red dots in Fig. 1). For larger clusters, for 50 geometries from PI+GLE/q-TIP4P/F simulations, we observe that the value of range-separation parameter converges with increasing the cluster size to a value close to 0.40. We can see at the same time that there is a spread between different geometries with a standard deviation around 0.03. The mean values of the range-separation parameter and standard deviations are given in Supporting Information as Table S1. For further work, we used the value of ω=0.45 to describe the convergence of the photoemission spectra with increasing of the cluster size (for cluster sizes n=10–100). We opted for a single value of range separated parameter for cluster sizes n=10–100 as we aim to demonstrate the changes of spectral shape with increasing the cluster size rather

ACS Paragon Plus Environment

9

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 34

than predict the precise positions of the peaks. We used a converged value of ω=0.40 for the largest (H2O)100 cluster, ω=0.41 for (H2O)35 and ω=0.56 for isolated water for quantitative spectra modeling. We expect that the error caused by using a geometry-independent rangeseparation parameter to be minor with respect to other sources of errors, as discussed below.

Figure 1. Cluster size dependence of the range-separation parameter for the ωPBE functional. The values of optimally tuned range-separation parameter for cluster structures taken from Cambridge cluster database45 are shown as red dots, while black dots indicate average values and standard deviations of range-separation parameter for clusters containing from 10 to 100 waters for 50 geometries sampled from PI+GLE simulations using q-TIP4P/F force field.

The optimal tuning should lead us to the correct values of the ionization energies within a given electronic structure method and basis set, i.e. if the ∆KS ionization energy for the HOMO orbital is erroneous, the error will propagate to the higher states. It is clearly needed to evaluate the quality of the electronic structure method used for the simulations. We tested the performance of the ωPBE functional for smaller water clusters where higher level calculations are feasible. The results of benchmark calculations are collected in Table 1. We can make several

ACS Paragon Plus Environment

10

Page 11 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

observations. First, the optimally tuned ωPBE/6-31g* ionization energies are surprisingly close to the ionization energies calculated with the high-level CCSD(T)/cc-pVTZ method. This agreement stems probably from some sort of error compensation; anyone working with (most of the) DFT methods should be comfortable with this aspect. It is also clearly seen that the rangeseparation parameter optimized for the larger clusters (ω=0.40) does not perform well for the small clusters, i.e. the functional is detuned. Note also that the basis set used is relatively modest to allow for extended dynamical calculations. Generally, we cannot expect accuracy better than 0.2–0.4 eV in all our spectra as the absolute position of the peaks is concerned. Within this range, we can assume the electronic structure level used in this work to be appropriate.

Table 1. Vertical Ionization Energies (VIE) in eV for a water cluster containing up to 20 water molecules using various methods. The cluster structures were taken from Cambridge cluster database.45 ωPBE/6-31g*(ω=0.40)

ωPBE/631g*(tuned)

ωPBE/631+g*(tuned)

CCSD(T)/ccpVTZb

systém

VIE(∆SCF) VIE(Koop) VIE

ω

VIE

ω

VIE

H2Oa

12.43

11.29

12.43

0.56

12.71

0.51

12.61

(H2O)3

11.75

11.30

12.13

0.52

12.29

0.46

12.18

(H2O)6

11.22

10.99

11.60

0.48

11.73

0.42

11.67

(H2O)10 11.08

10.71

11.07

0.45

11.32

0.40

-

(H2O)15 11.11

10.69

11.01

0.44

11.25

0.39

-

(H2O)20 11.25

11.25

11.25

0.40

11.06

0.31

-

a

Experimental VIE for the gas phase water molecule is 12.60 eV.6

b

s and p diffuse functions from aug-cc-pVTZ basis set were added for oxygen atom.

ACS Paragon Plus Environment

11

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 34

Photoelectron Spectra Modelling for Finite-Size Models. Let us now focus on the simulated photoemission spectra. Electrons in water are placed into five orbitals: core orbital 1a1 (oxygen 1s orbital, 539.5 eV),47 inner valence orbital 2a1 (oxygen 2s orbital, 32.62 eV)8 and outer valence orbitals 1b2 orbital (18.78 eV),8 3a1 orbital (14.84 eV)8 and 1b1 orbitals (12.60 eV),6 corresponding to lone pairs and sigma bonds. Electronic properties of water are strongly affected by intermolecular interactions and water’s dielectric function. In liquid water, the position of the HOMO peak moves to 11.16 eV,8 i.e. the solvent shift is 1.44 eV.8 This is an unusually high value for liquids which can be rationalized by water’s dielectric properties. Importantly, the photoelectron spectrum of liquid water has also a rather different shape with respect to its gas phase analogue. The HOMO peak is narrow and high, while the HOMO−1 (3a1) peak is much wider and lower in intensity. Fig. 2A shows the photoelectron spectrum calculated for a (H2O)100 cluster along the PI+GLE dynamical simulation on the ωPBE(ω=0.40)/6-31g* potential energy surface. We assume that this cluster size can reasonably well represent the liquid structure. In order to compare the simulated spectrum with the experiment, we fitted both the experimental and theoretical spectrum into a set of Gaussian peaks. As we can see from Fig. 2 this fit may not be ideal for the experimental curve. On the one hand, the merged 1b1/3a1 peak can be better fitted with three Gaussian peaks;8 on the other hand, it would obscure the assignment of the peaks to particular states and made the analysis ambiguous. The visual comparison between the experimental and theoretical curves shows a good agreement which is confirmed by quantitative analysis provided below.

ACS Paragon Plus Environment

12

Page 13 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Figure 2. A The comparison of modeled photoemission spectrum from PI+GLE simulation with ωPBE(ω=0.40)/6-31g* of (H2O)100 cluster (black curve) and experiment (red dots). The fit of the experimental data is depicted as a red curve. Theoretical curves were fitted to four Gaussian functions, the experimental fit includes two additional Gaussian functions to cover effects beyond the single particle picture in the 2a1 state. The fit is depicted in pink and grey. The widths and positions of the peaks are collected in Table 2. B The theoretical spectrum was weighted with differential photoionization cross sections inferred from the liquid phase experiment. The spectrum was normalized so that the surface below experimental and theoretical curves is equal. The experimental curves were measured at 200 eV with a detector at 90o.48 In the original spectrum published by Winter et al.,8 the photoemission spectrum included secondary electrons. These electrons were now carefully subtracted.49 The x-axis shows binding energies, i.e. negative values of ionization energies.

ACS Paragon Plus Environment

13

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 34

The peak positions and widths for the experimental and theoretical spectra are collected in Table 2. Note here that the reported widths slightly depend on the fitting procedure; yet, it does not affect the comparison of theory and experiment. As can be seen in Fig. 2 and Table 2, the shape of the spectra obtained from the PI+GLE/ωPBE(ω=0.40)/6-31g* simulation is generally in an excellent agreement with the experiment at least for the outer valence electrons (the inner valence electrons will be discussed later). The peaks widths of the outer valence peaks (1b1, 3a1 and 1b2) collected in Table 2 almost coincide with the experimental spectra. Also the relative positions are in a good agreement. The theoretical spectra are shifted by 0.2–0.4 eV to lower binding energies compared with the experiment which is dominantly caused by missing the longrange polarization in our finite size calculations and approximate character of the electronic structure method. We discuss it in more details in section Possible sources of errors. For the inner valence (2a1) electron, the observed error is on the other hand dominated by the accuracy of the DFT method itself as the self-interaction error become more severe for more tightly bound electrons. In this case the peak is shifted to higher binding energy.

Table 2. Ionization energies and full widths at half maximum (FWHM) for experimental gas phase and liquid spectra and for theoretical spectra obtained from PI+GLE/ωPBE(ω=0.40)/631g* simulations of (H2O)100 cluster. ionization energies [eV] orbital

exp. gas

exp. liquid*

FWHM

PI+GLE exp. gas

exp. liquid*

(H2O)100 1b1

12.60

11.21

11.51

PI+GLE (H2O)100

0.30

1.51

1.46

ACS Paragon Plus Environment

14

Page 15 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

3a1

14.84

13.44

13.89

1.18

3.24

3.08

1b2

18.78

17.54

17.76

1.75

2.22

2.24

2a1

32.62

30.87

30.74

2.82

3.10

2.44

*Contrary to ref.8, we fitted each 1b1, 3a1 and 1b2 peaks to a single Gaussian function. The 2a1 peak was fitted to three Gaussian functions as in ref.8 to account for the presence of satellite states as described in the text.

The agreement with the experiment in Fig. 2A is not fully satisfactory concerning the relative intensities. This is caused by two interrelated aspects. In our calculations, we neglect the electronic factors, i.e. we ignore the differences in the transition dipole moments for the ionization from different molecular orbitals (MOs). The relative photoionization cross sections deviates from unity even for isolated water molecule.6 One could think of correcting the calculated spectra with differential cross sections measured (or calculated) for water monomer. However, we also do not account for anisotropy of the photoemission; the measurements are not performed at magic angle, and the measured intensities are thus affected by the anisotropy parameter. To make the experimental and simulated data comparable, we processed the modeled spectra with differential photoionization cross sections inferred from the liquid phase experiment for the experimental arrangement. The spectra processed in this way are shown in Fig. 2B. The whole spectrum now fits almost perfectly to the experiment in position, width and intensity for all peaks. The only disagreement is observed for the inner valence 2a1 electron. The 2a1 peak deserves a special attention. The calculated spectrum misses side peaks around the main 2a1 peak. The side signals represent the satellite peaks50 resulting from electron correlation. The satellite states clearly cannot be described within the one-electron model. We fitted the experimental peak to three Gaussian functions and we compared only the central one to our theoretical spectrum. Yet, the central peak is still wider in the experiment than in the

ACS Paragon Plus Environment

15

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 34

simulations. This again can stem from the contributions from satellite states, but can be also related to finite lifetime of the 2a1 states due to non-adiabatic processes or non-local autoionization process (the so called Intermolecular Coulomb Decay process).51 We shall also address a question whether the presented calculations are converged with respect to the cluster size. In another words, we ask at which point a molecular cluster becomes liquid as seen from the PE spectroscopy perspective. Fig. 3 shows the evolution of the spectrum with clusters size. The spectra were calculated from the orbital energies recorded during the PI+GLE/ωPBE/6-31g* simulations of clusters containing 1, 10, 20, 60 and 100 water molecules. For all clusters presented here we used the value of range-separation parameter ω=0.45. Our choice to compare the spectra obtained with a single value of range-separation parameter is meant to demonstrate the evolution of a spectral shape with increasing the cluster size, not to predict an accurate shift of the peak positions. As the range-separation parameter is not tuned, we do not account properly for orbital relaxation which contributes significantly to the solvent shift; the orbital relaxation can in other terms be expressed as electronic polarization of the environment. The spectra in which the range-separation parameter was tuned for each cluster size, e.g. which show the evolution of both spectral shape and peak positions, is presented in Supporting Information as Figure S2. Yet even more accurate tuning on geometry-by-geometry basis (not done in this work), would be very important when we focus specifically on the size evolution of the peak positions. We can observe the importance of proper tuning for absolute position of the spectra on the example of isolated water molecule. The position of the HOMO peak is 11.70 eV with the ω=0.45 while it amounts to 12.37 eV for optimally tuned rangeseparation parameter ω=0.56 (value close to the experimental value of 12.60 eV8). The tuning procedure, however, does not affect the width of the spectrum. We can observe that the

ACS Paragon Plus Environment

16

Page 17 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

simulations correctly reproduce the widening of the peaks due to the intermolecular interaction. We can thus perform a simulation with any reasonable choice of the range-separation parameter to get the overall structure of the spectrum while setting the absolute position of the spectrum via benchmark calculations on a limited number of structures. The gas phase spectrum and the liquid PE spectrum differ profoundly in the position, shape and relative intensities of the peaks. The monomer has well separated 1b1/3a1 peaks and these two peaks gradually merge. Within several tens of water molecules, the spectrum seems to be already converged with respect to its spectral shape not peak positions. We see, however, that the 1b1/3a1 peaks are still better separated than they are in liquid water. This stems from the finite size effects; even for the cluster with 100 water molecules, most of the molecules are placed on the surface.

Figure 3. Cluster size dependence of photoemission spectrum. The spectra were calculated from the orbital energies recorded during the PI+GLE/ωPBE/6-31g* simulation for clusters containing 1, 10, 20, 60 and 100 water molecules. The range-separation parameter was set to 0.45 in all cases, the spectra presented here were not fitted to Gaussian functions. The spectra were

ACS Paragon Plus Environment

17

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 34

normalized so that the surface below curves is equal. The spectra with optimally tuned rangeseparation parameter for each cluster size are presented in Supporting Information as Figure S2.

There are two cluster sizes for which we can directly compare our calculations with the experiment: isolated water molecule52 and cluster with a mean size of 35 water molecules.7 The simulated spectra for these cases are shown in Fig. 4A and Fig. 4B and the parameters of the spectra are summarized in Supporting Information Table S2. The spectra were again decomposed into Gaussian functions and the expansion coefficients were weighted by the experimental differential cross sections. With this, we can observe an almost perfect match of the spectral shape with the experiment.

Figure 4. A Valence photoelectron spectrum of isolated water obtained from PI+GLE simulation with ωPBE(ω=0.56)/6-31g* (black curve). The experimental spectrum is taken from ref.52 The

ACS Paragon Plus Environment

18

Page 19 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

experimental valence spectrum was obtained using 180 eV photon energy. B Outer valence photoelectron spectrum of (H2O)35 water cluster. The experimental values for average cluster size =35 are taken from ref.7 (red dots, the fit is depicted as a red curve). The theoretical spectrum for an isolated water molecule was normalized so that the area below the theoretical curve equals the area below the experimental one; for (H2O)35 the spectrum was corrected by using the valence cross sections taken from ref.7 The corresponding values are 0.68 for 3a1 orbital and 0.93 for 1b2 orbital. Both theoretical and experimental spectra were fitted to three Gaussian functions, the functions are depicted as pink (for experimental) and grey (for theoretical) curves. The spectrum was again normalized so that the surfaces below fitted curves are equal. Grey and pink curves show Gaussian functions used in the fit.

Nuclear Quantum Effects. Our simulation techniques relying on the PIMD simulations should account also for nuclear quantum effects, demonstrated, among other, by isotope effects. In their seminal paper on liquid photoemission, Winter et al.8 state that “Within the limits of the present experiment, the electronic structures of liquid H2O and D2O are indistinguishable.” The Winter group has later re-measured the spectrum for deuterated water and a small isotopic shift of 0.07 eV was observed for the gas phase peak.53 The measured peaks are also narrower for deuterated water by 0.1–0.2 eV. In our fit of the experimental data48 we see a very slight shift of peak positions of 0.03–0.07 eV (peaks for D2O are higher in energy). Most of the observed effect can be attributed to differences in zero point energy between water and its deuterated analogue. This effect is only partially covered in our simulations based on the reflection principle.54 To address the observed effect, we investigated water clusters and deuterated water clusters by PI+GLE simulations and compared them with the results from a classical simulation. We used here

ACS Paragon Plus Environment

19

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 34

smaller clusters, (H2O)60 and (D2O)60, that can still represent the liquid structure. The spectra are presented in Fig. 5 and the values of peak positions and widths (including experiment) are collected in Supporting Information Table S3. As we can see from Fig. 5 and Table S3, the shifts between the PI+GLE results for H2O and D2O are in a very good agreement with the measured experimental data48 processed with the same fitting procedure as the calculated data. The observed shifts are up to 0.07 eV. The change in FWHM is again in an excellent agreement with experimental data; peaks for D2O are about 0.07–0.3 eV narrower. The “classical” spectrum is shifted to lower energies and the peak widths are significantly lower than in the H2O and D2O case.

Figure 5. Comparison of photoemission spectra of (H2O)60 and (D2O)60 obtained from PI+GLE simulations at the ωPBE(ω=0.45)/6-31g* level. The spectra are compared to a spectrum obtained from a classical MD at the ωPBE(ω=0.45)/6-31g* level for (H2O)60. Note that the spectrum obtained from classical MD is the same for (H2O)60 and (D2O)60.

Possible Sources of Errors. Let us conclude with few notes on possible sources of errors in the above simulations. In all dynamical calculations, a modest 6-31g* basis set was used. Yet the

ACS Paragon Plus Environment

20

Page 21 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

peak positions and widths are only negligibly influenced when using larger basis set. The spectra obtained at PI+GLE/ωPBE level using 6-31g* and 6-31+g* basis set for H2O and (H2O)10 cluster are shown in Supporting Information as Figure S1. Electronic structure modelled by DFT also inherently includes an error which affects the peak positions. Based on studies for smaller molecules and also Table 1, the error amounts to 0.2 eV in increasing the binding energies. Optimal tuning of the range-separation parameter can be also done on geometry-by-geometry basis.18 In this case even more accurate peak positions can be achieved, however, the overall effect on the spectra is in the case of water clusters small. The comparison is demonstrated for (H2O)35 in Figure S3. Another difficulty is connected to the fact that even large clusters do not fully represent liquids. One could be tempted to use clusters generated from simulations with periodic boundary conditions. In this case we can cut smaller clusters and recalculate the ionization energies. Such clusters are left, however, with unphysically dangling O–H bonds. This leads to spurious results, in agreement with previous observations.55 Another option is to use QM/MM simulations. We performed simulations with a system containing 500 water molecules, the inner 60 molecules were treated at the PI+GLE/ωPBE(ω=0.45)/6-31g* level, the outer part was simulated by TIP3P potential. The obtained spectra are presented in Supporting Information (Figure S4 and Table S4). The results show that the peaks are shifted towards lower energies by 0.3–0.4 eV. This is not surprising as we miss the effect of the electronic polarization of the MM part of the cluster. The general shape of the spectra is essentially the same as for the pure water cluster. Finally, we estimate the error caused for the absolute position of the ionization energy by using finite-size models. In our model, we miss the polarization contribution from more distant water molecules. Solvent polarization has two parts: (i) Slow nuclear polarization and (ii) fast

ACS Paragon Plus Environment

21

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 34

electronic polarization. Fully in the spirit of Born–Oppenheimer approximation, we shall consider only the electronic polarization as the atomic nuclei move too slowly to adjust their position to a changing electronic structure upon photoionization. In our calculations, the polarization should be included for the water molecules buried in the cluster. This is not a completely obvious statement for ionization energies calculated within the Koopmans’ theorem, the missing relaxation effect in the Koopmans’ theorem based on HF orbital energies in fact suggests a lack of electronic polarization. With the OT-RSH, the effect is included. The polarization of solvent molecules beyond the first few solvation shells explicitly included in the simulations might be still important. As at this distance the effects are non-specific, using the dielectric models should be safe, i.e. we can estimate the solvent shift by considering the electrical work needed for charging a sphere in vacuum and in an environment of a certain dielectric constant (Born theory). Note that only the optical part of the dielectric constant should be included as to account only for the electronic polarization. The shift is then given as ∆VIE =

1 1 /& 11 − 5, 8 234

where ε0 is the permittivity of vacuum, εopt is optical dielectric constant of water (we used the value of 1.8 from ref56), and R is the radius of the cluster. We estimated the radius from the molecular volume calculated along the trajectory with the MSMS code.57 The correction obtained in this way was 0.4–0.9 eV for water cluster sizes ranging from 100 to 10 water molecules. The results indicate that including the long-range polarization via the Born model yields a shift in peak position towards higher binding energies. In the case of (H2O)100 cluster, the estimated shift of 0.4 eV improves the agreement with the experimental values for 1b1, 3a1 and 1b2 peaks.

ACS Paragon Plus Environment

22

Page 23 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Conclusions Koopmans’ theorem has been used for a long time for interpreting photoelectron spectra of molecules and solids, but it is also used for simulating e.g. electrochemical quantities58,59 in MD simulations. Here we addressed the question whether a modified Koopmans’ theorem based on OT-RSH functionals can be used in conjunction with PIMD simulations to quantitatively model liquid photoemission spectra. The results obtained here for water seem to be rather promising in this respect. The main finding of the present study can be summarized as follows: (i) The structure of the liquid photoemission spectrum is well captured within the OT-RSH scheme combined with the reflection principle method. (ii) Including nuclear quantum effects leads to the broadening of the photoelectron peaks. (iii) Small isotope effect consistent with the experimental results was observed in our calculations. (iv) PE spectra of water clusters converge rapidly to the spectra of liquid water. The calculated PE spectra are in good agreement with the experimental data. While the cluster arrangement used in our study is not fully representative of the liquid structure, good agreement between the present simulations and experiment points to a leading role of the local hydrogen bond structure on the shape of the spectrum. Capping the cluster with QM/MM interface did not change the shape of the spectrum. However, the ionization energies were somewhat shifted; more realistic model within the QM/MM scheme would require employing polarizable QM/MM model. The simulations could also be executed with ab initio MD using periodic boundary conditions. The agreement would also improve if the electronic factors (i.e. transition dipole moments for ionization of electrons from different orbitals) are included into simulations. While these

ACS Paragon Plus Environment

23

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 34

quantities were not calculated in this work, it is in principle possible to calculate them to make the protocol fully non-empirical. The simple solution to account for different ionization probabilities of electrons from different orbitals is based on the so called Gelius model.60 Here, the ionization probabilities are constructed from atomic contributions based on population analysis of the molecular orbitals. Such model is not too accurate, but it can be quite useful for modeling more complex liquids containing heavier elements. The Koopmans’ theorem has some other advantages which are planned to be investigated in the future. The possibility of spatial projection of the calculated spectra allows for investigation of bulk liquids as well as liquid/air interfaces. This is interesting from several perspectives. First, liquid photoemission spectroscopy is in principle a surface sensitive technique and it would be useful to know to what extent the measured spectra depend on the distance from the surface. Second, there are some indications on different redox properties of molecules in bulk and on the interfaces61 which can have implications, e.g. for redox processes on atmospheric aerosols. The Koopmans’ based protocols can also easily reveal the structure of molecules at the interface between the liquid and solid phase (e.g. electrode) and can facilitate detailed understanding of chemical reactions taking place at the interface. Recent advances in photoelectron spectroscopy have also made it possible to study such chemical reactions on interfaces in situ experimentally.5,62 Combining both theoretical and experimental approaches is central for rational engineering of more efficient photoelectrochemical cell design or for energy-conserving materials.

ASSOCIATED CONTENT Supporting Information. Additional details on the theoretical photoemission spectra obtained within classical, PI+GLE molecular dynamics and QM/MM scheme are provided in the

ACS Paragon Plus Environment

24

Page 25 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION Corresponding Author *[email protected] Author Contributions The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Funding Sources This work was supported by the Czech Science Foundation GA CR project No. 13-34168S and MSMT project KONTAKT II LH15081. DH acknowledges financial support from specific university research MSMT No 20-SVV/2016. ACKNOWLEDGMENT We would like to thank to Dr. Bernd Winter and Dr. Robert Seidel (Helmholtz Zentrum Berlin) for providing unpublished data and for valuable discussions. ABBREVIATIONS PIMD, OT-RSH, PI+GLE, PE, PES, RP-PIMD, VIE, GLE, MD REFERENCES (1)

Hüfner, S. Photoelectron Spectroscopy, 3rd revised and enlarged ed. 2003 Springer:

Berlin; New York, 2003.

ACS Paragon Plus Environment

25

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(2)

Page 26 of 34

Winter, B.; Faubel, M. Photoemission from Liquid Aqueous Solutions. Chem. Rev. 2006,

106, 1176–1211. (3)

Šištík, L.; Ončák, M.; Slavíček, P. Simulations of Photoemission and Equilibrium Redox

Processes of Ionic Liquids: The Role of Ion Pairing and Long-Range Polarization. Phys. Chem. Chem. Phys. 2011, 13, 11998–12007. (4)

Seidel, R.; Thürmer, S.; Winter, B. Photoelectron Spectroscopy Meets Aqueous Solution:

Studies from a Vacuum Liquid Microjet. J. Phys. Chem. Lett. 2011, 2, 633–641. (5)

Starr, D. E.; Liu, Z.; Hävecker, M.; Knop-Gericke, A.; Bluhm, H. Investigation of

Solid/Vapor Interfaces Using Ambient Pressure X-Ray Photoelectron Spectroscopy. Chem. Soc. Rev. 2013, 42, 5833–5857. (6)

Banna, M. S.; McQuaide, B. H.; Malutzki, R.; Schmidt, V. The Photoelectron Spectrum

of Water in the 30 to 140 eV Photon Energy Range. J. Chem. Phys. 1986, 84, 4739–4744. (7)

Barth, S.; Ončák, M.; Ulrich, V.; Mucke, M.; Lischke, T.; Slavíček, P.; Hergenhahn, U.

Valence Ionization of Water Clusters: From Isolated Molecules to Bulk. J. Phys. Chem. A 2009, 113, 13519–13527. (8)

Winter, B.; Weber, R.; Widdra, W.; Dittmar, M.; Faubel, M.; Hertel, I. V. Full Valence

Band Photoemission from Liquid Water Using EUV Synchrotron Radiation. J. Phys. Chem. A 2004, 108, 2625–2632. (9)

Thürmer, S.; Ončák, M.; Ottosson, N.; Seidel, R.; Hergenhahn, U.; Bradforth, S. E.;

Slavíček, P.; Winter, B. On the Nature and Origin of Dicationic, Charge-Separated Species Formed in Liquid Water on X-Ray Irradiation. Nat. Chem. 2013, 5, 590–596.

ACS Paragon Plus Environment

26

Page 27 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(10) Gaiduk, A. P.; Govoni, M.; Seidel, R.; Skone, J. H.; Winter, B.; Galli, G. Photoelectron Spectra of Aqueous Solutions from First Principles. J. Am. Chem. Soc. 2016, 138, 6912–6915. (11) do Couto, P. C.; Estácio, S. G.; Cabral, B. C. The Kohn-Sham Density of States and Band Gap of Water: From Small Clusters to Liquid Water. J. Chem. Phys. 2005, 123, 54510. (12) Segarra-Martí, J.; Merchán, M.; Roca-Sanjuán, D. Ab Initio Determination of the Ionization Potentials of Water Clusters (H2O)n (n = 2−6). J. Chem. Phys. 2012, 136, 244306. (13) 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, 34109. (14) Svoboda, O.; Ončák, M.; Slavíček, P. Simulations of Light Induced Processes in Water Based on Ab Initio Path Integrals Molecular Dynamics. II. Photoionization. J. Chem. Phys. 2011, 135, 154302. (15) Coe, J. V.; Earhart, A. D.; Cohen, M. H.; Hoffman, G. J.; Sarkas, H. W.; Bowen, K. H. Using Cluster Studies to Approach the Electronic Structure of Bulk Water: Reassessing the Vacuum Level, Conduction Band Edge, and Band Gap of Water. J. Chem. Phys. 1997, 107, 6023–6031. (16) Öhrwall, G.; Fink, R. F.; Tchaplyguine, M.; Ojamäe, L.; Lundwall, M.; Marinho, R. R. T.; Brito, A. N. de; Sorensen, S. L.; Gisselbrecht, M.; Feifel, R.; Rander, T.; Lindblad, A.; Schulz, J.; Sæthre, L. J.; Mårtensson, N.; Svensson, S.; Björneholm, O. The Electronic Structure of Free Water Clusters Probed by Auger Electron Spectroscopy. J. Chem. Phys. 2005, 123, 54310.

ACS Paragon Plus Environment

27

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 34

(17) Lee, S. Y.; Brown, R. C.; Heller, E. J. Multidimensional Reflection Approximation: Application to the Photodissociation of Polyatomics. J. Phys. Chem. 1983, 87, 2045–2053. (18) do Couto, P. C.; Hollas, D.; Slavíček, P. On the Performance of Optimally Tuned RangeSeparated Hybrid Functionals for X-Ray Absorption Modeling. J. Chem. Theory Comput. 2015, 11, 3234–3244. (19) Ončák, M.; Šištík, L.; Slavíček, P. Can Theory Quantitatively Model Stratospheric Photolysis? Ab Initio Estimate of Absolute Absorption Cross Sections of ClOOCl. J. Chem. Phys. 2010, 133, 174303. (20) Law, Y. K.; Hassanali, A. A. Role of Quantum Vibrations on the Structural, Electronic, and Optical Properties of 9-Methylguanine. J. Phys. Chem. A 2015, 119, 10816–10827. (21) Kaczmarek, A.; Shiga, M.; Marx, D. Quantum Effects on Vibrational and Electronic Spectra of Hydrazine Studied by “On-the-Fly” Ab Initio Ring Polymer Molecular Dynamics. J. Phys. Chem. A 2009, 113, 1985–1994. (22) Crespo-Otero, R.; Barbatti, M. Spectrum Simulation and Decomposition with Nuclear Ensemble: Formal Derivation and Application to Benzene, Furan and 2-Phenylfuran. Theor. Chem. Acc. 2012, 131, 1–14. (23) Svoboda, O.; Ončák, M.; Slavíček, P. Simulations of Light Induced Processes in Water Based on Ab Initio Path Integrals Molecular Dynamics. I. Photoabsorption. J. Chem. Phys. 2011, 135, 154301.

ACS Paragon Plus Environment

28

Page 29 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(24) Schwartz, C. P.; Uejio, J. S.; Saykally, R. J.; Prendergast, D. On the Importance of Nuclear Quantum Motions in near Edge X-Ray Absorption Fine Structure Spectroscopy of Molecules. J. Chem. Phys. 2009, 130, 184109. (25) Della Sala, F.; Rousseau, R.; Görling, A.; Marx, D. Quantum and Thermal Fluctuation Effects on the Photoabsorption Spectra of Clusters. Phys. Rev. Lett. 2004, 92, 183401. (26) Ceriotti, M.; Fang, W.; Kusalik, P. G.; McKenzie, R. H.; Michaelides, A.; Morales, M. A.; Markland, T. E. Nuclear Quantum Effects in Water and Aqueous Systems: Experiment, Theory, and Current Challenges. Chem. Rev. 2016, 116, 7529–7550. (27) Ceriotti, M.; Manolopoulos, D. E.; Parrinello, M. Accelerating the Convergence of Path Integral Dynamics with a Generalized Langevin Equation. J. Chem. Phys. 2011, 134, 84104. (28) Maršálek, O.; Markland, T. E. Ab Initio Molecular Dynamics with Nuclear Quantum Effects at Classical Cost: Ring Polymer Contraction for Density Functional Theory. J. Chem. Phys. 2016, 144, 54112. (29) John, C.; Spura, T.; Habershon, S.; Kühne, T. D. Quantum Ring-Polymer Contraction Method: Including Nuclear Quantum Effects at No Additional Computational Cost in Comparison to Ab Initio Molecular Dynamics. Phys. Rev. E 2016, 93, 43305. (30) Ceriotti, M.; Bussi, G.; Parrinello, M. Colored-Noise Thermostats à la Carte. J. Chem. Theory Comput. 2010, 6, 1170–1180. (31) Dolgounitcheva, O.; Zakrzewski, V. G.; Ortiz, J. V. Electron Propagator Calculations on Uracil and Adenine Ionization Energies. Int. J. Quantum Chem. 2000, 80, 831–835.

ACS Paragon Plus Environment

29

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 34

(32) Roca-Sanjuán, D.; Rubio, M.; Merchán, M.; Serrano-Andrés, L. Ab Initio Determination of the Ionization Potentials of DNA and RNA Nucleobases. J. Chem. Phys. 2006, 125, 84302. (33) Ghosh, D.; Roy, A.; Seidel, R.; Winter, B.; Bradforth, S.; Krylov, A. I. First-Principle Protocol for Calculating Ionization Energies and Redox Potentials of Solvated Molecules and Ions: Theory and Application to Aqueous Phenol and Phenolate. J. Phys. Chem. B 2012, 116, 7269–7280. (34) Govoni, M.; Galli, G. Large Scale GW Calculations. J. Chem. Theory Comput. 2015, 11, 2680–2696. (35) Fang, C.; Li, W.-F.; Koster, R. S.; Klimeš, J.; Blaaderen, A. van; Huis, M. A. van. The Accurate Calculation of the Band Gap of Liquid Water by Means of GW Corrections Applied to Plane-Wave Density Functional Theory Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2014, 17, 365–375. (36) Baer, R.; Livshits, E.; Salzner, U. Tuned Range-Separated Hybrids in Density Functional Theory. Annu. Rev. Phys. Chem. 2010, 61, 85–109. (37) Salzner, U.; Baer, R. Koopmans’ Springs to Life. J. Chem. Phys. 2009, 131, 231101. (38) Kronik, L.; Stein, T.; Refaely-Abramson, S.; Baer, R. Excitation Gaps of Finite-Sized Systems from Optimally Tuned Range-Separated Hybrid Functionals. J. Chem. Theory Comput. 2012, 8, 1515–1531. (39) Stein, T.; Kronik, L.; Baer, R. Reliable Prediction of Charge Transfer Excitations in Molecular Complexes Using Time-Dependent Density Functional Theory. J. Am. Chem. Soc. 2009, 131, 2818–2820.

ACS Paragon Plus Environment

30

Page 31 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(40) Egger, D. A.; Weissman, S.; Refaely-Abramson, S.; Sharifzadeh, S.; Dauth, M.; Baer, R.; Kümmel, S.; Neaton, J. B.; Zojer, E.; Kronik, L. Outer-Valence Electron Spectra of Prototypical Aromatic Heterocycles from an Optimally Tuned Range-Separated Hybrid Functional. J. Chem. Theory Comput. 2014, 10, 1934–1952. (41) Ufimtsev, I. S.; Martínez, T. J. Quantum Chemistry on Graphical Processing Units. 3. Analytical Energy Gradients, Geometry Optimization, and First Principles Molecular Dynamics. J. Chem. Theory Comput. 2009, 5, 2619–2628. (42) Henderson, T. M.; Janesko, B. G.; Scuseria, G. E. Generalized Gradient Approximation Model Exchange Holes for Range-Separated Hybrids. J. Chem. Phys. 2008, 128, 194105. (43) Slavíček, P.; Ončák, Milan; Hollas, D.; Svoboda, O. ABIN, version 1.0, https://github.com/PHOTOX/ABIN (accessed 8/29/2016) (44) GLE4MD

Project:

Automatic

Input

Generation.

http://epfl

cosmo.github.io/gle4md/index.html?page=matrix (accessed 8/29/2016) (45) The Cambridge Cluster Database. http://www-wales.ch.cam.ac.uk/CCD.html (accessed 8/29/2016) (46) Habershon, S.; Markland, T. E.; Manolopoulos, D. E. Competing Quantum Effects in the Dynamics of a Flexible Water Model. J. Chem. Phys. 2009, 131, 24501. (47) Sankari, R.; Ehara, M.; Nakatsuji, H.; Senba, Y.; Hosokawa, K.; Yoshida, H.; De Fanis, A.; Tamenori, Y.; Aksela, S.; Ueda, K. Vibrationally Resolved O 1s Photoelectron Spectrum of Water. Chem. Phys. Lett. 2003, 380, 647–653.

ACS Paragon Plus Environment

31

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 32 of 34

(48) We thank to Dr. B. Winter and Dr. R. Seidel for providing us the unpublished experimental data. These data are very similar to those published in ref. 53. (49) The secondary electrons were subtracted by Dr. R. Seidel. (50) Cederbaum, L. S.; Domcke, W.; Schirmer, J.; Niessen, W. V. Correlation Effects in the Ionization of Molecules: Breakdown of the Molecular Orbital Picture. In Advances in Chemical Physics; Prigogine, I., Rice, S. A., Eds.; John Wiley & Sons, Inc., 1986; pp 115–159. (51) Jahnke, T. Interatomic and Intermolecular Coulombic Decay: The Coming of Age Story. J. Phys. B At. Mol. Phys. 2015, 48, 82001. (52) Grell, G.; Bokarev, S. I.; Winter, B.; Seidel, R.; Aziz, E. F.; Aziz, S. G.; Kühn, O. MultiReference Approach to the Calculation of Photoelectron Spectra Including Spin-Orbit Coupling. J. Chem. Phys. 2015, 143, 74104. (53) Slavíček, P.; Kryzhevoi, N. V.; Aziz, E. F.; Winter, B. Relaxation Processes in Aqueous Systems upon X-Ray Ionization: Entanglement of Electronic and Nuclear Dynamics. J. Phys. Chem. Lett. 2016, 7, 234–243. (54) Prakash, M. K.; Weibel, J. D.; Marcus, R. A. Isotopomer Fractionation in the UV Photolysis of N2O: Comparison of Theory and Experiment. J. Geophys. Res. Atmospheres 2005, 110, D21315. (55) Isborn, C. M.; Mar, B. D.; Curchod, B. F. E.; Tavernelli, I.; Martínez, T. J. The Charge Transfer Problem in Density Functional Theory Calculations of Aqueously Solvated Molecules. J. Phys. Chem. B 2013, 117, 12189–12201.

ACS Paragon Plus Environment

32

Page 33 of 34

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

(56) Ellison, W. J.; Lamkaouchi, K.; Moreau, J.-M. Water: A Dielectric Reference. J. Mol. Liq. 1996, 68, 171–279. (57) Sanner, M. F.; Olson, A. J.; Spehner, J.-C. Reduced Surface: An Efficient Way to Compute Molecular Surfaces. Biopolymers 1996, 38, 305–320. (58) Cheng, J.; Liu, X.; VandeVondele, J.; Sulpizi, M.; Sprik, M. Redox Potentials and Acidity Constants from Density Functional Theory Based Molecular Dynamics. Acc. Chem. Res. 2014, 47, 3522–3529. (59) Adriaanse, C.; Cheng, J.; Chau, V.; Sulpizi, M.; VandeVondele, J.; Sprik, M. Aqueous Redox Chemistry and the Electronic Band Structure of Liquid Water. J. Phys. Chem. Lett. 2012, 3, 3411–3415. (60) Gelius, U.; Siegbahn, K. ESCA Studies of Molecular Core and Valence Levels in the Gas Phase. Faraday Discuss. Chem. Soc. 1972, 54, 257–268. (61) Martins-Costa, M. T. C.; Anglada, J. M.; Francisco, J. S.; Ruiz-Lopez, M. F. Reactivity of Atmospherically Relevant Small Radicals at the Air–Water Interface. Angew. Chem. Int. Ed. 2012, 51, 5413–5417. (62) Wu, C. H.; Weatherup, R. S.; Salmeron, M. B. Probing Electrode/Electrolyte Interfaces in Situ by X-Ray Spectroscopies: Old Methods, New Tricks. Phys. Chem. Chem. Phys. 2015, 17, 30229–30239.

ACS Paragon Plus Environment

33

Journal of Chemical Theory and Computation

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 34 of 34

Table of Contents Graphic

ACS Paragon Plus Environment

34