Optimal Tuning of Range-Separated Hybrids for Solvated Molecules

Article pubs.acs.org/JCTC

Optimal Tuning of Range-Separated Hybrids for Solvated Molecules with Time-Dependent Density Functional Theory Martina Rubešová, Eva Muchová, and Petr Slavíček* Department of Physical Chemistry, University of Chemistry and Technology, Technická 5, 16628 Prague, Czech Republic S Supporting Information *

ABSTRACT: The applicability range of density functional theory (DFT) can be improved with no additional parametrization by imposing some exact conditions. Enforcing equality between the orbital energy of the highest occupied Kohn−Sham orbital and ionization energy, determined from the total energy difference between neutral and ionized states (ΔKS), leads to the concept of optimally tuned range-separated hybrid functionals. Here, we present an alternative tuning scheme for range-separated hybrid functionals based on enforcing the equality between the ΔKS ionization energy and the ionization energy calculated by means of the time-dependent DFT using the concept of ionization as an excitation to the distant center (OT-IEDC scheme). The scheme can be naturally applied to solvated systems described either within the explicit solvation or dielectric continuum models. We test the performance of the scheme on a benchmark set of molecules. We further show that the scheme allows for reliably modeling liquid phase photoemission spectra.

■

INTRODUCTION Recent advances in liquid state photoemission spectroscopy1,2 pose challenges to theoretical models for ionization energies in the liquid state. To interpret the spectra, theoretical chemistry needs to quantitatively estimate binding energies for each electron in a system. This goal can be accomplished in a thrifty way within the Hartree−Fock (HF) method via the Koopmans’ theorem,3 equaling the ionization energy (IE) of the ith electron to the negative of the respective orbital energy(ε) IEi = −εi

Kohn−Sham density functional theory (KS-DFT) is an effective one-particle theory that includes electron correlation, and its use is typically considered as the most efficient improvement upon the HF theory. It is therefore natural to think of KS orbital energies in Koopmans’ terms. The ionization potential theorem for the HOMO orbital IE = E(N − 1) − E(N ) = −εHOMO

holds exactly for the KS-DFT with an exact energy functional.13−16 Chong et al.4 showed that the exact KS orbital energies are close to exact IEs even for outer-valence orbitals; the eigenvalues close to εHOMO deviated typically by 0.05−0.1 eV. For approximate functionals, the condition in eq 2 is, however, severely violated.17,18 The generalized gradient approximation (GGA) functionals tend to shift the HOMO energy by several eVs;4 hybrid functionals with a fraction of exact HF exchange perform slightly better (for example, the B3LYP HOMO energy is typically 70% of the experimental IEs17). Imposing the condition in eq 2 can be conveniently enforced for range-separated hybrid functionals (RSHs)19 in which the amount of exact exchange smoothly varies based on an inter electron distance r12

(1)

The quantitative performance of Koopmans’ theorem is, however, poor. It yields inaccurate IEs for valence electrons, and the error can become as large as tens of eVs for core electrons.4,5 There are two principal error sources: neglecting the correlation energy and, dominantly, ignoring orbital relaxation upon ionization. More accurate ionization energies can be obtained with advanced techniques, for instance, using Green function-based methods, e.g., the outer valence Green function approach6 or the so-called GW (Green function G with Coulomb interaction W) calculations.7 High quality/high cost approach for ionization energies can be achieved via coupled cluster-based approaches such as the equation of motion for ionization potential coupled cluster method with single and double excitations (EOM-IP-CCSD),8,9 symmetryadapted cluster/configuration interaction (SAC−CI),10 secondorder multireference perturbation theory (CASPT2),11 or algebraic diagrammatic construction (ADC(2)) methods.12 Koopmans’ theorem remains attractive, however, due to its conceptual simplicity and computational feasibility. © 2017 American Chemical Society

(2)

erf(ωr12) erfc(ωr12) 1 = + r12 r12 r12

(3)

The ratio between the nonlocal exact HF exchange (applied to the first term) and the semilocal exchange term (applied to Received: June 27, 2017 Published: September 5, 2017 4972

DOI: 10.1021/acs.jctc.7b00675 J. Chem. Theory Comput. 2017, 13, 4972−4983

Article

Journal of Chemical Theory and Computation

Figure 1. Two schemes for calculating ionization energies within TDDFT. (A) Calculation of ionization energies by summing the lowest ionization energy IEHOMO and excitation energies to a singly occupied molecular orbital (SOMO, depicted as a blue circle) denoted as ΔEi→k. Excitations in which the SOMO is not refilled with an electron, as those depicted in a lower part of panel A, should be projected out, i.e., we should avoid the satellite states that typically represent only a minor contribution to the photoemission spectrum. (B) Ionization as an excitation into a distant center (IEDC) approach.

into the energetics of the solute molecule. For extended systems, ω values can be very low, e.g., tuning drives the rangeseparated functional toward the corresponding GGA functional, which provides too small HOMO−LUMO gaps.18,41,47−49 The same effect arises for a solute embedded in a dielectric continuum. As a result, Boruah et al.43 concluded that “The performance of range-separated functionals in reproducing the HOMO and LUMO energies in the solvent medium is not encouraging.” Very recently, Zheng et al.50 presented a novel approach for the optimal tuning scheme in a polarizable continuum based on inclusion of the solvent screening effect directly into the optimally tuned range-separated hybrid functional (the so-called screened range-separated hybrid (SRSH) functional51). In the SRSH scheme, the dielectric constant enters directly the partitioning of the 1 operator,

the second term) is governed by the range-separation parameter ω. The division above applies for the long-range corrected functionals. Note that, for some flavors of RSH, the exact exchange might apply also for a short range; throughout the paper, we will further consider here only the long-range corrected RSHs. The RSHs have greatly improved the performance of DFT for thermochemistry, molecular geometries, barrier heights, and excited states with charge transfer character.20−26 We can impose the condition in eq 2 to be fulfilled by tuning ω to satisfy the ionization potential theorem IEω = E ω(N − 1) − E ω(N ) = −εHOMO(ω)

(4)

where Eω(N − 1) refers to the energy of the ionized system and Eω(N) to the energy of the neutral system computed with the RSH functional with the corresponding system-tuned rangeseparation parameter ω; their difference is denoted as ΔKS. Although ΔKS only weakly varies with ω, KS orbital energies change significantly. Enforcing eq 4 was suggested by Stein, Kronik, and Baer as an efficient route to treat charge transfer excitation.27 Salzner and Baer suggested28 that enforcing eq 4 guarantees the validity of eq 1 (Koopmans’-like theorem) even for lower lying orbital energies, see also ref 26. In this case, the KS orbital energies also include orbital relaxation (by virtue of eq 4) and correlation (to the extent to which the correlation energy is described by the correlation functional). Optimal tuning of the range-separated hybrid functionals (OT-RSHs) has since become a widely used procedure for obtaining ionization characteristics of molecular systems.17,29−31 It has also been shown that tuning of the RSH functionals helps tailoring DFT for calculations of absorption energies,17,30,32−36 including the X-ray absorption spectra.37−40 A question pertaining to the performance of the OT-RSH scheme for solvated molecules will now be addressed. Previously, the applicability of the OT-RSH concept in solution was studied by de Queiroz et al.,41,42 Boruah et al.,43 and Sun et al.;44 several other studies were based mainly on tuning RSH functionals in vacuum and a posteriori including solvent shifts to ionization energies.45,46 As has been shown by de Queiroz et al., in solution, the ionization energy decreases due to solvent electronic polarization while the orbital energy remains approximately constant for the given ω parameter.41 Because the orbital energy monotonically decreases with ω, fulfilling eq 4 yields artificially low ω values. In another words, the exact exchange part of the functional is suppressed even though the electronic structure of the solute system does not really change. Clearly, the electronic relaxation of the environment is packed

r12

whereas the ω parameter is fixed to its gas-phase value. We present below a new approach for modeling ionization energies in solution with the range-separated hybrid functionals using an alternative tuning scheme. We suggest a tuning scheme for range-separated hybrid functionals working consistently for molecules both in the gas phase and in the solution. The scheme is based on formulating the exact condition not in terms of orbital energies but via well-defined ionization energies based on time-dependent density functional theory (TDDFT52) with restricted occupations. This approach represents a suitable basis for an extension of the Koopmans’ theorem within the DFT. Furthermore, dielectric continuum models can easily be incorporated. Finally, even higher ionization energies are accurately described. We provide details of the proposed protocol in the next section. Theoretical Approach. Our general aim is to tune RSH functionals so that they provide a reliable electronic structure of molecules both in the gas phase and in solution. The quality of the electronic structure is tested by calculating a complete spectrum of ionization energies, or a significant part of the spectrum (e.g., valence or core electron part), on equal theoretical footing. We use TDDFT to calculate the electron binding energies in a system. One possible approach for calculating ionization energies with TDDFT is depicted in Figure 1A. In this approach, HOMO ionization energy is calculated via ΔKS (IEHOMO), and subsequently, higher ionization energies are estimated as IEi = IE HOMO + ΔEi → k

