Solvatochromic Shifts in UV–Vis Absorption Spectra: The Challenging

Mar 11, 2016 - ... and of a High Performance Computing Center of the Matej Bel University in Banska Bystrica using the HPC infrastructure acquired in ...
1 downloads 0 Views 1MB Size
Article pubs.acs.org/JCTC

Solvatochromic Shifts in UV−Vis Absorption Spectra: The Challenging Case of 4‑Nitropyridine N‑Oxide

Šimon Budzák,† Adéle D. Laurent,† Christian Laurence,† Miroslav Medved’,‡ and Denis Jacquemin*,†,¶ †

CEISAM, UMR CNRS 6230, BP 92208, 2 Rue de la Houssinière, 44322 Nantes, Cedex 3, France Department of Chemistry, Faculty of Natural Sciences, Matej Bel University, Tajovského 40, SK-97400 Banská Bystrica, Slovak Republic ¶ Institut Universitaire de France, 103 bd Saint-Michel, F-75005 Paris, Cedex 05, France ‡

S Supporting Information *

ABSTRACT: 4-Nitropyridine N-oxide is a well-known molecular probe for which the experimental UV/vis absorption spectrum has been measured in a large number of solvents. Previous measurements and their analyses suggest a dominant role of the solvent hydrogen-bond donation (HBD) capability in the solvatochromic shifts measured for the absorption spectra. Herein, we analyze these solvatochromic effects using a series of complementary approaches, including empirical solvent parameters, high-level calculation of the excited-state dipole and polarizability, several flavors of the polarizable continuum model, as well as dynamics using an effective fragment potential (EFP) description of the solvent molecules. First, applying a recently proposed set of solvent parameters, we show the importance of dispersion interactions for non-HBD solvents. This statement confronts advanced coupled-cluster and multireference calculations of dipole moments and polarizabilities of both the ground and excited states in gas phase. We further address the pros and cons of implicit solvent models combined to time-dependent density functional theory (TD-DFT) in describing the solvents effects for all (HBD and non-HBD) media, the simplest linearresponse approach turning out to be the most adequate. Finally, we show that the explicit TD-DFT/EFP2 models work correctly for HBD molecules and allow for restoration of the main experimental trends. in a polymer film or frozen glass. Alternatively, UV−vis spectroscopic measurements of charge-transfer maxima of probe solutes in solution can be used to provide indirect information on the excess dipole moment and excess polarizability through an analysis of solvent−solute interactions by means of the Kamlet−Taft3 or similar equations, and such a procedure was used for 4-nitropyrine N-oxide4 studied here. Recently new sets of parameters describing solvents were developed by Catalán5 and our group.6 These new parametrizations enabled accurate estimations of dispersion and electrostatic solvent−solute effects to be reached in several cases.7−9 Besides experimental studies, one can also determine ES electrical properties by computational chemistry tools. However, in contrast with GS properties, several limitations appear because (i) ESs are often complex due to quasidegeneracy problems; (ii) the ordering of ESs can change during ES geometry optimization, which can be problematic, e.g., for the evaluation of so-called adiabatic excess properties; (iii) the accurate determination of electric excess properties generally requires larger basis sets compared to the GS properties; and (iv) the ES properties are often sensitive to contributions from double excitations. Consequently, the electric excess properties have mostly been computed only

1. INTRODUCTION The exploration of the electrical properties of electronically excited states (ES) with theoretical and experimental tools is a field of active research, which can be ascribed to several reasons. Let us emphasize that (i) the ES electrical properties largely influence solvatochromic effects that tune absorption and emission energies because they determine the strength of intermolecular interactions between an excited chromophore (a solute) and the surrounding solvent molecules; (ii) the knowledge of electric dipole moments and polarizability of the ES, often expressed as the difference between the excited and ground state values, i.e., the excess dipole moment (Δμ) and excess polarizability (Δα), is also vital for understanding several processes, e.g., the exciton localization/delocalization in conjugated systems.1 These conjugated molecules are indeed very intensively studied because their semiconducting and lightemitting properties make them attractive candidates for applications in both polymer or molecular electronic devices, such as light-emitting diodes and field effect transistors.2 Direct experimental measurements of ES properties are rather difficult to achieve except for tiny molecules. In fact, the main sources of experimental data for ES dipole moments and polarizabilities are the Stark effect and electroabsorption/ electrofluorescence measurements, where the shift of the absorption/fluorescence band is studied as a function of the applied electric field strength. These measurements can be carried out for molecules in gas phase, fluid solution, embedded © 2016 American Chemical Society

Received: February 10, 2016 Published: March 11, 2016 1919

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation for small systems, including water,10 acetone,11 and 2cyclopenten-1-one,12 which were examined thoroughly using advanced quantum chemical methods including electron correlation effects, such as the complete active space selfconsistent field method (CASSCF), second-order multiconfigurational perturbation theory (CAS-PT2), or equation of motion (EOM)/response coupled cluster (CC) methods. The main conclusion that can be drawn from these studies is that electric excess properties are strongly state-dependent and can reach high values. Descriptions of the electronic structures of ESs based on time-dependent density functional theory (TD-DFT) calculations are also increasingly popular due to their advantageous compromise between accuracy and computational cost. During the past decade, TD-DFT has indeed become the method of choice for evaluating the electronic excitation energies in large molecules for which advanced wave function-based methods cannot be applied in practice. It was shown that the combination of TD-DFT and standard numerical differentiations with respect to the external electric fieldalso known as the finite field techniqueyields reasonable estimates of the electric excess properties of uracil,13 benzonitrile,14 conjugated oligomers,15 and several other chromophores.16 TD-DFT also allows one to directly model the excited state in solutions when combined with implicit solvent models like the polarizable continuum model (PCM)17 or the solvent model with density (SMD).18 Recent years have brought significant developments in this area with methods going from a standard linear response (LR) method toward approaches explicitly accounting for electron density of excited states, e.g., corrected LR and self-consistent methods (see Computational Details). Nevertheless, the benchmarking of TD-DFT with respect to accurate wave function methods, including the control of the single excitation character of the studied ESs, remains necessary to reach reliable estimation of the excess properties. A more complex but computationally demanding class of models for describing solvent−solute interactions is constituted of the explicit approaches. In these theories, the solute is usually embedded inside a sufficiently large cluster of solvent molecules. Molecular dynamics (MD) or Monte Carlo (MC) techniques are then used to sample possible solvent−solute configurations that are used for excited state calculations. Because full ab initio MD/MC would require extremely large computational resources, molecular mechanics (MM) can be used for the MD/MC step while a QM/MM hybrid approach is applied for calculating excited state energies on the different snapshots in a reasonable amount of time.19−23 Here, we employed the effective fragment potential (EFP) method24,25 as an effective nonempirical force field. In this work, we use several approaches to investigate the solvatochromism of 4-nitropyridine N-oxide (Figure 1). This allows us to identify different difficulties connected with the description of ES electric properties by theoretical approaches. The solvent effects for this probe have mainly been studied in the ultraviolet part of the spectra, between 330 and 355 nm, and the related shifts have been used to probe hydrogen bonds.4 First, we reanalyze these experimental results with a modern set of solvent parameters to study the relative influences of the dispersion interactions and the change of dipole moment upon excitation. To this end, we calculate the GS and ES properties in gas phase by means of TD-DFT and more advanced techniques, such as CC and CAS-PT2. We next attempt to reproduce the experimentally observed absorption maxima shifts using several approaches based on both

Figure 1. The studied molecular probe, 4-nitropyridine N-oxide (top), and its frontier orbitals (HOMO, left; LUMO, right).

continuum and explicit models. Compared to previous works in this field, we offer a comprehensive assessment of several complementary theoretical and experimental viewpoints for a challenging molecular probe.

2. COMPUTATIONAL DETAILS As stated above, we have applied several approaches to model the solvatochromic shifts of the π → π* transition for 4nitropyridine N-oxide. 2.1. Empirical Approaches. Solvation capabilities of different solvents can be empirically described in terms of solvent dipolarity, polarizability, and hydrogen bond (HB) acidity and basicity. To study the influence of the solvent on the observed absorption spectra (or any physicochemical property A), one can transform the rather intuitive understanding of solvent−solute interactions into the following Kamlet−Taft− Abboud−Abraham equation3,26 A = A 0 + dδ + eπ * + aα + bβ (1) where A0 corresponds to the gas phase value of the property, and δ, π*, α, and β describe the solvent family-dependent polarizability, dipolarity/polarizability, HB acidity, and HB basicity, respectively. Experimental absorption spectra measurements in a large set of solvents are followed by fitting A of eq 1 to obtain the values of the corresponding regression coefficients (d, e, a, b) in linear solvation energy relationship (LSER) and, consequently, to quantify their importance. Catalán5 developed analogous LSER with differently defined solvent parameters. His SP, SdP, SA, and SB refer to the solvent’s polarizability, dipolarity, HB acidity, and HB basicity, respectively. To make the process of introducing a new solvent more efficient and accessible, some of us reparametrized6,27 the multiparameter expansion in eq 1. Parameter δ with discrete values was substituted by a new dispersion interaction (DI) parameter defined through the solvent’s refractive index. Next, new α1, β1, and ELS28 parameters as a measure of a solvent’s HB acidity, HB basicity, and solvent−solute electrostatic interactions, respectively, were introduced based on a combination of experimental measurements and theoretical calculations. The parameters of solvents used in this study are summarized in the Supporting Information (SI). 2.2. Gas Phase Calculations. The dipole moment and polarizability in the ground and excited states were calculated at the RHF, CIS, CIS(D), CC2, and CCSD levels of theory in the gas phase. In the case of coupled cluster methods, the response theory was applied as implemented in the Dalton program.29 In addition, CASSCF and CAS-PT2 calculations were performed with an active space of 14 electrons in 12 orbitals using Molcas 1920

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation software.30 We used the POL atomic basis set31 as a good compromise between accuracy and computational demands. The theoretical values were directly compared with the outcomes of the experimental analysis. 2.3. Implicit Models. Alternatively, one can use the calculated gas phase electric properties in equations based on either the Onsager model (eq 2) and/or the state-specific polarizability (SMSSP) approach (eq 3)32 to compare solvent shifts obtained with experimental results. In the Onsager model, the solvent shift (Δωsolv) is determined from the dipole moment change between the ground (μ0) and excited (μi) states 1 Δωsolv = − gd (μi − μ0 )2 2

nonequilibrium (fast process) regime. VEM calculations were performed with the VEMGauss software.38 2.4. Explicit Models. Explicit solvent molecules were considered in the framework of the EFP method, which is a nonempirical polarizable force field recognized for its accurate description of noncovalent interactions39 and for its efficiency in determining the excited energies of solvated systems.40,41 EFP2 molecular dynamics followed by the calculations of excitation energies with QM/EFP2 were performed to acquire deeper insights into the solvatochromism of the studied molecule. We used MD-EFP2 as implemented in Gamess42 and the QChem software43 was applied for QM/EFP2 calculations. In the EFP2 polarizable force field, the interactions between rigid molecules (fragments) are mediated through Coulombic, induction, exchange-repulsion, and dispersion terms. Components for individual fragments are calculated at the Hartree−Fock level in gas phase for individual solvent and solute molecules. The Coulomb interaction between the fragments is described using the distributed multipolar expansion (up to octopoles). The expansion centers are taken to be atom centers and bond midpoints. The induction part of the total interaction is represented by the interaction of the induced dipole on one fragment with the permanent dipole on another fragment and expressed in terms of the dipole polarizability. The molecular polarizability is expressed here as a tensor sum of localized molecular orbital (LMO) polarizabilities. Thus, the number of polarizable points is equal to the number of bonds and lone pairs in the molecule. This dipole-induced dipole interaction is solved iteratively to self-consistency and consequently partly includes many-body effects. The exchange repulsion term is obtained as an expansion in the intermolecular overlap between frozen LMOs on each fragment. Finally, by calculating dynamic polarizabilities over a range of frequencies for a fragment and integrating, one obtains C6 coefficients and interfragment dispersion contributions. Parameters for dichloromethane (DCM) and carbon tetrachloride were obtained according to a recommended setup,44 first performing geometry optimization at the HF/6311++G(3df,2p) level followed by an electrostatic calculation (the Stone localization) with the smaller 6-31+G(d) atomic basis set. Water, methanol, and pyridine parameters were taken from the library distributed within Gamess software. In the case of the probe molecule, we used the slightly larger aug-cc-pVTZ atomic basis set for the determination of all EFP2 parameters. This choice of a very flexible atomic basis set is justified because the gas phase dipole moment of the probe is close to zero but is composed of two counteracting groups with possibly different sensitivities to the solvent used. For short-range screening of polarization potential, we used the Tang−Toennis Gaussian formula with parameter a equal to 0.25 for the probe, whereas the default value of 0.6 was used for all solvents. This choice was based on test calculations. Namely, the SCS-MP2 geometry optimization of small clusters of the probe and two waters predicted N−O···H2O distances of ∼2.0 Å. The EFP2 optimization of this cluster with a default screening parameter (0.6) led to significantly shorter equilibrium distances (∼1.7 Å), and MD simulations for a larger cube of water molecules ended in overpolarization and nonphysical short distances between the solvent and the solute (∼1.1 Å) corresponding to the socalled polarization catastrophe. A smaller value of the screening parameter was able to dump this behavior while offering structures comparable with the SCS-MP2 results. Further data

(2)

where gd is a polarity function of the considered solvent. However, this model is limited because for nonpolar solvents the electrostatic interactions become less important, whereas dispersion interactions, not taken into account in eq 2, may play a major role. Marenich et al.32 developed a simple empirical correction for this interaction based on the difference between the ground (αGS) and excited state polarizability (αES) of the solute. An increase in polarizability induces solvent “dispersion” red shifts Δωsolv = −σd

n2 − 1 ES (α − αGS) n2 + 2

(3)