(5)

where ΔE is the excitation energy of the radical cation formed upon the HOMO electron ionization.53−57 This i→k

4973


Article


energies, which is not fully guaranteed, especially for large systems.66 The scheme can be easily adapted to include the effect of the environment; the critical point is to properly account for the screening effect of the environment on the excitation energies. We employ the concept of nonequilibrium solvation, which has been developed and widely tested for excitation processes.67−72 Both within the ionization (left-hand side of eq 7) and charge transfer excitation (right-hand side of eq 7), the electronic structure changes take place on a subfemtosecond time scale. The solvent responds to a new charge distribution in the solute, yet the solvent’s microscopic motions have different characteristic times to reach a polarization response. Although the response of the solvent’s electrons is immediate (it takes place simultaneously with the excitation or ionization), the solvent nuclei relax at a much slower rate fully within the spirit of the Born−Oppenheimer approximation. The ionized/excited solute appears in a nonequilibrium environment at the moment of the ionization/excitation. This effect is captured within a dielectric model provided that we know the frequency dependence of the dielectric constant. For the nonequilibrium effects to be included, the reaction field is separated into two parts. The electronic part of the solvent response is modeled by the high frequency limit of the dielectric constant, e.g., instead of the full dielectric constant (80.4 for water at 25 °C), only the high-frequency or optical part (εopt = 1.78 for water) applies for the final state.73,74 The slow component of ε is taken as ε − εopt, where ε is the static dielectric constant. By the slow/fast division, we consider consistently the vertical character of the excitation/ionization for both the solute and solvent. Two schemes to account for the solvation energetics in the excited states have been designed: state-specific (SS)74−77 and linear response (LR) schemes.76,78,79 In our simulations, we adopted the SS scheme, which is more suitable for processes involving large reorganizations of electronic densities.69,80 In the SS approach, it is assumed that, upon reorganization of the solvent’s electronic structure (excitation or ionization), the “fast” solvent degrees of freedom respond immediately to the new electron density of a solute. A new wave function |ψi⟩ is obtained from a nonlinear Schrödinger equation

approach has certain drawbacks. First, the HOMO ionization and higher ionizations are not treated on equal footing. Second, a radical cation often suffers from spin contamination (albeit less than for wave function methods58−60) and TDDFT excitations for open-shell systems are generally less reliable than the respective calculations for ground state closed-shell systems. Finally, not all excited states correspond to single electron ionizations, the TDDFT excitations lead also to satellite states of the two hole-single particle type (see panel A of Figure 1). These states have to be removed by a posteriori projection based on the inspection of the character of the transitions. Here, we present a second option in which we treat ionization as an excitation into a distant center (IEDC, see Figure 1B). A similar approach was previously used for wave function-based methods.61 We employed this approach previously in the context of core-level spectroscopy of iron complexes.62 We provide a more detailed exploration, profound testing of the method, and discussion of its applicability in the context of solvation modeling. Let us assume a target system with ionization energy >5 eV. We can place a sodium cation (with electron affinity around 5 eV63), for example, at a large distance from the investigated molecule. Using the TDDFT method, we calculate the charge transfer excitation from the target system to the sodium cation. The ionization energy is then given as +

i → 3s(Na ) IE i = ΔECT + EA Na+ + ΔEcorr

(6)