where σd is an empirical constant (1.850 kcal mol−1 Å−3), and n is solvent’s refractive index. The second family of methods used here relies on an implicit solvation picture based on the polarizable continuum model (PCM) going far beyond the simple Onsager model. For determining solvatochromic shifts, several PCM-TD-DFT formalisms can be applied: the linear response (LR) approach,33 which is the most widely used, the perturbative corrected LR (cLR) scheme,34 and the self-consistent statespecific models, i.e., IBSF (according to the names of the authors Improta, Barone, Scalmani, and Frisch)35 and the vertical excitation model (VEM).36 These formulations differ by the treatment of the change of the cavity polarization between the GS and the ES: (i) in the LR scheme, the transition densities are used to determine the change of the PCM charges upon excitation; (ii) in the cLR scheme, the one-particle TDDFT density matrix that accounts for the orbital relaxation contribution is used in a perturbation approach, and (iii) in both IBSF and VEM; an iterative self-consistent procedure is applied to evaluate the excited-state reaction field. In IBSF, the polarization effects are introduced into the calculation implicitly through the modification of the GS molecular orbitals and energies, whereas the VEM method treats such effects without modifying the GS wave function. Here, we evaluated the performance of the LR, cLR, and VEM methods. All DFT/TDDFT calculations were performed with the Gaussian software37 using the M06-2X functional, applying a so-called ultraf ine DFT integration grid (99,590 points), and the aug-cc-pVDZ atomic basis set. The possible basis set effects were assessed by using the larger aug-cc-pVTZ for several solvents (see the SI). No significant differences between solvent shifts calculated with the two basis sets were found. For each solvent, a geometry optimization in the selected solvent was first performed before the TD-DFT calculation of the excitation energies in the 1921

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation

and t-test value. The same holds for LSER based on the sole ELS parameter for this subgroup (for more details on the constructed LSER, see the SI). The relatively small value of the correlation coefficient for dispersion-based LSERs can be partly explained by the small range of variation (less than 0.2) of the DI parameter. We underline that solvents with high values of the DI parameter like CS2 and halogenated aromatic solvents (e.g., 1-iodonapthalene) as well as non-HBD solvents with small DI values like polyfluoroalkanes were not considered in the original experimental measurements. Next, we analyzed the full set of solvents. Only acetic acid was excluded from our analysis because the reference dye for establishing α1, betaine 30, is protonated in this solvent, and thus its α1 parameter is still unknown. Lagalante et al. also excluded acetic acid as well as water from their LSER constructions.4 In this group, the largest variance in solvent shifts is 3500 cm−1, corresponding to the difference between mesitylene and water. In HBD solvents, a hypsochromic displacement is observed. After constructing a correlation matrix for the full set of solvents, the highest correlation was observed for the α1 and SA parameters. Both our α16 and Catalán’s SA5 describe hydrogen bond donation capability of the solvent. LSER constructed purely using α1 has a regression coefficient of r = 0.928. The inclusion of the second parameter, DI, increases the correlation to r = 0.955, which covers 91.1% of the total variance in the different solvents. After normalizing LSER regression coefficients, we estimated the HB effect to account for 79% and dispersion to explain 21% of the total variance, keeping in mind that both effects act in opposite directions. Further addition of other terms brings no statistically significant improvement. To conclude the analysis of experimental results, we observed (i) a bathochromic effect for increasing the DI parameter in non-HBD solvents, suggesting an increase of polarizability upon excitation, (ii) strong hypsochromic effects for large α1 parameters, hinting that HBD solvents strongly interact with the probe, and (iii) no sensitivity to the ELS parameter, indicating a small/negligible change of the dipole moment upon excitation. 3.2. Ab Initio Results in Gas Phase. We have used several levels of theory to quantify the GS and ES dipoles and polarizabilities. Our results are collected in Table 1. The 4nitropyridine N-oxide molecule carries two electron-withdrawing groups on the pyridine ring in relative para position (see Figure 1). As a result, the GS dipole moment is small. The nitro group is, according to our calculations, a slightly stronger electron acceptor than the N-oxide; therefore, the total dipole moment vector points from the N-oxide moiety of the molecule toward the nitro group. The GS dipole moment has this direction for all methods, and its magnitude ranges from 0.277 D (RHF) to 1.119 D (CC2) (see Table 1). Our best estimate is 0.577 D calculated by the finite field approach on the basis of CCSD(T) energies. Because similar values are obtained with CAS-PT2, one can safely conclude that the GS dipole is indeed small. The first allowed vertically excited state belongs to the A1 symmetry representation, and it lies higher in energy than the forbidden B1, B2, and A2 states. The exact ordering of the states depends on the method used. The A1 transition is of π → π* character (see Figure 1). The electron is excited from an orbital partially localized on the oxygen atom of the N-oxide moiety and moves into an orbital localized on the opposite side of the molecule. Franck−Condon electron density redistribution thus

supporting this choice can be found in the SI (SAPT and EFP2 comparison). The EFP2 molecular dynamics simulations were performed using the canonical ensemble (NVT) at 300 K using the Nose− Hoover thermostat with a 1 fs time step. Each system was placed into a solvent box, and the periodic boundary conditions were imposed. The box size and corresponding number of solvent molecules was chosen to fit the experimental density of the solvent (1.00, 0.792, 1.33, 0.98, and 1.59 g cm−3 for water, methanol, DCM, pyridine (PYR), and tetrachloromethane, respectively). For exact box sizes and numbers of solvent molecules, see the SI. The system was first equilibrated (usually 10 and 20 ps for water and methanol, respectively) followed by a 10 ps production run during which 42 snapshots where taken for further excited state calculations at the Tamm−Damcoff (TDA)45 M06-2X/6-31G(d)/EFP2 level of theory. It is worth noting that the EFP MD simulations are usually undertaken with a rather small number of solvent molecules. In the present article, we decided to increase the box size (between 100 and 200 solvent molecules) to better reproduce solvent effects.

3. RESULTS AND DISCUSSION 3.1. Analysis of the Experimental Data with Empirical Approaches. Lagalante et al.4 measured the UV−vis spectra of 4-nitropyridine N-oxide in 48 solvents ranging from nonpolar n-heptane to polar and protic solvents, such as water and aliphatic alcohols. They used the empirical Kamlet−Taft parameters3 to establish an LSER (eq 1) for the charge transfer absorption in the tested solvents. The most significant parameter identified in this analysis was α, which measures the hydrogen bond donation ability of the solvent. This was shown to be sufficient to obtain a simple satisfying LSER. The correlation obtained with α only (r = 0.967) is not significantly improved when including the π* parameter (describing the solvent dipolarity and polarizability) nor the δ term specific to each solvent class (1 for aromatics, 0.5 for polychlorinated solvents, and 0 for all other solvents). As expected from the structure of the dye, the β parameter was statistically insignificant. Considering the hypsochromic shift for HBD solvents, the authors concluded that the dipole moment of the vertically excited state should be smaller than its ground state counterpart.4 As a first step in our analysis, we performed LSER for nonHBD solvents only (23 solvents, see the SI). We used the two sets of parameters compiled by Catalán5 and some of us.6 In this non-HBD subgroup, the solvent-induced shift attains 660 cm−1 between the probe in n-hexane (28920 cm−1) and in mesitylene (28260 cm−1), the former being a solvent with one of the smallest dispersion (DI) parameters,6 whereas the latter has the largest DI coefficient value in our solvent set. The redshift between these two solvents suggests a stronger dispersion interaction in the (Franck−Condon) excited state than in the ground state. Using only this parameter for construction of LSER leads to a moderate correlation coefficient of r = −0.784. The extrapolated gas phase transition energy obtained from this dependence is 30583 cm−1. Catalán’s solvent polarizability (SP) parameter5 has a similar physical meaning and yields a similar correlation coefficient of r = −0.771. The studied subgroup includes both aromatic (5) and nonaromatic (18) solvents but no significant difference between them could be observed. Our attempts to establish better LSER by adding the ELS parameter6 showed no statistical importance of ELS in terms of the standard error 1922

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation

that Δμ attains ca. 4 D, which is a relatively large change. This contradiction with the experimental assessment of the low importance of the ELS parameter cannot be explained using gas phase models only; thus, we addressed it with implicit and explicit solvent models (see the next sections). A common choice for calculations of excited states of larger molecules is TD-DFT. In the following section, we compare several approaches for the calculation of solvatochromic shifts using PCM-TD-DFT and the M06-2X functional. As can be seen in Table 1, the M06-2X GS and ES values of the dipole moments agree well with our best estimates, hinting that this functional is well-suited here. The DFT polarizability of the GS is also in excellent agreement with both CCSD(T) and CASPT2 results, whereas the respective value for ES is slightly smaller than the CCSD (unrelaxed) value. The difference comes from the most sensitive αzz, which is significantly smaller than the CC result. In the ES, the quality of the correlation energy treatment induces significant differences, when, for example, the change between the CC2 and CCSD value is 80 au3. Independently of the level of theory, an increase of polarizability can be observed upon excitation, which is consistent with the conclusion of the previous section. 3.3. Theoretical Results with Implicit Solvent Models. PCM-TD-DFT studies of excited states have become common in the field of computational chemistry, and they are also often used to support many investigations studies.47−49 Nevertheless, such works usually consider only a limited number of solvents (e.g., polar DCM and nonpolar cyclohexane), and in such cases, one often looks at large differences between two or three media. Here we consider a much more extended panel of solvents. Following our discussion of the experimental data, we first comment on the excitation energies in different non-HBD solvents (see Figure 2 and Table 2). Experimental analysis has shown that the dispersion energy represents for this subgroup

Table 1. Ground and Excited State Dipole Moments and Polarizabilities Calculated at Various Levels of Theory in Gas Phase method ground/excited state

ground state

excited state (A1)

a

Dipole Moment [D] 0.276 0.651 1.119 0.488 0.393 0.577 0.619 0.643 0.765 Polarizabilitya,c [au3] RHF/CIS 89.6 MP2/CIS(D) 97.2 CC2 (unrelaxed)b 105.7 CCSD (unrelaxed)b 97.0 (146.1) CCSD(T) (relaxed)b (140.6) CASSCF(14,12) (183) CAS-PT2(14,12) (150) M06-2X/TD-M06-2X 94.8 (142.6) RHF/CIS MP2/CIS(D) CC2 (unrelaxed)b CCSD (unrelaxed)b CCSD (relaxed)b CCSD(T) (relaxed)b CASSCF(14,12) CAS-PT2(14,12) M06-2X/TD-M06-2X

2.771 4.390 6.076 4.381

9.255 5.363 4.519 110.8 (167.7) (211.5) 166.7 (321.4) 140.1 (241.4) (165) (326) 117.1 (193.4)

a In previous work, Lazzaretti et al.46 calculated GS dipole and polarizability of 0.09 D and 87.6 au3, respectively. bRelaxed/unrelaxed corresponds to orbital relaxation during electric property calculations. c Values of αzz are reported in parentheses (z axis is parallel to the dipole moment vector).

leads to an increase of the dipole moment in the excited state with the same direction as in the GS. For the ES, the value of the dipole moment is more sensitive to the method used than in the GS. Our best estimates are 4.381 and 5.363 D at the CCSD and CAS-PT2 levels of theory, respectively. This means

Figure 2. Correlation between experimental absorption maximum and vertical excitation energies in different non-HBD solvents. Theoretical data were calculated at the M06-2X/aug-cc-pVDZ level of theory in combination with different flavors of the PCM model. 1923

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation

None of the theoretical approaches predicts mesitylene as being the most red-shifted. LR, SMSSP, and VEM-D favor pyridine as the refractive index in mesitylene is slightly smaller than the one of pyridine. Nevertheless experimental error bars overlap for pyridine and mesitylene. The highest excitation energy is measured in n-hexane and this result is accurately restored by both LR and VEM-D, whereas with SMSSP, diethyl ether yields the highest excitation energy. SMSSP is based only on the refractive index, whereas in both LR and VEM-D, the electrostatics also plays a role, which bathochromically shifts the n-hexane response. The spread of the values in this group is very small experimentally (660 cm−1 = 0.08 eV), and both LR and SMSSP predict an even smaller value (see Table 2). When including all solvents (both non-HBD and HBD) in the analysis, the correlation with experimental data gets worse. As can be seen in Table 2, no approach is able to follow the experimental trends with VEM/VEM-D even predicting the opposite trend compared to the experimental data. According to the experimental analysis, the formation of hydrogen bonds plays an important role in the blue shifting of excitation energy, which can explain the failure of implicit models that lack an explicit description of specific interactions, in this case, hydrogen bonds. 3.4. Theoretical Results Using Explicit Solvent Molecules. As shown above, implicit solvent models yield significant errors and do not reproduce several experimental trends. A natural next step in such cases would be to add several explicit solvent molecules to study their influence. There are arguments against such a static supermolecular approach because it does not catch the dynamics of the solvent, and more importantly, in our case, the probe contains two polar groups that need to be solvated on a balanced footing. This is why we decided to use molecular dynamics with each fragment described by its effective potential, covering repulsive exchange repulsion as well as electrostatic, polarization, and dispersion interactions. Our choice of media was made to cover important representatives of each subgroup of experimentally used solvents. From the HBD group, we have chosen water and methanol, the latter being sufficiently shifted from the water absorption maximum.52 DCM was included as a “bridge” between HBD and non-HBD groups of solvent (ELS = 0.6, α1 = 0.1), whereas PYR and CCl4 were chosen as representatives of non-HBD solvents. They present a series with decreasing polarity; moreover, one can possibly observe stacking interactions between PYR molecules and the solute. PYR is also a good representative of systems with dispersion-controlled interactions and is only slightly blue-shifted compared to mesitylene, which is a larger molecule and thus computationally more demanding. CCl4 represents the opposite side of the dispersion-dominated series of non-HBD, and contrary to nhexane, it is a very rigid molecule, which is ideal for the EFP methodology. After performing an equilibration step of 10−20 ps in EFP2-MD, we collected 42 snapshots for which the excitation energy was calculated at the TDA-M06-2X/631G(d)/EFP2 level (Table 3). The resulting histograms are depicted in Figure 3. In contrast with the gas phase where only one bright state is observed, intensity borrowing takes place in solution due to symmetry breaking such that the intensity of the third, fourth, and fifth excited state can become significant. If only one bright state is present, it exhibits a high oscillator strength of ∼0.5, which is similar to that in the gas phase. However, in 25−50% of the snapshots, we see several states with oscillator strengths exceeding 0.1. To cover such a