+

) where ΔEi→3s(Na is the charge transfer excitation from the CT molecule to the sodium cation, calculated at the TDDFT level, EANa+ is the electron affinity of the sodium cation, and ΔEcorr is assumed to be a constant correction term to be discussed shortly. The sodium cation is employed just for convenience; in fact, we can use an arbitrary electron-accepting center: either a real molecular system or a suitably designed artificial system with a required electron affinity. The electron affinity of the sodium cation must be calculated for each particular method and basis set. Here, the method was presented for TDDFT, but other methods like EOM-CCSD can be used as well. The IEDC method then provides a consistent set of ionization energies as all the excitations are treated equally. Moreover, with certain restrictions of the orbital space, IEDC gives a whole orbital energy spectrum. The obvious deficiency of the approach is connected mainly to the TDDFT description of charge transfer states. It is well-known that approximate functionals underestimate the energies of the CT states (spurious charge transfer states).52,64,65 One option to overcome the problem is to include the ΔEcorr correction factor, which hopefully has a constant value for various systems or at least for various states. An alternative approach is to employ the optimal tuning concept; in fact, the IECD method has a big potential to improve the optimal tuning (OT-IECD). We can set the ω parameter so that the equation

slow fast SS vac Ĥ i |ψi ⟩(Ĥ + V0̂ + Vî )|ψi ⟩ = EiSS|ψi ⟩

where Ĥ vac is a standard Hamiltonian for a solute in vacuum. A reaction field operator is separated into two parts denoted as ̂ fast V̂ slow 0 (i = 0 represents the ground state) and Vi , e.g., the reaction field operator creates two electrostatic potentials of the apparent surface charge densities that describe fast and slow response of the solvent. To solve eq 8, the polarization of a solute should be solved self-consistently with respect to the solvent’s reaction field. This nonequilibrium solvent description is straightforward for the ΔKS calculations (the left-hand side of eq 7).53,81−85 The situation is a little more complex for the nonequilibrium vertical excitations67,80 within the TDDFT, partially because of nonorthogonality of the wave functions. To circumvent the problems, one can calculate the excitation energies within a zero-order perturbative approximation to SS67,68,80 with the ptSS method. In this approach, first order corrections to the ith state energy, polarization energy, or vertical excitation energy can be expressed with the zerothorder excited state wave function. This wave function is constructed from relaxed molecular orbitals for the ground state

+

HOMO → 3s(Na ) E ω(N − 1) − E ω(N ) = ΔECT + EA Na+

(8)

(7)

is either satisfied or the difference between the left- and righthand side of the equation is minimized. We argue below that this type of tuning is more stable and consistent compared to tuning per eq 4, which is formulated for the HOMO orbital energy. Note, however, that the approach still relies on the ability of the DFT method to calculate the ΔKS ionization 4974


Article


photoemission process. In the photoemission experiment, the intensity of emitted electrons is recorded as a function of electron binding energies (or alternatively, kinetic energies of ejected electrons) for a constant photon energy. Considering constant oscillatory strengths, we record the abundance of binding energies within a certain energy range along the PIMD trajectory. The PIMD simulations cover both thermal and nuclear quantum effects. The PIMD simulations were accelerated with the so-called quantum thermostat based on the generalized Langevin equation (GLE) formalism (PI + GLE scheme).94,95 The GLE thermostat parameters were taken from an online database.96 The total length of the PI + GLE simulations was 50 ps. The temperature in all simulations was set to 300 K, and the time step was set to 20 au (∼0.5 fs). For ∼500 snapshots obtained from the PIMD (more specifically PI + GLE) simulations, we calculated the ionization energies by the IEDC approach and performed energy binning. Solvent spectral broadening was also modeled within the reflection principle with an additional broadening (RP-AB) scheme.97 In the AB approach, each ith snapshot is characterized by a reorganization energy λi given as a difference between the energy of the ionized system using the nonequilibrium (Eion,noneq) and equilibrium solvation (Eion,eq)

in the PCM. The Q-Chem 4.1 code used in the present work includes additional corrections via the perturbative linearresponse terms (ptLR), and the excitation energies are in good agreement with the experimental data.67,80 Technical Details: Benchmark Data Set and Computational Methods. We tested the approach on a set of molecules for which accurate values of electron binding energies in solution are available. The test set included anions, neutral solutes, and a cation considering both organic molecules and inorganic systems, particularly the imidazolium cation, imidazole, adenine, phenol, aniline, veratrole, dimethyl sulfide (DMS), the water molecule, hydrogen peroxide, the phenolate anion, and the phosphate anions. For selected molecular systems (phenol, uridine, and water clusters), we modeled the whole valence photoemission spectrum. Phenol spectra comprise two closely lying ionization states, e.g., we could test the performance of the IEDC approach for higher ionizations. Uridine is a multichromophoric system; it contains two chromophore units, and its spectrum includes many ionized states. Water clusters represent an extreme case; here, the system contains a large number of units with almost degenerate ionization energies. For electronic structure calculations, we employed hybrid and RSH functionals with the 6-31+g* basis set. Specifically, we used the BMK functional,86 which contains a higher fraction of exact exchange (42%) and performs well for ionization energies.53 Next, we used the global hybrid Minnesota M06HF functional87 containing 100% of HF exchange. Dominantly, we focused on the RSH functionals where the optimal tuning is possible, here represented by the LC-ωPBE functional.88 Local or GGA functionals were not tested as they are clearly destined to fail for charge transfer excitations. The values were compared to benchmark calculations at the CCSD/aug-cc-pVDZ level. The ionization energies at the respective levels were obtained both as the ΔCCSD values and by the IEDC method with the EOM-CCSD/aug-cc-pVDZ method. We assessed the performance of the IEDC method in the gas phase as well as in water environment. The solvent effects were included via a dielectric continuum model (C-PCM)89−91 with a nonequilibrium model of solvation (NEPCM) for the ionized states, e.g., we used only the optical part of dielectric constant for water (1.78). Because the studied electronic processes involve large reorganization of electron densities, we employed the state specific (SS) scheme as implemented in the Q-Chem 4.1 code.80 To model the photoemission spectra in solution, we have to account for inhomogeneous broadening caused by interactions between a solute and solvent in both the ground and ionized/ excited states. For phenol, phenolate, and water, we simulated the photoemission spectra with a semiclassical reflection principle combined with path integral molecular dynamics (PIMD) simulations,29,92,93 i.e., we calculated the electron binding energies along a PIMD trajectory of the isolated molecule

λi = E ion,noneq (R⃗ i) − E ion,eq (R⃗ i)

The reorganization energy classically relates to a spectral variance Si as Si2 = 2kBTλi

(11)

where kB is Boltzmann constant and T is temperature. The photoemission cross section is then given as a sum (Nsnapshot) (Nelectrons)

σ (E ) ≈

∑

∑

i=1

n=1

2 2 1 e−(E − VIEi ,n) /2Si Si 2π

(12)

where Nsnapshot denotes the number of snapshots taken from the molecular simulation and Nelectrons denotes the number of ionized electrons. VIEi,n stands for the vertical ionization energy of the ith snapshot and nth electron. In other words, in the RPAB scheme, each snapshot is broadened with a Gaussian function with a variance based on electrostatic arguments. In the extreme case, we can take only a single structure that is empirically broadened. This approach was used for uridine using a variance of 0.34 eV. Phenol, phenolate, and water were treated using the RP-AB scheme with the initial ground state nuclear density sampled with the PI + GLE simulations. Most electronic structure calculations were performed either in the Q-Chem 4.1 package,98 the CCSD, M06, and BMK calculations were performed in the Gaussian09 code, revision D.01.99 Initial sampling of the ground state densities with the PI + GLE scheme was performed in our in-house build program ABIN.100

■

(Nelectrons)

σ (E B ; E , E k ) ≈

(10)

RESULTS AND DISCUSSION The results section is divided into three parts. First, we show the performance of the IEDC approach on a benchmark set of molecules; we present the HOMO ionization energies (or HOMO−n ionization energies, if the experimental data are available) for molecules in the gas phase and in solution. Second, we model the valence photoemission spectra of phenol and uridine. Finally, we present a computationally difficult case

∑ ∫ ρ(R⃗)fn (E)δ(E − EB,n − E k ) dR⃗ n=1

(9)

where ρ(R⃗ ) is the nuclear density evaluated with the PIMD simulations, EB,n is the binding energy of the ejected electron, Ek is the kinetic energy of an ejected electron, E is the incident photon energy, and f n(E) is the oscillatory strength for the 4975


Article

Journal of Chemical Theory and Computation Table 1. Ionization Energies for a Benchmark Set of Moleculesa LC- ωPBE

CCSD experiment

gas phase

diel. cont.

gas phase

diel. cont.

compound

orbital

gas phase

liquid

IEDC

ΔCCSD

ΔCCSD

ω

IEDC

ω

IEDC

imidazolium imidazole

HOMO HOMO HOMO−1 HOMO HOMO HOMO−1 HOMO HOMO−1 HOMO HOMO−1 HOMO HOMO HOMO−1 HOMO−2 HOMO HOMO−1 HOMO HOMO HOMO HOMO

NA 8.78−9.12101 10.3896 8.45−8.49103 8.56−8.75101 9.28104 8.02−8.10101 9.12104 NA NA 8.67−8.71101 12.62106 14.84107 18.76107 11.65108 13.02108 NA NA NA NA

8.9681 8.26,858.51102 9.80102 7.5753b 7.80105 8.60105 7.49102 8.62102 7.68102 8.63102 NA 11.16107 13.50107 17.34107 9.95108 10.67108 7.00−7.20105 8.30−8.5056 8.80−9.0056 9.30−9.7056

14.80 8.88 10.19 8.39 8.63 9.37 8.02 9.14 8.13 9.07 8.50 12.34 14.60 18.80 11.19 12.50

15.10 8.80

9.45 7.84

8.25 8.51

7.50 7.68

7.87

7.08

7.79

7.15

8.51 12.43

7.61 11.03

11.44

10.20

0.372 0.326 0.326 0.295 0.307 0.307 0.300 0.300 0.277 0.277 0.344 0.498 0.498 0.498 0.489 0.489

15.36 9.04 10.12 8.42 8.62 9.55 7.95 9.24 7.81 8.75 8.73 12.70 14.83 18.74 11.86 13.18

0.325 0.352 0.352 0.324 0.306 0.306 0.330 0.330 0.312 0.312 0.376 0.512 0.512 0.512 0.509 0.509 0.368 0.334 0.358 0.373

9.41 7.98 9.35 7.60 7.66 8.54 7.07 8.37 7.12 7.99 7.74 11.18 13.34 17.01 10.50 11.78 5.83 6.42 7.51 9.26

adenine phenol aniline veratrole DMS H2O

H2O2 phenolate PO43− HPO42− H2PO4−

5.61 6.03 7.15 8.91

The first column shows experimental ionization energies; then, we present the calculated values using the CCSD method with the aug-cc-pVDZ basis set (using either the IEDC method for ionization energy or the difference between the ground and ionized state ΔCCSD) and optimally-tuned IEDC with LC-ωPBE (6-31+g* basis set) for both the gas phase and dielectric continuum environments. For LC-ωPBE, the ω parameter is tuned so that the condition in eq 7 is satisfied. Blank cells refer to calculations that could not be performed (e.g., the ionization energies are below 5 eV). We did not perform the IEDC calculations within a dielectric continuum as the EOM-CCSD method could not be presently used together with a consistent dielectric continuum model describing the non-equilibrium solvation. bBecause of the low solubility of adenine, we show for comparison the experimental data of adenosine instead. a

between ΔCCSD and OT-IEDC approaches is still very good; however, both values significantly deviate from the experiment, especially for multiply charged ions; this result is not surprising. The polarized continuum models are developed for neutral molecules, and the performance is significantly worse for multiply charged species.56 The IEDC approach clearly cannot exceed the quality of the ΔKS or ΔCCSD calculations. We demonstrate the difference between the newly suggested OT-IEDC scheme and the tuning scheme based on orbital energies on the example of imidazole in Figure 2. Let us start with discussing the tuning according to eq 4. We observe that the ΔKS energies depend only very slightly on the ω parameter, whereas the orbital energies show a rather strong variation; the two curves intersect at around 0.3 a−1 0 , and the value is very close to the optimal range-separation parameter within the OT-IEDC approach (eq 7). The orbital energies, however, do not change significantly upon solvation, whereas ΔKS values change significantly, especially for the equilibrium model of solvation. As a result, the ω parameter incorrectly shifts to very low values (0.05 a−1 0 for equilibrium solvation model, 0.17 a−1 0 for the nonequilibrium solvation model. An analogous drop is observed for all other molecules and follows findings from previously published results by de Quieroz et al.42 Such a large shift leads to a drastic shrinking of the HOMO− LUMO gap upon solvation. By contrast, the IEDC method does not yield significantly different results for the ω parameter in the gas and liquid phases, e.g., for imidazole it only slightly increases from 0.33 to 0.35. This is not surprising because the electronic structure of neutral molecules is not expected to be affected by solvation for neutral molecules. The IEDC tuning

of a system with many almost degenerate ionization states, e.g., photoemission from water clusters. For water clusters, we also show the optimal ω parameter converges with respect to cluster size. Let us start with the performance of IEDC with LC-ωPBE/631+g* (OT-IEDC). The ionization energies are presented in Table 1. We compare the values to experimental data and to the CCSD/aug-cc-pVDZ values obtained both as ΔCCSD approach and by IEDC/EOM-CCSD method. In the gas phase, the experimental results are in a good agreement with ab initio calculations at all levels. We can thus infer that (i) the IEDC concept is reasonable as the IEs obtained at the CCSD level closely match for the ΔCCSD and IEDC approaches and (ii) the optimal tuning with the LC-ωPBE functional provides consistent results. We also observe a close match between the CCSD-IEDC and OT-IEDC results for lower-lying electrons. As for the simulations in the liquid phase, the most relevant comparison is between the ΔCCSD calculation in the dielectric continuum and the OT-IEDC scheme, again considering the nonequilibrium solvation model. The difference here is typically less than 0.3 eV for the neutral and cationic species. The difference is slightly larger when comparing the values with the experimental data, but it is a problem of the dielectric continuum models per se. Isolated anions in the gas phase can be electronically unstable (for instance, doubly and triply charged phosphates), and they are only stabilized in solution (however, a singly charged dihydrogen phosphate anion is stable in the gas phase). We can consistently discuss the performance of the OT-IEDC method for ionization only for the solvated anions. As can be seen in Table 1, the agreement 4976


Article


Figure 2. Dependence of the ionization and orbital energies for imidazole in the gas phase and for imidazole in a dielectric continuum. The graph shows the HOMO orbital energy in the gas phase (brown) and in the dielectric continuum (red), the IEDC energies in the gas phase (gray) and in the dielectric continuum (green), the ΔKS energies in the gas phase (blue), in the dielectric continuum using the nonequilibrium solvation model (black) and in the dielectric continuum using the equilibrium solvation model (purple). The IEDC optimal range-separation parameter ω for the gas phase is 0.33 a−1 0 (intercept of the gray and blue curves) and for the dielectric continuum 0.35 a−1 0 (intercept of the black and green curves). The difference is significantly higher for the optimal tuning using orbital energies where the optimal ω parameter is for the gas phase almost the same as for the IECD scheme (0.33 a−1 0 , intercept of the brown and blue curves) but amounts to only 0.17 a−1 0 in the dielectric continuum (intercept of the red and black curves) with the nonequilibrium solvation (denoted as noneq). The equilibrium solvation model (denoted as eq) provides optimal ω too low for all discussed cases.

Figure 3. Photoemission spectrum of a bare phenol embedded in a dielectric continuum model (from the PIMD simulations) compared with the experimental data.105 The dashed line corresponds to spectra obtained within the additional broadening scheme.97 The C-PCM model was used for modeling of water solvent. The photoemission spectrum was evaluated along the PI + GLE trajectory. Optimal tuning was performed for each geometry with a mean ω value of 0.304.

IEDC method is in reasonable agreement with the experiment. Unlike in the experiment, the peaks corresponding to ionization of different electrons are well resolved. The agreement with the experiment in the fwhms is, however, improved by assuming additional broadening due to fluctuations of the solvent molecules. We fitted the spectra to two Gaussians with peak positions and fwhm presented in Table 2; the peak positions Table 2. Full Width at Half Maximum of the Two Peaks in the Phenol Photoemission Spectruma

clearly avoids artificial lowering of the range-separation parameter. Had we used the ΔKS energy difference within the equilibrium solvation model, the range-separation parameter would have moved to very low values again. The physically correct behavior is achieved by virtue of nonequilibrium solvation accounting for the electronic polarization of the solvent. This can be contrasted with the recently introduced scheme of Zheng et al.50 in which the value of the ω parameter is assumed to be equal for the gas phase and for the solution. The performance of the IECD approach combined with hybrid functionals is significantly worse. The ionization energies obtained for the same set of molecules with the BMK method (42% of the exact exchange) and with the M06-HF method (100% of the exact exchange) are presented in Table S1. In the ideal case, the ΔEcorr term in eq 6 would remain constant for all systems. This is, however, not the case. One could expect that the M06-HF method with a full exact exchange part in particular could provide consistent results, yet this expectation was not met. Apparently, the M06-HF approach suffers from an imbalance between the correlation and exchange part. The optimal tuning thus appears as an essential prerequisite for the application of the IEDC method within the DFT framework. Next, we present the valence photoemission spectra with the OT-IEDC approach. Figure 3 shows the calculated spectrum of solvated phenol compared with the experiment.105 The experimental spectrum can be fitted to two Gaussian functions; the experimental peaks are centered at 7.8 ± 0.1 and 8.6 ± 0.1 eV (at photon energy of 200 eV) with a full width at half maximum (fwhm) of ∼1 eV. The spectrum obtained with the

fwhm [eV]

ionization energy [eV]

peak no.

1st

2nd

1st

2nd

phenol phenol (AB) Pphenol exp.b

0.44 0.73 1.07

0.36 0.68 1.15

7.84 7.84 7.77

8.50 8.50 8.62

a

Comparison of spectra modeled for bare phenol embedded in a dielectric continuum with and without using additional broadening (AB) scheme and the bexperimental data.105 The C-PCM model was used for modeling of water solvent. The photoemission spectrum was evaluated along the PI + GLE trajectory. Optimal tuning was performed for each geometry with a mean ω value of 0.304.

are 7.84 and 8.50 eV. As was previously reported by Ghosh et al.,105 almost no electron density is localized on the surrounding water molecules in the final ionized state of phenol. Microsolvation in such a case does not bring any significant improvement to the calculated spectra concerning their positions; it could, however, improve the agreement of the peak widths.105 The experimental widths are unusually high, which suggests that specific solvent effects play a role. The second photoemission spectrum was calculated for a nucleoside uridine. Photoemission spectrum for this molecule covers several ionization states and poses a bigger challenge for accurate spectra modeling. Here, we calculated the spectrum with empirical broadening of the stick spectrum in minimum energy structures to compare the results with previous calculations.109 The results obtained by the OT-IEDC method with empirical broadening are presented in Figure 4. As can be 4977


Article


Finite size clusters in the gas phase containing 6 and 35 water molecules do not represent the liquid structure, and their spectra are predictably shifted toward higher energies compared with the experimental values from liquid photoemission. A more relevant comparison can be performed for a cluster of 35 water molecules for which a gas phase cluster experiment exists.110 The values agree well, yet the predicted widths are slightly lower. For clusters containing 6 and 35 water molecules embedded in the polarized continuum, the peak positions are in excellent agreement with experimental data,111 and the peaks widths are again somewhat lower. The spectra can be improved in several directions (e.g., including ionization cross sections) as was discussed earlier,29 but these improvements are not connected to the performance of the OT-IEDC method and would be beyond the scope of the present paper. We can again conclude here that the OT-IEDC method predicts accurate results for ionization energies even for such a difficult system. For water clusters, we also studied the dependence of the ω parameter in the range-separated functional LC-ωPBE on the size of water clusters. We performed the calculations either for optimized structures taken from the Cambridge cluster database112 or a set of structures generated by the PIMD simulations (typically using 100 frames). The results are shown in Figure 6. As can be seen, the optimal omega parameter (ω) for the isolated optimized water clusters decreases exponentially with cluster size. An exponential decay has been previously reported by de Queiroz and Kümmel41 for tuning of ω with increasing system size of oligothiophenes both in the gas phase and in a dielectric continuum or by Sosa Vazquez et al.114 for tuning of ω for one ethane molecule surrounded by an increasing number of explicit water molecules. The curve can be fitted as

Figure 4. Photoemission spectrum of uridine embedded in a dielectric continuum obtained from the minimal energy structure. Positions of individual binding energies in the minimum energy structure are indicated by red bars. The spectrum is broadened with the σ parameter equal to 0.34. The dashed lines correspond with the experimental data.109 The C-PCM model was used for modeling of water solvent. Optimal tuning resulted in an ω value of 0.303.

inferred from the spectra, the agreement with the experimental data is very good. The calculated peak positions are shifted by 0.2 eV compared to those of a previous study by Schroeder et al.,109 and the overall shape of the spectra has better agreement. The Koopman’s-like theorem seems to be fulfilled in the present case. Finally, arguably the most challenging system is the photoemission from liquid water. We tested the OT-IEDC approach on smaller finite-size water clusters (6−35 water molecules) similarly to our previous study.29 We included the effect of more distant solvent molecules via the polarized continuum. The simulated spectra are presented in Figure 5. The simulated spectra were fitted to four Gaussian functions; the peak widths and positions are collected in Table 3 and compared to the experimental values.29

ω = 0.14e−0.125x + 0.364

(13)

where x is the number of water molecules in the cluster. All of the values have approximately the same variance. The “saturation” is achieved for cluster sizes larger than 30. Note that the value of the range-separation parameter does not converge to the value found for an isolated water molecule embedded in a dielectric continuum. An ionized water cluster is clearly not an isolated water molecule surrounded by the other molecules acting as spectators as all the water units have approximately the same ionization energies.

■

DISCUSSION AND CONCLUSIONS We introduced and tested an alternative tuning scheme for the range-separated hybrid functionals. The protocol is based on enforcing the equality between the TDDFT-based HOMO ionization energy calculated as an excitation into a distant center (IEDC) and ΔKS energy. The functional is tuned to retrieve TDDFT excitation energy for a well-defined charge transfer state where the energetics can be independently calculated within the ΔKS framework. The OT-IEDC can be considered as an extension of the Koopmans’ theorem toward DFT where correlation effects are included by virtue of the TDDFT method. In this approach, all ionization energies are treated on equal footing. The accuracy of the ionization energy for the HOMO electron is determined by the quality of a given functional, and HOMO−n energies are estimated with reasonable accuracy. Using the OT-LC-ωPBE functional, we typically achieved an accuracy of ionization energies of several tenths of eV for neutral systems when compared to the CCSD or experimental reference. Further improvement of accuracy

Figure 5. Various approaches for photoemission spectra modeling of liquid water by the OT-IEDC method in comparison with experimental spectra.113 The photoemission spectrum was evaluated along the PI + GLE trajectory. Optimal tuning was performed for each geometry. Modeling of the spectra in the dielectric continuum was performed with the additional broadening scheme. 4978


Article


Table 3. Full Width at Half Maximum and Peak Positions of the Four Peaks in Photoemission Spectrum of Water Compared for Gas Phase and Dielectric Continuum Representing Bulk Water and Using the C-PCM Modela fwhm [eV]

ionization energy [eV]

peak no.

1st

2nd

3rd

4th

1st

2nd

3rd

4th

(H2O)6 (H2O)35 (H2O)6 in diel. cont. (H2O)35 in diel. cont. (H2O)35 exp.b liquid water exp.b

1.44 1.67 1.21 1.14 1.35 1.51

2.47 2.55 2.22 2.23 3.06 3.24

2.29 2.35 2.01 1.83 2.25 2.22

2.32 2.49 2.05 2.07

12.25 11.77 11.09 11.22 11.72 11.21

14.46 14.08 13.41 13.61 13.95 13.43

18.47 17.93 17.07 17.07 18.17 17.54

31.83 30.91 30.40 30.08

3.10

30.87

a

The photoemission spectrum was evaluated along the PI + GLE trajectory. Optimal tuning was performed for each geometry. Modeling of the spectra in the dielectric continuum was performed with the additional broadening scheme. bExperimental data were taken from ref 111. and fitted as discussed in ref 29.

to electronic excitations, thermochemical kinetics, and other applications. We tested the performance of the IEDC scheme also with hybrid functionals without optimal tuning using the BMK functional with a moderate share of the exact exchange and M06-HF functional with the full HF exchange. The IEDC approach with hybrid functional needs to be complemented by a relatively large ΔEcorr shift because hybrid functionals typically fail to describe charge transfer states. The shift cannot, however, be considered as a universal functional-dependent constant. Even the M06-HF functional fails to provide accurate numerical values of the ionization energies despite the full HF exchange. This finding is in accordance with other studies on systems with the important role of exact exchange.118 The benchmark set of molecules used in this work focuses on photoemission in the valence space. It can be, however, easily extended toward lower-lying electrons. Enforcing the energy condition in eq 7 for the HOMO electron might not be the best choice here, a better condition would be to enforce one of the core ionization energies calculated by the ΔKS method combined with, e.g., the maximum overlap method119 approach be equal to the TDDFT calculation. Such an approach could also appear necessary for more complex systems with multiple weakly interacting chromophoric units. It is by no means guaranteed however that the range-separation parameter optimized for one unit works equally well for the other units. The simplicity of the approach allows, for example, to efficiently model a whole spectrum of electron binding energies, e.g., to model photoemission spectra or redox potentials by simply saving molecular orbital energies during the ab initio molecular dynamics (MD) run.120,121 It also supports a simple spatial projection of calculated density of electronic states; in other words, it is possible to identify ionization within a certain molecular fragment or volume. This can be advantageous for modeling photoemission spectra of liquids by cluster models where the focus is only on the molecules inside the cluster or for modeling surface phenomena. The choice of sodium as the distant center might seem somewhat arbitrary. Indeed, we can use any other center with an unpaired electron having the HOMO ionization energy lower than the target system and HOMO−1 energy above the LUMO of the target system. Such a center can be, in principle, designed for each particular system, e.g., using pseudopotential technology. The present calculations can be conveniently performed with available quantum chemistry codes. Various ingredients are needed (range-separated functionals, state specific solvation, nonequilibrium solvation for ionization), but several computational codes can perform them. As the

Figure 6. Convergence of the ω parameter with water cluster size using the OT-IEDC method for every cluster size separately along with a comparison of the approach with the minimum energy geometry112 with the statistical value from the PIMD structure distributions.

would be connected to TDDFT itself, e.g., via including an orbital relaxation in advanced schemes.115,116 The OT-IEDC scheme was coupled with the C-PCM dielectric-based solvent model. The nonequilibrium solvation with a state-specific approach for the excited state needs to be consistently used for both the TDDFT and ΔKS calculations. In our implementation, we coupled the DFT/TDDFT calculations with the C-PCM model, and the suggested method emerged as a new efficient and accurate tool for modeling photoemission from liquid solutions besides existing approaches. Such quantitative approaches are highly required as investigating the structure and dynamics of liquid solutions has only recently been enabled by liquid microjet technology.1,2,117 We have to note here that the results for ions in the liquid phase show larger deviations from experimental values; this is however caused mostly by insufficiencies of polarized continuum models. Calculations of ionization energies for relatively small systems embedded in a dielectric continuum have practical as well as conceptual advantages. Although we could in principle perform DFT calculations of solvents in large solvent clusters, we would ultimately face the problem of size consistency.66 In conclusion, the OT-IEDC method represents a general tuning scheme that can be consistently applied for calculations of ionization energies in both the gas phase and solution. Moreover, the OT-IEDC approach yields an optimally tuned range-separated functional that can be straightforwardly applied 4979


Article


(14) Almbladh, C.-O.; von Barth, U. Exact results for the charge and spin densities, exchange-correlation potentials, and density-functional eigenvalues. Phys. Rev. B: Condens. Matter Mater. Phys. 1985, 31, 3231. (15) Perdew, J. P.; Parr, R. G.; Levy, M.; Balduz, J. L., Jr Densityfunctional theory for fractional particle number: derivative discontinuities of the energy. Phys. Rev. Lett. 1982, 49, 1691. (16) Perdew, J. P.; Levy, M. Comment on “Significance of the highest occupied Kohn-Sham eigenvalue. Phys. Rev. B: Condens. Matter Mater. Phys. 1997, 56, 16021. (17) Baer, R.; Livshits, E.; Salzner, U. Tuned range-separated hybrids in density functional theory. Annu. Rev. Phys. Chem. 2010, 61, 85. (18) Tozer, D. J.; De Proft, F. Computation of the hardness and the problem of negative electron affinities in density functional theory. J. Phys. Chem. A 2005, 109, 8923. (19) Leininger, T.; Stoll, H.; Werner, H.-J.; Savin, A. Combining long-range configuration interaction with short-range density functionals. Chem. Phys. Lett. 1997, 275, 151. (20) Vydrov, O. A.; Scuseria, G. E. Assessment of a long-range corrected hybrid functional. J. Chem. Phys. 2006, 125, 234109. (21) Jacquemin, D.; Perpete, E. A.; Scuseria, G. E.; Ciofini, I.; Adamo, C. TD-DFT performance for the visible absorption spectra of organic dyes: conventional versus long-range hybrids. J. Chem. Theory Comput. 2008, 4, 123. (22) Rohrdanz, M. A.; Martins, K. M.; Herbert, J. M. A long-rangecorrected density functional that performs well for both ground-state properties and time-dependent density functional theory excitation energies, including charge-transfer excited states. J. Chem. Phys. 2009, 130, 054112. (23) Iikura, H.; Tsuneda, T.; Yanai, T.; Hirao, K. A long-range correction scheme for generalized-gradient-approximation exchange functionals. J. Chem. Phys. 2001, 115, 3540. (24) Heyd, J.; Scuseria, G. E.; Ernzerhof, M. Hybrid functionals based on a screened Coulomb potential. J. Chem. Phys. 2003, 118, 8207. (25) Tawada, Y.; Tsuneda, T.; Yanagisawa, S.; Yanai, T.; Hirao, K. A long-range-corrected time-dependent density functional theory. J. Chem. Phys. 2004, 120, 8425. (26) Chai, J.-D.; Head-Gordon, M. Systematic optimization of longrange corrected hybrid density functionals. J. Chem. Phys. 2008, 128, 084106. (27) 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. (28) Salzner, U.; Baer, R. Koopmans’ springs to life. J. Chem. Phys. 2009, 131, 231101. (29) Hollas, D.; Muchová, E.; Slavíček, P. Modeling Liquid Photoemission Spectra: Path-Integral Molecular Dynamics Combined with Tuned Range-Separated Hybrid Functionals. J. Chem. Theory Comput. 2016, 12, 5009. (30) Egger, D. A.; Weissman, S.; Refaely-Abramson, S.; Sharifzadeh, S.; Dauth, M.; Baer, R.; Kü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. (31) Bokareva, O. S.; Grell, G.; Bokarev, S. I.; Kühn, O. Tuning Range-Separated Density Functional Theory for Photocatalytic Water Splitting Systems. J. Chem. Theory Comput. 2015, 11, 1700. (32) Autschbach, J.; Srebro, M. Delocalization Error and “Functional Tuning” in Kohn−Sham Calculations of Molecular Properties. Acc. Chem. Res. 2014, 47, 2592. (33) 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. (34) Jacquemin, D.; Moore, B.; Planchat, A. l.; Adamo, C.; Autschbach, J. Performance of an optimally tuned range-separated hybrid functional for 0−0 electronic excitation energies. J. Chem. Theory Comput. 2014, 10, 1677. (35) Laurent, A. D.; Jacquemin, D. TD-DFT benchmarks: A review. Int. J. Quantum Chem. 2013, 113, 2019.

orbital space included in the TDDFT calculations can be reduced, the computational cost of the OT-IEDC does not significantly exceed the cost of ground state calculations.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.7b00675. Additional data on the performance of IEDC with hybrid functionals (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. ORCID

Petr Slavíček: 0000-0002-5358-5538 Funding

This work was supported by the Czech Science Foundation GA Č R project No. 13-34168S. In addition, this project was supported by MŠMT project KONTAKT II LH15081. M.R. is a student supported by a specific university research project (MSMT No. 20-SVV/2017). Notes

The authors declare no competing financial interest.

■

REFERENCES

(1) Seidel, R.; Winter, B.; Bradforth, S. E. Valence Electronic Structure of Aqueous Solutions: Insights from Photoelectron Spectroscopy. Annu. Rev. Phys. Chem. 2016, 67, 283. (2) Winter, B.; Faubel, M. Photoemission from liquid aqueous solutions. Chem. Rev. 2006, 106, 1176. (3) Koopmans, T. Ü ber die Zuordnung von Wellenfunktionen und Eigenwerten zu den einzelnen Elektronen eines Atoms. Physica 1934, 1, 104. (4) Chong, D.; Gritsenko, O.; Baerends, E. Interpretation of the Kohn−Sham orbital energies as approximate vertical ionization potentials. J. Chem. Phys. 2002, 116, 1760. (5) Bellafont, N. P.; Illas, F.; Bagus, P. S. Validation of Koopmans’ theorem for density functional theory binding energies. Phys. Chem. Chem. Phys. 2015, 17, 4015. (6) Ortiz, J. V. Electron propagator theory: an approach to prediction and interpretation in quantum chemistry. Wiley Interdiscip Rev. Comput. Mol. Sci. 2013, 3, 123. (7) Govoni, M.; Galli, G. Large Scale GW Calculations. J. Chem. Theory Comput. 2015, 11, 2680. (8) Stanton, J. F.; Gauss, J. Analytic energy derivatives for ionized states described by the equation-of-motion coupled cluster method. J. Chem. Phys. 1994, 101, 8938. (9) Pieniazek, P.; Arnstein, S.; Bradforth, S.; Krylov, A.; Sherrill, C. Benchmark full configuration interaction and EOM-IP-CCSD results for prototypical charge transfer systems: noncovalent ionized dimers. J. Chem. Phys. 2007, 127, 164110. (10) Farrokhpour, H.; Ghandehari, M. Photoelectron spectra of some important biological molecules: symmetry-adapted-cluster configuration interaction study. J. Phys. Chem. B 2013, 117, 6027. (11) Roca-Sanjuán, D.; Rubio, M.; Merchán, M.; Serrano-Andrés, L. Ab initio determination of the ionization potentials of DNA and RNA nucleobases. J. Chem. Phys. 2006, 125, 084302. (12) Starcke, J. H.; Wormit, M.; Schirmer, J.; Dreuw, A. How much double excitation character do the lowest excited states of linear polyenes have? Chem. Phys. 2006, 329, 39. (13) Levy, M.; Perdew, J. P.; Sahni, V. Exact differential equation for the density and ionization energy of a many-particle system. Phys. Rev. A: At., Mol., Opt. Phys. 1984, 30, 2745. 4980


Article

Journal of Chemical Theory and Computation (36) Milanese, J. M.; Provorse, M. R.; Alameda, E., Jr; Isborn, C. M. Convergence of Computed Aqueous Absorption Spectra with Explicit Quantum Mechanical Solvent. J. Chem. Theory Comput. 2017, 13, 2159. (37) Cabral do Couto, P.; Hollas, D.; Slavíček, P. On the Performance of Optimally Tuned Range-Separated Hybrid Functionals for X-ray Absorption Modeling. J. Chem. Theory Comput. 2015, 11, 3234. (38) Verma, P.; Bartlett, R. J. Increasing the applicability of density functional theory. V. X-ray absorption spectra with ionization potential corrected exchange and correlation potentials. J. Chem. Phys. 2016, 145, 034108. (39) Song, J.-W.; Watson, M. A.; Nakata, A.; Hirao, K. Coreexcitation energy calculations with a long-range corrected hybrid exchange-correlation functional including a short-range Gaussian attenuation (LCgau-BOP). J. Chem. Phys. 2008, 129, 184113. (40) Besley, N. A.; Peach, M. J.; Tozer, D. J. Time-dependent density functional theory calculations of near-edge X-ray absorption fine structure with short-range corrected functionals. Phys. Chem. Chem. Phys. 2009, 11, 10350. (41) de Queiroz, T. B.; Kümmel, S. Charge-transfer excitations in low-gap systems under the influence of solvation and conformational disorder: Exploring range-separation tuning. J. Chem. Phys. 2014, 141, 084303. (42) de Queiroz, T. B.; Kümmel, S. Tuned range separated hybrid functionals for solvated low bandgap oligomers. J. Chem. Phys. 2015, 143, 034101. (43) Boruah, A.; Borpuzari, M. P.; Kawashima, Y.; Hirao, K.; Kar, R. Assessment of range-separated functionals in the presence of implicit solvent: Computation of oxidation energy, reduction energy, and orbital energy. J. Chem. Phys. 2017, 146, 164102. (44) Sun, H.; Ryno, S.; Zhong, C.; Ravva, M. K.; Sun, Z.; Körzdörfer, T.; Bredas, J.-L. Ionization Energies, Electron Affinities, and Polarization Energies of Organic Molecular Crystals: Quantitative Estimations from a Polarizable Continuum Model (PCM)-Tuned Range-Separated Density Functional Approach. J. Chem. Theory Comput. 2016, 12, 2906. (45) Zheng, S.; Geva, E.; Dunietz, B. D. Solvated charge transfer states of functionalized anthracene and tetracyanoethylene dimers: A computational study based on a range separated hybrid functional and charge constrained self-consistent field with switching Gaussian polarized continuum models. J. Chem. Theory Comput. 2013, 9, 1125. (46) Phillips, H.; Zheng, Z.; Geva, E.; Dunietz, B. D. Orbital gap predictions for rational design of organic photovoltaic materials. Org. Electron. 2014, 15, 1509. (47) Zhang, G.; Musgrave, C. B. Comparison of DFT methods for molecular orbital eigenvalue calculations. J. Phys. Chem. A 2007, 111, 1554. (48) Allen, M. J.; Tozer, D. J. Eigenvalues, integer discontinuities and NMR shielding constants in KohnSham theory. Mol. Phys. 2002, 100, 433. (49) Tozer, D. J. Relationship between long-range charge-transfer excitation energy error and integer discontinuity in Kohn−Sham theory. J. Chem. Phys. 2003, 119, 12697. (50) Zheng, Z.; Egger, D. A.; Brédas, J.-L.; Kronik, L.; Coropceanu, V. Effect of Solid-State Polarization on Charge-Transfer Excitations and Transport Levels at Organic Interfaces from a Screened RangeSeparated Hybrid Functional. J. Phys. Chem. Lett. 2017, 8, 3277. (51) Refaely-Abramson, S.; Sharifzadeh, S.; Jain, M.; Baer, R.; Neaton, J. B.; Kronik, L. Gap renormalization of molecular crystals from density-functional theory. Phys. Rev. B: Condens. Matter Mater. Phys. 2013, 88, 081204. (52) Ullrich, C. A. Time-dependent density-functional theory: concepts and applications; OUP: Oxford, U.K., 2011. (53) Pluhařová, E.; Slavíček, P.; Jungwirth, P. Modeling Photoionization of Aqueous DNA and Its Components. Acc. Chem. Res. 2015, 48, 1209.

(54) Svoboda, O.; Ončak, M.; Slavíč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. (55) Pluhařová, E.; Schroeder, C.; Seidel, R.; Bradforth, S. E.; Winter, B.; Faubel, M.; Slavíček, P.; Jungwirth, P. Unexpectedly small effect of the DNA environment on vertical ionization energies of aqueous nucleobases. J. Phys. Chem. Lett. 2013, 4, 3766. (56) Pluhařová, E.; Ončaḱ , M.; Seidel, R.; Schroeder, C.; Schroeder, W.; Winter, B.; Bradforth, S. E.; Jungwirth, P.; Slavíček, P. Transforming anion instability into stability: contrasting photoionization of three protonation forms of the phosphate ion upon moving into water. J. Phys. Chem. B 2012, 116, 13254. (57) Pluhařová, E.; Jungwirth, P.; Bradforth, S. E.; Slavíček, P. Ionization of Purine Tautomers in Nucleobases, Nucleosides, and Nucleotides: From the Gas Phase to the Aqueous Environment. J. Phys. Chem. B 2011, 115, 1294. (58) Holthausen, M.; Koch, W. A Chemist’s guide to density functional theory; Wiley-VCH: Weinheim, Germany, 2000. (59) Ipatov, A.; Cordova, F.; Doriol, L. J.; Casida, M. E. Excited-state spin-contamination in time-dependent density-functional theory for molecules with open-shell ground states. J. Mol. Struct.: THEOCHEM 2009, 914, 60. (60) Menon, A. S.; Radom, L. Consequences of spin contamination in unrestricted calculations on open-shell species: Effect of Hartree− Fock and Møller− Plesset contributions in hybrid and double-hybrid density functional theory approaches. J. Phys. Chem. A 2008, 112, 13225. (61) Stanton, J. F.; Gauss, J. A simple scheme for the direct calculation of ionization potentials with coupled-cluster theory that exploits established excitation energy methods. J. Chem. Phys. 1999, 111, 8785. (62) Ogi, Y.; Obara, Y.; Katayama, T.; Suzuki, Y.-I.; Liu, S.; Bartlett, N.-M.; Kurahashi, N.; Karashima, S.; Togashi, T.; Inubushi, Y. Ultraviolet photochemical reaction of [Fe (III)(C2O4) 3] 3− in aqueous solutions studied by femtosecond time-resolved X-ray absorption spectroscopy using an X-ray free electron laser. Struct. Dyn. 2015, 2, 034901. (63) Baugh, J.; Burkhardt, C.; Leventhal, J.; Bergeman, T. Precision Stark spectroscopy of sodium 2 P and 2 D states. Phys. Rev. A: At., Mol., Opt. Phys. 1998, 58, 1585. (64) Magyar, R.; Tretiak, S. Dependence of spurious charge-transfer excited states on orbital exchange in TDDFT: large molecules and clusters. J. Chem. Theory Comput. 2007, 3, 976. (65) Dreuw, A.; Head-Gordon, M. Failure of time-dependent density functional theory for long-range charge-transfer excited states: the zincbacteriochlorin− bacteriochlorin and bacteriochlorophyll− spheroidene complexes. J. Am. Chem. Soc. 2004, 126, 4007. (66) Whittleton, S. R.; Sosa Vazquez, X. A.; Isborn, C. M.; Johnson, E. R. Density-functional errors in ionization potential with increasing system size. J. Chem. Phys. 2015, 142, 184106. (67) You, Z.-Q.; Mewes, J.-M.; Dreuw, A.; Herbert, J. M. Comparison of the Marcus and Pekar partitions in the context of non-equilibrium, polarizable-continuum solvation models. J. Chem. Phys. 2015, 143, 204104. (68) Caricato, M.; Mennucci, B.; Tomasi, J.; Ingrosso, F.; Cammi, R.; Corni, S.; Scalmani, G. Formation and relaxation of excited states in solution: A new time dependent polarizable continuum model based on time dependent density functional theory. J. Chem. Phys. 2006, 124, 124520. (69) Guido, C. A.; Jacquemin, D.; Adamo, C.; Mennucci, B. Electronic Excitations in Solution: The Interplay between State Specific Approaches and a Time-Dependent Density Functional Theory Description. J. Chem. Theory Comput. 2015, 11, 5782. (70) Mewes, J.-M.; Herbert, J. M.; Dreuw, A. On the accuracy of the general, state-specific polarizable-continuum model for the description of correlated ground-and excited states in solution. Phys. Chem. Chem. Phys. 2017, 19, 1644. (71) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum mechanical continuum solvation models. Chem. Rev. 2005, 105, 2999. 4981


Article


(92) Ončaḱ , M.; Šištík, L.; Slavíček, P. Can theory quantitatively model stratospheric photolysis? Ab initio estimate of absolute absorption cross sections of ClOOCl. J. Chem. Phys. 2010, 133, 174303. (93) Della Sala, F.; Rousseau, R.; Görling, A.; Marx, D. Quantum and thermal fluctuation effects on the photoabsorption spectra of clusters. Phys. Rev. Lett. 2004, 92, 183401. (94) Ceriotti, M.; Manolopoulos, D. E.; Parrinello, M. Accelerating the convergence of path integral dynamics with a generalized Langevin equation. J. Chem. Phys. 2011, 134, 084104. (95) Ceriotti, M.; Bussi, G.; Parrinello, M. Colored-noise thermostats à la carte. J. Chem. Theory Comput. 2010, 6, 1170. (96) Klasinc, L.; Rušcǐ č, B.; Kajfež, F.; Šunjić, V. Photoelectron spectroscopy of the heterocycles imidazole and methylimidazoles. Int. J. Quantum Chem. 1978, 14, 367. (97) Rubešová, M.; Jurásková, V.; Slavíček, P. Efficient modeling of liquid phase photoemission spectra and reorganization energies: Difficult case of multiply charged anions. J. Comput. Chem. 2017, 38, 427. (98) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X. Advances in molecular quantum chemistry contained in the Q-Chem 4 program package. Mol. Phys. 2015, 113, 184. (99) Frisch, M.; Trucks, G.; Schlegel, H.; Scuseria, G.; Robb, M.; Cheeseman, J.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. Gaussian 09, revision D. 01; Gaussian, Inc., Wallingford, CT, 2009. (100) Slavíček, P.; Ončaḱ , M.; Hollas, D.; Svoboda, O. ABIN, version 1.0. https://github.com/PHOTOX/ABIN. (101) Linstrom, P. J.; Mallard, W. NIST Chemistry webbook; NIST standard reference database No. 69, 2001. (102) Tentscher, P. R.; Seidel, R.; Winter, B.; Guerard, J. J.; Arey, J. S. Exploring the aqueous vertical ionization of organic molecules by molecular simulation and liquid microjet photoelectron spectroscopy. J. Phys. Chem. B 2014, 119, 238. (103) Trofimov, A.; Schirmer, J.; Kobychev, V.; Potts, A.; Holland, D.; Karlsson, L. Photoelectron spectra of the nucleobases cytosine, thymine and adenine. J. Phys. B: At., Mol. Opt. Phys. 2006, 39, 305. (104) Debies, T.; Rabalais, J. Photoelectron spectra of substituted benzenes: II. Seven valence electron substituents. J. Electron Spectrosc. Relat. Phenom. 1972, 1, 355. (105) 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. (106) Page, R. H.; Larkin, R. J.; Shen, Y.; Lee, Y.-T. High-resolution photoionization spectrum of water molecules in a supersonic beam. J. Chem. Phys. 1988, 88, 2249. (107) Banna, M.; McQuaide, B.; Malutzki, R.; Schmidt, V. The photoelectron spectrum of water in the 30 to 140 eV photon energy range. J. Chem. Phys. 1986, 84, 4739. (108) Thürmer, S.; Seidel, R.; Winter, B.; Ončaḱ , M.; Slavíček, P. Flexible H2O2 in Water: Electronic Structure from Photoelectron Spectroscopy and Ab Initio Calculations. J. Phys. Chem. A 2011, 115, 6239. (109) Schroeder, C. A.; Pluhařová, E.; Seidel, R.; Schroeder, W. P.; Faubel, M.; Slavíček, P.; Winter, B.; Jungwirth, P.; Bradforth, S. E. Oxidation half-reaction of aqueous nucleosides and nucleotides via photoelectron spectroscopy augmented by ab initio calculations. J. Am. Chem. Soc. 2014, 137, 201. (110) Barth, S.; Ončaḱ , M.; Ulrich, V.; Mucke, M.; Lischke, T.; Slavíček, P.; Hergenhahn, U. Valence ionization of water clusters: From isolated molecules to bulk. J. Phys. Chem. A 2009, 113, 13519. (111) Winter, B.; Weber, R.; Widdra, W.; Dittmar, M.; Faubel, M.; Hertel, I. Full valence band photoemission from liquid water using EUV synchrotron radiation. J. Phys. Chem. A 2004, 108, 2625. (112) Wales, D. J.; Hodges, M. P. Global minima of water clusters (H2O) n, n≤ 21, described by an empirical potential. Chem. Phys. Lett. 1998, 286, 65.

(72) Ren, S.; Harms, J.; Caricato, M. An EOM-CCSD-PCM Benchmark for Electronic Excitation Energies of Solvated Molecules. J. Chem. Theory Comput. 2016, 13, 117. (73) Mennucci, B.; Cammi, R.; Tomasi, J. Excited states and solvatochromic shifts within a nonequilibrium solvation approach: A new formulation of the integral equation formalism method at the selfconsistent field, configuration interaction, and multiconfiguration selfconsistent field level. J. Chem. Phys. 1998, 109, 2798. (74) Cossi, M.; Barone, V. Separation between fast and slow polarizations in continuum solvation models. J. Phys. Chem. A 2000, 104, 10614. (75) Improta, R.; Barone, V.; Scalmani, G.; Frisch, M. J. A statespecific polarizable continuum model time dependent density functional theory method for excited state calculations in solution. J. Chem. Phys. 2006, 125, 054103. (76) Cammi, R.; Mennucci, B. Linear response theory for the polarizable continuum model. J. Chem. Phys. 1999, 110, 9877. (77) Tomasi, J.; Persico, M. Molecular interactions in solution: an overview of methods based on continuous distributions of the solvent. Chem. Rev. 1994, 94, 2027. (78) Cammi, R.; Corni, S.; Mennucci, B.; Tomasi, J. Electronic excitation energies of molecules in solution: State specific and linear response methods for nonequilibrium continuum solvation models. J. Chem. Phys. 2005, 122, 104513. (79) Cossi, M.; Barone, V. Time-dependent density functional theory for molecules in liquid solutions. J. Chem. Phys. 2001, 115, 4708. (80) Mewes, J.-M.; You, Z.-Q.; Wormit, M.; Kriesche, T.; Herbert, J. M.; Dreuw, A. Experimental Benchmark Data and Systematic Evaluation of Two a Posteriori, Polarizable-Continuum Corrections for Vertical Excitation Energies in Solution. J. Phys. Chem. A 2015, 119, 5446. (81) Jagoda-Cwiklik, B.; Slavíček, P.; Nolting, D.; Winter, B.; Jungwirth, P. Ionization of aqueous cations: Photoelectron spectroscopy and ab initio calculations of protonated imidazole. J. Phys. Chem. B 2008, 112, 7355. (82) Šištík, L.; Ončaḱ , M.; Slavíč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. (83) Duchemin, I.; Jacquemin, D.; Blase, X. Combining the GW formalism with the polarizable continuum model: A state-specific nonequilibrium approach. J. Chem. Phys. 2016, 144, 164106. (84) Muñoz-Losa, A.; Markovitsi, D.; Improta, R. A State-Specific PCM−DFT method to include dynamic solvent effects in the calculation of ionization energies: Application to DNA bases. Chem. Phys. Lett. 2015, 634, 20. (85) Jagoda-Cwiklik, B.; Slavíček, P.; Cwiklik, L.; Nolting, D.; Winter, B.; Jungwirth, P. Ionization of imidazole in the gas phase, microhydrated environments, and in aqueous solution. J. Phys. Chem. A 2008, 112, 3499. (86) Boese, A. D.; Martin, J. M. Development of density functionals for thermochemical kinetics. J. Chem. Phys. 2004, 121, 3405. (87) Zhao, Y.; Truhlar, D. G. Density functional for spectroscopy: no long-range self-interaction error, good performance for Rydberg and charge-transfer states, and better performance on average than B3LYP for ground states. J. Phys. Chem. A 2006, 110, 13126. (88) Rohrdanz, M. A.; Herbert, J. M. Simultaneous benchmarking of ground-and excited-state properties with long-range-corrected density functional theory. J. Chem. Phys. 2008, 129, 034107. (89) Cossi, M.; Rega, N.; Scalmani, G.; Barone, V. Energies, structures, and electronic properties of molecules in solution with the C-PCM solvation model. J. Comput. Chem. 2003, 24, 669. (90) Barone, V.; Cossi, M. Quantum calculation of molecular energies and energy gradients in solution by a conductor solvent model. J. Phys. Chem. A 1998, 102, 1995. (91) Truong, T. N.; Stefanovich, E. V. A new method for incorporating solvent effect into the classical, ab initio molecular orbital and density functional theory frameworks for arbitrary shape cavity. Chem. Phys. Lett. 1995, 240, 253. 4982


Article

Journal of Chemical Theory and Computation (113) Grell, G.; Bokarev, S. I.; Winter, B.; Seidel, R.; Aziz, E. F.; Aziz, S. G.; Kühn, O. Multi-reference approach to the calculation of photoelectron spectra including spin-orbit coupling. J. Chem. Phys. 2015, 143, 074104. (114) Sosa Vazquez, X. A.; Isborn, C. M. Size-dependent error of the density functional theory ionization potential in vacuum and solution. J. Chem. Phys. 2015, 143, 244105. (115) Pastore, M.; Assfeld, X.; Mosconi, E.; Monari, A.; Etienne, T. Unveiling the nature of post-linear response Z-vector method for timedependent density functional theory. J. Chem. Phys. 2017, 147, 024108. (116) Ronca, E.; Angeli, C.; Belpassi, L.; De Angelis, F.; Tarantelli, F.; Pastore, M. Density relaxation in time-dependent density functional theory: combining relaxed density natural orbitals and multireference perturbation theories for an improved description of excited states. J. Chem. Theory Comput. 2014, 10, 4014. (117) Seidel, R.; Thürmer, S.; Winter, B. Photoelectron spectroscopy meets aqueous solution: studies from a vacuum liquid microjet. J. Phys. Chem. Lett. 2011, 2, 633. (118) Cheng, X.; Zhang, Y.; Jónsson, E.; Jónsson, H.; Weber, P. M. Charge localization in a diamine cation provides a test of energy functionals and self-interaction correction. Nat. Commun. 2016, 7, 11013. (119) Gilbert, A. T.; Besley, N. A.; Gill, P. M. Self-consistent field calculations of excited states using the maximum overlap method (MOM). J. Phys. Chem. A 2008, 112, 13164. (120) Opalka, D.; Pham, T. A.; Sprik, M.; Galli, G. Electronic Energy Levels and Band Alignment for Aqueous Phenol and Phenolate from First Principles. J. Phys. Chem. B 2015, 119, 9651. (121) Liu, X.; Cheng, J.; Sprik, M. Aqueous Transition-Metal Cations as Impurities in a Wide Gap Oxide: The Cu2+/Cu+ and Ag2+/Ag+ Redox Couples Revisited. J. Phys. Chem. B 2014, 119, 1152.

4983


Optimal Tuning of Range-Separated Hybrids for Solvated Molecules

Recommend Documents