Table 2. Correlation Coefficients between Experimental and Calculated Excitation Energies Using Different Implicit and Empirical Models along with the Maximal Difference in the Transition Energies between the Solvents in the Group model

corr. coeff. r

max. diff. [cm−1]

Non-HBD Solvents Models Based on Continuum Models experiment 660 PCM/LR 0.798 310 PCM/cLR 0.434 1220 VEM 0.379 670 VEM-D 0.736 970 Models Based on Gas Phase Electric Properties Only Onsagera 0.552 390 SMSSPb 0.751 340 All Solvents Models Based on Continuum Models experiment 3500 PCM/LR 0.337 460 PCM/cLR 0.280 1240 VEM −0.674 1520 VEM-D −0.431 1920 a

Solute−solvent interactions based on eq 2. bSolute−solvent dispersion contribution based on eq 3 is assumed only.

the main contribution to solvent−solute interactions. In such a case, one would expect a rather low performance of models, such as PCM, mainly based on electrostatic considerations. The most technically advanced method used here is the VEM-D approach, which combines an iterative equilibration of the reaction field (VEM) with an empirically calculated dispersion correction obtained through eq 3. The Δα value was taken from the CCSD calculation in gas phase. Surprisingly, the VEM-D approach does not provide the most accurate results (see Figure 2 and Table 2). The correlation coefficient between the experimental and VEM-D excitation data is 0.736 only. This partial failure can mainly be ascribed to the electrostatic part of the total solvent shift, as can be seen from the poorer performance of the purely electrostatic VEM (r = 0.379); the SMSSP correction can only partly compensate the original trend. cLR, which can be considered as the first iteration of VEM, performs worse than VEM-D (due to the lack of dispersion) but slightly better than VEM. Globally, the LR model provides the best correlation with the experimental data (r = 0.798). There is an explanation for this apparently surprising result. Indeed, as pointed by Corni et al.,50 LR is successful due to a partial inclusion of dispersion interactions between the solvent and the solute: the dynamic term “−gdμ201” present in the LR expression comes from the response of the solvent to the electron density of the solute oscillating at the Bohr frequency. Such a term, having an origin in the dynamical solute−solvent interactions, might be classified as dispersionlike.51 The SMSSP model based on eq 3 is the second most successful approach with a correlation coefficient of r = 0.751, which is comparable to that of the LR method. This is consistent with the experimental analysis, indicating that the dispersion term dominates the interaction in non-HBD solvents. Both PCM as well as Onsager empirical methods based on pure electrostatics give significantly smaller correlation with the experimental data. The largest experimentally observed difference in this non-HBD group is between the excitation energies in n-hexane and mesitylene. 1924

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation

consisting in considering only the most intensive transition in all snapshots. One can observe in Figure 3 the correct ordering of strong HBD solvents with water being the most blue-shifted followed by methanol. The difference between the absorption maximum in water and in methanol is 2260 cm−1, which is in qualitative agreement with the experimental value (1300 cm−1). When using solvents that induce weaker specific interactions with the probe, the width of the spectra decreases, as expected. This change is already visible for the water−methanol pair. Because the hydrogen bond between the pyridine N−O of the probe and water generates a large solvent shift, the excitation energy is very sensitive to the parameters of this HB. During the MD simulations for the probe in water, the shortest hydrogen bond length (pyridine N−O···H2O) varies in the interval of 1.5−2.1 Å, explaining the larger width of the spectra. Although the variation of the hydrogen bond length for methanol is comparable, a slightly tighter spectra is obtained due to a weaker HB strength. Similar conclusions were obtained for the comparable case of p-nitroaniline in water, dioxane and cyclohexane.40 We underline that the probe as well as the solvent molecules were kept rigid during the MD simulation (normal EFP procedure), which implies that the

Table 3. Excitations Energies Calculated Using the Explicit TDA-M06-2X/6-31G(d)/MD/EFP2 Methodologya solvent

experiment

gas phase water methanol dichloromethane CCl4 pyridine

3.800b 3.938 3.776 3.581 3.569 3.509

method 1 4.858 5.281 4.995 4.824 4.800 4.837

(0.424) (0.138) (−0.034) (−0.058) (−0.021)

method 2 5.104 4.824 4.658 4.614 4.686

(−0.177) (−0.171) (−0.166) (−0.186) (−0.151)

a

See the text for an explanation of methods 1 and 2. Direct and indirect solvent effects are shown in parentheses. For method 1, the value is defined as the difference between the vertical excitation energy in solvent and in the gas phase. For method 2, we show only direct effects, i.e., resulting in the differences between methods 1 and 2. bThe gas phase value is extrapolated from LSER based on α1 and DI parameters (see SI, LSER 5).

scenario, for each snapshot, we first calculated the oscillator strength weighted excitation energy, and this value was used for further averaging over the snapshots. Nevertheless, averaging over the oscillator strength makes only a negligible difference (e.g., 0.002 eV for methanol) with the simpler approach

Figure 3. Simulated absorption spectra (histogram) of the probe molecule in (a) water, methanol, and pyridine and (b) dichloromethane and tetrachloromethane. 1925

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation

Figure 4. Comparison of solvent shifts calculated using explicit TDA-M06-2X/EFP2 and implicit solvation methods for selected solvents (PYR, CCl4, DCM, methanol, and water). SM stands for supermolecular approach.

other, and thus the final calculated (and experimentally observed) excitation energy is very close to the gas phase value despite the strong interaction between the solute and the solvent. Looking back to the experimental data, one can pinpoint a few pairs of solvents directly supporting the above analysis. For example, hexane and ethyl acetate are non-HBD solvents (α1 = 0) presenting the same value of the DI parameter, and the probe UV−vis spectrum in ethyl acetate is red-shifted (−180 cm−1) with respect to hexane as expected due to a stronger electrostatic interaction with the solvent. The same can be found for the pair of pentane-methyl acetate and cyclohexane-dimethylformamide (redshift of −250 cm−1). This shows that the electrostatic effect is experimentally present but unseen in the fits. During MD simulations with pyridine, we typically found two PYR molecules in stacking interactions with the probe, whereas two other molecules form very weak hydrogen bonds with the pyridine N−O oxygen atom. A typical length for this short contact is 2.2 Å, whereas the stacking occurs at distances of ∼3.6−4.0 Å. According to the experimental observations, the excitation energy in pyridine is redshifted compared to CCl4, whereas our calculations show the opposite trend. From a theoretical point of view, this is related to several different effects. First, method 1 already brings the excitation energy in pyridine higher than in CCl4, which is connected with the formation of a very weak hydrogen bond between pyridine and the solute. Similarly, dichloromethane is also able to form a weak HB with the solute, and therefore, according to method 1, both solvents yield excitation energies that are higher than that of CCl4, a completely aprotic solvent. The perturbative correction for polarization of solvent brings all the values down but does not change the ordering. A stronger redshift can be caused only by solvent−solute dispersion interactions, but those are not accounted for by TDA-DFT/EFP2. A rough estimate of the dispersion shift by SMSSP (eq 3) brings CCl4 and PYR even closer to each other, though this does not invert the ordering. For the case of pyridine, we also tested the supermolecular approach. Because specific interactions seem to be correctly described by EFP2, we kept the close contact pyridine molecules described by EFP2. In this way, we defined a QM region as the probe plus the pyridine molecules with the

results were not influenced by any possible (intramolecular) vibronic effects and, in this sense, are directly comparable with the continuum models results. The QM/EFP2 interaction contrasts with the full EFP2/ EFP2 only from Coulombic and polarization contributions with the former coming from the interaction of permanent multipoles of solvent molecules and the latter representing induced dipole moments of the solvent interacting with the solute wave function. The polarization term is solved until selfconsistence for each fragment and wave function. TDA-DFT/ EFP2 calculations include both direct and indirect contributions from the solvent. The indirect term comes from the orbital relaxation of the solute in the field generated by the permanent multipoles and induced dipoles of the solvent, whereas the direct term is the response of the polarizable environment to the change in solute electronic density upon excitation.53 To be consistent with the literature (e.g., ref 40), an approximation considering only indirect contribution is referred to as method 1, whereas the inclusion of both indirect and direct terms corresponds to method 2. For the studied molecule, we found that indirect and direct contributions have opposite effects. The indirect effect is significant for H-bonded water and methanol. In the case of water, this effect strongly blueshifts the excitation energy compared to the gas phase (+0.424 eV) and also remains large for methanol (+0.138 eV). For DCM, PYR, and CCl4, it acts in the opposite way and is smaller (max of −0.058 eV). The direct effect lowers the gas phase excitation energy by 0.15−0.19 eV. This effect is in line with the empirical trend (eq 2) and also with the larger calculated values of the excited state dipole moment, leading to a redshift upon solvation. Our values of dipole moments for the ground and excited states are thus not in contradiction to experimental findings but show two counteracting effects that cannot be distinguished from experimental fitting. To summarize, the indirect effect increases the excitation energy due to the formation of specific hydrogen bonds, whereas the direct effect leads to its decrease due to more favorable solvent−solute electrostatic interactions in the excited state. In practice, it means that the excitation energy in water, where the first effect is dominant, is blue-shifted with respect to the gas phase value. For methanol, both effects nearly cancel each 1926

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Journal of Chemical Theory and Computation



distance between the center of the probe and the center of the solvent molecules being less than 5 Å. In most snapshots, the QM region is constituted by the probe and the two stacked pyridines. The rest of the system was kept at the EFP2 level. Despite the increased size of the system, we did not find significant improvement with respect to the solvent described by EFP2 (see Figure 4). Method 1 with a larger QM region gives an excitation energy that is lower than that of pure EFP2 and now properly follows the experimental trend. Conversely, the result of method 2 for this partitioning of the system shows almost no change with respect to pure EFP2 for the solvent.

Article

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00149. Details of MD simulations, parameters used in empirical approaches and corresponding LSER equations, information on the basis set effects for PCM calculations, and comparisons between EFP2 and SAPT2 interaction energy components (PDF)



AUTHOR INFORMATION

Corresponding Author

4. CONCLUSIONS

*E-mail: [email protected].

In this work, we have explored several approaches to model solvatochromic effects for a probe. In the original experimental work, it was concluded from LSER analysis that 4-nitropyridine N-oxide was suitable for probing hydrogen bonds.4 It was also suggested that the ground state dipole moment was larger than the excited state dipole moment. However, although our CASPT2 and CCSD calculations reproduced the expected increase of polarizability (large dispersion effect), they systematically predicted a relatively large increase in dipole moment upon excitation in contrast with the LSER analysis. This inspired us to reanalyze the experimental data with other sets of empirical parameters for describing the solvent. These new analyses confirmed the importance of the ability of the solvent to act as a hydrogen bond donor but did not indicate any importance of the electrostatic interactions. This result is rather consistent with our analysis of different continuum models, hinting at significant dispersion effects, but is obviously not in line with the gas phase wave function results. To understand the origins of these divergent results, we simulated solvent−solute interactions using the MD/EFP2 technique and further used TDA-M06-2X/EFP2 for calculating excitation energies on the snapshots. We found good correlation with experimental data for the sequence water-methanol-dichloromethane-tetrachloromethane. In the case of strong HBDs (i.e., water and methanol), our analysis suggests that the excitation energy results from two counteracting effects. On the one hand, the formation of hydrogen bonds leads to stabilization of the HOMO, and this blue-shifts the spectra. On the other hand, the polarization response of the solvent to the density change of the solute following excitation systematically yields bathochromic shifts and is in line with the stronger electrostatic solvent−solute interactions in the excited state predicted by the dipole moment calculation. This explains the original contradictions. Taking a closer look back at the experimental data, one can find a few pairs of solvents that have no hydrogen bonding ability, the same dispersion effects, and still display bathochromic shifts with respect to the ELS parameter, which is consistent with theory. This work comes as a warning: small “fitted” electrostatic effects cannot be safely interpreted as negligible dipole moment change. This study has shown also that the available implicit/explicit solvent models nicely complement each other. The TDA-DFT/EFP2 approach works sufficiently well for systems with hydrogen bond interactions, whereas for solvents dominated by dispersion, LR (or SMSSP) satisfactorily reproduces the experimental trends.

Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS S.B. acknowledges the support of the European Research Council (ERC - Marches 278845 Grant) for his postdoctoral grant. D.J. acknowledges the ERC for financial support in the framework of a Starting Grant (Marches 278845). S.B., A.L., M.M., and D.J. are grateful to the Campus France and the Slovak Research and Development Agency for financial support of the bilateral collaboration in the framework of the PHC STEFANIK 2016 (35646SE) and project No. SK-FR-20150003, respectively. The authors thank the Région des Pays de la Loire for the LUMOMAT RFI project. This research used resources of the Centre de Calcul Intensif des Pays de Loire (CCIPL), a local Troy cluster, and of a High Performance Computing Center of the Matej Bel University in Banska Bystrica using the HPC infrastructure acquired in projects ITMS 26230120002 and 26210120002 (Slovak infrastructure for high performance computing) supported by the Research and Development Operational Programme funded by the ERDF.



REFERENCES

(1) She, C.; Easwaramoorthi, S.; Kim, P.; Hiroto, S.; Hisaki, I.; Shinokubo, H.; Osuka, A.; Kim, D.; Hupp, J. T. J. Phys. Chem. A 2010, 114, 3384−3390. (2) Friend, R. H.; Gymer, R. W.; Holmes, A. B.; Burroughes, J. H.; Marks, R. N.; Taliani, C.; Bradley, D. D. C.; Santos, D. A. D.; Bredas, J. L.; Logdlund, M.; Salaneck, W. R. Nature 1999, 397, 121−128. (3) Kamlet, M. J.; Abboud, J. L.; Taft, R. W. J. Am. Chem. Soc. 1977, 99, 6027−6038. (4) Lagalante, A. F.; Jacobson, R. J.; Bruno, T. J. J. Org. Chem. 1996, 61, 6404−6406. (5) Catalán, J. J. Phys. Chem. B 2009, 113, 5951−5960. (6) Laurence, C.; Legros, J.; Chantzis, A.; Planchat, A.; Jacquemin, D. J. Phys. Chem. B 2015, 119, 3174−3184. (7) Giordano, L.; Shvadchak, V. V.; Fauerbach, J. A.; Jares-Erijman, E. A.; Jovin, T. M. J. Phys. Chem. Lett. 2012, 3, 1011−1016. (8) Krystkowiak, E.; Koput, J.; Maciejewski, A. Phys. Chem. Chem. Phys. 2012, 14, 8842−8851. (9) El Seoud, O. A.; Ramadan, A. R.; Sato, B. M.; Pires, P. A. R. J. Phys. Chem. C 2010, 114, 10436−10443. (10) Palenikova, J.; Kraus, M.; Neogrády, P.; Kelloe, V.; Urban, M. Mol. Phys. 2008, 106, 2333−2344. (11) Pašteka, L. F.; Melicherčík, M.; Neogrády, P.; Urban, M. Mol. Phys. 2012, 110, 2219−2237. (12) Fišanova, J.; Č ernušaḱ , I.; Kelloe, V. J. Mol. Model. 2012, 18, 4751−4759. (13) Pluta, T.; Kolaski, M.; Medved’, M.; Budzák, Š. Chem. Phys. Lett. 2012, 546, 24−29. 1927

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation (14) Medved’, M.; Budzák, Š.; Pluta, T. Theor. Chem. Acc. 2015, 134, 78. (15) Grozema, F. C.; Telesca, R.; Snijders, J. G.; Siebbeles, L. D. A. J. Chem. Phys. 2003, 118, 9441−9446. (16) Budzák, Š.; Medved’, M.; Mennucci, B.; Jacquemin, D. J. Phys. Chem. A 2014, 118, 5652−5656. (17) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. Rev. 2005, 105, 2999−3094. (18) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2009, 113, 6378−6396. (19) Cerezo, J.; Avila Ferrer, F. J.; Prampolini, G.; Santoro, F. J. Chem. Theory Comput. 2015, 11, 5810−5825. (20) Pedone, A.; Prampolini, G.; Monti, S.; Barone, V. Phys. Chem. Chem. Phys. 2011, 13, 16689−16697. (21) Schwabe, T. J. Phys. Chem. B 2015, 119, 10693−10700. (22) Eilmes, A.; Kubisiak, P. J. Phys. Chem. B 2015, 119, 13185− 13197. (23) Etienne, T.; Assfeld, X.; Monari, A. Comput. Theor. Chem. 2014, 1040, 360−366. (24) Gordon, M. S.; Freitag, M. A.; Bandyopadhyay, P.; Jensen, J. H.; Kairys, V.; Stevens, W. J. J. Phys. Chem. A 2001, 105, 293−307. (25) Gordon, M. S.; Smith, Q. A.; Xu, P.; Slipchenko, L. V. Annu. Rev. Phys. Chem. 2013, 64, 553−578. (26) Taft, R. W.; Abboud, J.-L. M.; Kamlet, M. J.; Abraham, M. H. J. Solution Chem. 1985, 14, 153−186. (27) Cerón-Carrasco, J. P.; Jacquemin, D.; Laurence, C.; Planchat, A.; Reichardt, C.; Sraïdi, K. J. Phys. Chem. B 2014, 118, 4605−4614. (28) We use the acronym ELS for electrostatic interactions in this work instead of ES, as was introduced originally,6 to avoid confusion with the acronym of “excited state”. (29) Aidas, K.; Angeli, C.; Bak, K. L.; Bakken, V.; Bast, R.; Boman, L.; Christiansen, O.; Cimiraglia, R.; Coriani, S.; Dahle, P.; Dalskov, E. K.; Ekström, U.; Enevoldsen, T.; Eriksen, J. J.; Ettenhuber, P.; Fernández, B.; Ferrighi, L.; Fliegl, H.; Frediani, L.; Hald, K.; Halkier, A.; Hättig, C.; Heiberg, H.; Helgaker, T.; Hennum, A. C.; Hettema, H.; Hjertenæs, E.; Høst, S.; Høyvik, I.-M.; Iozzi, M. F.; Jansík, B.; Jensen, H. J. Aa.; Jonsson, D.; Jørgensen, P.; Kauczor, J.; Kirpekar, S.; Kjærgaard, T.; Klopper, W.; Knecht, S.; Kobayashi, R.; Koch, H.; Kongsted, J.; Krapp, A.; Kristensen, K.; Ligabue, A.; Lutnæs, O. B.; Melo, J. I.; Mikkelsen, K. V.; Myhre, R. H.; Neiss, C.; Nielsen, C. B.; Norman, P.; Olsen, J.; Olsen, J. M. H.; Osted, A.; Packer, M. J.; Pawlowski, F.; Pedersen, T. B.; Provasi, P. F.; Reine, S.; Rinkevicius, Z.; Ruden, T. A.; Ruud, K.; Rybkin, V. V.; Sałek, P.; Samson, C. C. M.; de Merás, A. S.; Saue, T.; Sauer, S. P. A.; Schimmelpfennig, B.; Sneskov, K.; Steindal, A. H.; Sylvester-Hvid, K. O.; Taylor, P. R.; Teale, A. M.; Tellgren, E. I.; Tew, D. P.; Thorvaldsen, A. J.; Thøgersen, L.; Vahtras, O.; Watson, M. A.; Wilson, D. J. D.; Ziolkowski, M.; Ågren, H. WIREs Comput. Mol. Sci. 2014, 4, 269−284. (30) Aquilante, F.; De Vico, L.; Ferré, N.; Ghigo, G.; Malmqvist, P.a.; Neogrady, P.; Pedersen, T. B.; Pitoňaḱ , M.; Reiher, M.; Roos, B. O.; Serrano-Andrés, L.; Urban, M.; Veryazov, V.; Lindh, R. J. Comput. Chem. 2010, 31, 224−247. (31) Sadlej, A. J. Collect. Czech. Chem. Commun. 1988, 53, 1995− 2016. (32) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Chem. Theory Comput. 2013, 9, 3649−3659. (33) Cammi, R.; Mennucci, B. J. Chem. Phys. 1999, 110, 9877−9886. (34) Caricato, M.; Mennucci, B.; Tomasi, J.; Ingrosso, F.; Cammi, R.; Corni, S.; Scalmani, G. J. Chem. Phys. 2006, 124, 124520. (35) Improta, R.; Scalmani, G.; Frisch, M. J.; Barone, V. J. Chem. Phys. 2007, 127, 074504. (36) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G.; Guido, C. A.; Mennucci, B.; Scalmani, G.; Frisch, M. J. Chem. Sci. 2011, 2, 2143− 2161. (37) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A.; Nakatsuji, H.; Caricato, M.; Li, X.; Hratchian, H. P.; Izmaylov, A. F.; Bloino, J.; Zheng, G.; Sonnenberg, J. L.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima,

T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M.; Heyd, J. J.; Brothers, E.; Kudin, K. N.; Staroverov, V. N.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Rega, N.; Millam, J. M.; Klene, M.; Knox, J. E.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Zakrzewski, V. G.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Dapprich, S.; Daniels, A. D.; Farkas, O.; Foresman, J. B.; Ortiz, J. V.; Cioslowski, J.; Fox, D. J. Gaussian09, revision E.01; Gaussian Inc.: Wallingford, CT, 2009. (38) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. VEMGAUSS, version 2013. http://comp.chem.umn.edu/vemgauss, University of Minnesota: Minneapolis, 2013. (39) Ghosh, D.; Kosenkov, D.; Vanovschi, V.; Williams, C. F.; Herbert, J. M.; Gordon, M. S.; Schmidt, M. W.; Slipchenko, L. V.; Krylov, A. I. J. Phys. Chem. A 2010, 114, 12739−12754. (40) DeFusco, A.; Minezawa, N.; Slipchenko, L. V.; Zahariev, F.; Gordon, M. S. J. Phys. Chem. Lett. 2011, 2, 2184−2192. (41) Slipchenko, L. V. J. Phys. Chem. A 2010, 114, 8824−8830. (42) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S.; Windus, T. L.; Dupuis, M.; Montgomery, J. A. J. Comput. Chem. 1993, 14, 1347−1363. (43) Shao, Y.; Gan, Z.; Epifanovsky, E.; Gilbert, A. T.; Wormit, M.; Kussmann, J.; Lange, A. W.; Behn, A.; Deng, J.; Feng, X.; Ghosh, D.; Goldey, M.; Horn, P. R.; Jacobson, L. D.; Kaliman, I.; Khaliullin, R. Z.; Kuś, T.; Landau, A.; Liu, J.; Proynov, E. I.; Rhee, Y. M.; Richard, R. M.; Rohrdanz, M. A.; Steele, R. P.; Sundstrom, E. J., III; Woodcock, H. L.; Zimmerman, P. M.; Zuev, D.; Albrecht, B.; Alguire, E.; Austin, B.; Beran, G. J. O.; Bernard, Y. A.; Berquist, E.; Brandhorst, K.; Bravaya, K. B.; Brown, S. T.; Casanova, D.; Chang, C.-M.; Chen, Y.; Chien, S. H.; Closser, K. D.; Crittenden, D. L.; Diedenhofen, M., Jr.; DiStasio, R. A.; Do, H.; Dutoi, A. D.; Edgar, R. G.; Fatehi, S.; Fusti-Molnar, L.; Ghysels, A.; Golubeva-Zadorozhnaya, A.; Gomes, J.; Hanson-Heine, M. W.; Harbach, P. H.; Hauser, A. W.; Hohenstein, E. G.; Holden, Z. C.; Jagau, T.-C.; Ji, H.; Kaduk, B.; Khistyaev, K.; Kim, J.; Kim, J.; King, R. A.; Klunzinger, P.; Kosenkov, D.; Kowalczyk, T.; Krauter, C. M.; Lao, K. U.; Laurent, A. D.; Lawler, K. V.; Levchenko, S. V.; Lin, C. Y.; Liu, F.; Livshits, E.; Lochan, R. C.; Luenser, A.; Manohar, P.; Manzer, S. F.; Mao, S.-P.; Mardirossian, N.; Marenich, A. V.; Maurer, S. A.; Mayhall, N. J.; Neuscamman, E.; Oana, C. M.; Olivares-Amaya, R.; O’Neill, D. P.; Parkhill, J. A.; Perrine, T. M.; Peverati, R.; Prociuk, A.; Rehn, D. R.; Rosta, E.; Russ, N. J.; Sharada, S. M.; Sharma, S.; Small, D. W.; Sodt, A.; Stein, T.; Stück, D.; Su, Y.-C.; Thom, A. J.; Tsuchimochi, T.; Vanovschi, V.; Vogt, L.; Vydrov, O.; Wang, T.; Watson, M. A.; Wenzel, J.; White, A.; Williams, C. F.; Yang, J.; Yeganeh, S.; Yost, S. R.; You, Z.-Q.; Zhang, I. Y.; Zhang, X.; Zhao, Y.; Brooks, B. R.; Chan, G. K.; Chipman, D. M.; Cramer, C. J., III; Goddard, W. A.; Gordon, M. S.; Hehre, W. J.; Klamt, A., III; Schaefer, H. F.; Schmidt, M. W.; Sherrill, C. D.; Truhlar, D. G.; Warshel, A.; Xu, X.; Aspuru-Guzik, A.; Baer, R.; Bell, A. T.; Besley, N. A.; Chai, J.-D.; Dreuw, A.; Dunietz, B. D.; Furlani, T. R.; Gwaltney, S. R.; Hsu, C.-P.; Jung, Y.; Kong, J.; Lambrecht, D. S.; Liang, W.; Ochsenfeld, C.; Rassolov, V. A.; Slipchenko, L. V.; Subotnik, J. E.; Van Voorhis, T.; Herbert, J. M.; Krylov, A. I.; Gill, P. M.; Head-Gordon, M. Mol. Phys. 2015, 113, 184−215. (44) Gordon, M. S.; Slipchenko, L.; Li, H.; Jensen, J. H. In Chapter 10 The Effective Fragment Potential: A General Method for Predicting Intermolecular Teractions; Spellmeyer, D., Wheeler, R., Eds.; Annual Reports in Computational Chemistry; Elsevier, 2007; Vol. 3, pp 177− 193. (45) Dreuw, A.; Head-Gordon, M. Chem. Rev. 2005, 105, 4009− 4037. (46) Lazzeretti, P.; Malagoli, M.; Turci, L.; Zanasi, R. J. Mol. Struct.: THEOCHEM 1993, 288, 255−259. (47) Zakrzewska, A.; Zalesny, R.; Kolehmainen, E.; Osmialowski, B.; Jedrzejewska, B.; Agren, H.; Pietrzak, M. Dyes Pigm. 2013, 99, 957− 965. 1928

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929

Article

Journal of Chemical Theory and Computation (48) Benelhadj, K.; Muzuzu, W.; Massue, J.; Retailleau, P.; CharafEddin, A.; Laurent, A. D.; Jacquemin, D.; Ulrich, G.; Ziessel, R. Chem. Eur. J. 2014, 20, 12843−12857. (49) Carlotti, B.; Benassi, E.; Cesaretti, A.; Fortuna, C. G.; Spalletti, A.; Barone, V.; Elisei, F. Phys. Chem. Chem. Phys. 2015, 17, 20981− 20989. (50) Corni, S.; Cammi, R.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 2005, 123, 134512. (51) McRae, E. G. J. Phys. Chem. 1957, 61, 562−572. (52) Aliphatic alcohols with longer carbon chains were not selected because, in EFP, we use rigid fragments and the dynamics of longer chains of the apolar part would be neglected (53) The direct term is evaluated as a perturbative correction, which is added to account to the solvent response to excited-state density. This is accomplished by recalculating the solvent-induced dipoles corresponding to this excited-state density. In some sense, this approximation is similar to the cLR method. A positive feature of such an approach is that the calculated excited states remain orthogonal to each other. We underline that such an approach does not account for the excited state density relaxation response to induced dipoles, which is expected to be a relatively small contribution in most cases. This small correction implies a demanding calculation54 that is methodologically hard to establish in the TD-DFT framework. (54) Arora, P.; Slipchenko, L. V.; Webb, S. P.; DeFusco, A.; Gordon, M. S. J. Phys. Chem. A 2010, 114, 6742−6750.

1929

DOI: 10.1021/acs.jctc.6b00149 J. Chem. Theory Comput. 2016, 12, 1919−1929