Combining Explicit Quantum Solvent with a Polarizable Continuum

Article pubs.acs.org/JPCB

Cite This: J. Phys. Chem. B 2017, 121, 10105-10117

Combining Explicit Quantum Solvent with a Polarizable Continuum Model Makenzie R. Provorse Long† and Christine M. Isborn*,‡ †

Department of Chemistry, University of Central Arkansas, Conway, Arkansas 72035, United States Chemistry and Chemical Biology, University of California Merced, Merced, California 95343, United States

‡

S Supporting Information *

ABSTRACT: A promising approach for accurately modeling both short-range and long-range solvation effects is to combine explicit quantum mechanical (QM) solvent with a classical polarizable continuum model (PCM), but the best PCM for these combined QM/classical calculations is relatively unexplored. We find that the choice of the solvation cavity is very important for obtaining physically correct results since unphysical double counting of solvation effects from both the QM solvent and the classical dielectric can occur with a poor choice of cavity. We investigate the dependence of electronic excitation energies on the definition of the PCM cavity and the self-consistent reaction field method, comparing results to large-scale explicit QM solvent calculations. For excitation energies, we identify the difference between the ground and excited state dipole moments as the key property determining the sensitivity to the PCM cavity. Using a linear response PCM approach combined with QM solvent, we show that excitation energies are best modeled by a solvent excluded surface or a scaled van der Waals surface. For the aqueous solutes studied here, we find that a scaled van der Waals surface defined by universal force field radii scaled by a factor of 1.5 gives reasonable excitation energies. When using an external iteration state-specific PCM approach, however, the excitation energies are most accurate with a larger PCM cavity, such as a solvent accessible surface.

■

INTRODUCTION The reliable prediction of condensed phase reactivity is important for modeling a variety of chemical processes, including catalysis and energy conversion. One approach to model both the short-range and long-range effects of the surrounding medium on a condensed phase reaction is to combine a quantum mechanical (QM) treatment of explicit solvent with a classical solvent model, such as molecular mechanical (MM) point charges or a polarizable continuum model (PCM).1−3 The QM/MM method typically uses electrostatic embedding to polarize the QM region with respect to fixed MM point charges. The PCM method encloses the QM charge density in a molecular cavity and uses apparent surface charges to compute the electronic response (i.e., polarization) of a structure-less dielectric continuum with respect to the QM charge density evaluated at points on the cavity surface. Treating explicit solvent with QM should capture solute−solvent interactions such as charge transfer and hydrogen bonding not accounted for via a classical solvent model. Including explicit QM solvent in the PCM molecular cavity gives a QM/PCM discrete-continuum solvation model that is promising for capturing short- and long-range solvation effects. The solvent dielectric constant represents the average ability of the solvent to stabilize charges at a given temperature, or in © 2017 American Chemical Society

the case of the PCM, the ability of the continuum solvent to stabilize the QM electron density. Because the solvent dielectric constant used in PCM calculations is an experimental value obtained from bulk solvent, it implicitly accounts for thermal averaging and other solvent effects not necessarily included in the QM electronic structure method, such as dispersion. On the other hand, explicit QM solvent gives the solvent response for a specific molecular configuration. Therefore, combining explicit QM solvent with a PCM represents a mismatch in solvent models at the QM/PCM interface. If bulk solvation effects are desired, combined QM/PCM solvent models require averaging over many QM solvent configurations. Despite this mismatch in configuration specific and bulk solvent effects, the combination of short-range specific effects modeled with explicit QM solvent and bulk solvent polarization modeled with a PCM can improve the agreement between experimental and predicted values of electronic excitation energies and solvatochromic shifts, especially for protic, polar solvents such as water4−8 but also for nonpolar solvents that interact with the solute via π-stacking.9 A discrete-continuum solvation model is also recommended to accurately predict solvation free energies Received: July 7, 2017 Revised: September 20, 2017 Published: October 9, 2017 10105

DOI: 10.1021/acs.jpcb.7b06693 J. Phys. Chem. B 2017, 121, 10105−10117

Article

The Journal of Physical Chemistry B and pKa values,10−15 interaction free energies,16 and rotational free energy barriers.17 However, the accuracy of some molecular properties, including solvatochromic shifts,4,7,18 proton affinities,19 solvation free energies and pKa values,13,20,21 NMR shifts,22 and metal ion−water distances,23 depends on the definition of the molecular cavity used to define the boundary between the discrete molecular system (i.e., QM region) and the continuum. Many PCM cavity definitions exist, but benchmarking of these cavities usually does not include explicit QM solvent and the reliability of PCM methods with explicit QM solvent is relatively unexplored. Given that condensed phase reactions are an important class of chemical transformations, a clear understanding of which PCM methods work best with explicit QM solvent would be extremely useful for the chemistry community. The definition of the PCM cavity is not unique. The simplest approach is to define the cavity as a sphere or ellipsoid, but this lacks the molecular detail necessary to accurately model local solute−solvent interactions.2,23,24 Many molecular-based cavities have been proposed and a schematic of common PCM molecular cavities is shown in Figure 1. The most common

Alternative definitions of the molecular cavity account for the physical size of solvent molecules. The solvent excluded surface (SES), also known as a Connolly surface, is defined by rolling a solvent sphere (black circle in Figure 1c) along the unscaled VDW surface.28 Using the GEPOL algorithm,29 an SES is approximated by adding spheres of different sizes (red and orange circles in Figure 1c) to fill in the space between the unscaled VDW surface and the solvent sphere.30 This surface is physically reasonable because it retains the molecular structure of the QM region, but excludes dielectric continuum from areas inaccessible to the solvent sphere. The solvent accessible surface (SAS) is defined by adding the radius of the solvent sphere to unscaled atomic radii to give an enlarged VDW-like surface (dashed black line in Figure 1d).31 In general, an SAS cavity has a larger volume and more surface area than an SES cavity. An SES is recommended for computing electronic properties without explicit QM solvent,2,30 but it is computationally more expensive than an SAS or a VDW surface because the additional spheres used to define an SES significantly increase the number of surface points at which the response of the dielectric continuum is computed. Some PCM methods use an isosurface to define the PCM cavity. For example, the isodensity PCM (IPCM) and selfconsistent IPCM methods use an isosurface based on the electron density to define the PCM cavity.32 The construction of a cavity based on an electronic isodensity surface has shown some success, including for the computation of quantum Monte Carlo/PCM excited state geometries using a PCM method that includes both surface and volume charges,33 but isodensity methods can often be plagued by convergence issues.34,35 The formulation of Fattebert and Gygi36,37 is the basis for other isodensity-based continuum solvation models.38,39 Continuum solvation models depend not only on the definition of the cavity, but also on the model of the electrostatic interaction between the QM region and the solvent dielectric. A PCM is one class of self-consistent reaction field (SCRF) methods that iteratively determine the QM-PCM electrostatic interaction in a self-consistent manner; therefore, the functional form of the effective Hamiltonian term used to model the QM-PCM electrostatic interaction will affect the final QM density. A widely used SCRF method is the integral equation formalism of PCM (IEF-PCM),40 which, in its implementation in Gaussian09,41 uses a continuous surface charge formalism42 to create a smooth and continuous reaction field. Conductor-like continuum solvation models use a screened conductor-like continuum to model bulk solvent; these methods include C-PCM, 43,44 COSMO, 45 and others.46−48 In addition to electrostatic interactions, nonelectrostatic contributions such as dispersion and cavitation energies may be included in continuum solvation models.4,49−52 The SMD SCRF method uses semiempirical functions to model the nonelectrostatic interactions in addition to the electrostatic interactions modeled by the IEF-PCM.49 Other SCRF methods include nonelectrostatic terms in the effective Hamiltonian to explicitly account for nonelectrostatic effects during the SCRF procedure, thus affecting the final QM density.50,52 Although the SMD nonelectrostatic terms do not alter the electron density from that obtained using the IEF-PCM, using a different set of atomic radii to define the VDW surface will change the electron density and the resulting excitation energies. Different sets of atomic radii are available to define

Figure 1. Schematic of molecular PCM cavities. Atomic spheres are blue, solvent spheres are black, and added spheres of different sizes are red and orange.

definition of a PCM molecular cavity is given by overlapping spheres centered on atomic nuclei (blue circles in Figure 1a), where the size of each sphere is defined by a van der Waals radius (i.e., VDW surface).2 The atomic radii used to define the VDW surface can be scaled by a parameter (α) to control the overlap of the nuclear spheres (see Figure 1b). This scaling approach is often used in continuum-only solvent models (i.e., no explicit QM solvent) to correct for local solute−solvent interactions not accurately modeled by a structure-less dielectric continuum,20,24−26 but is not typically used with mixed discrete-continuum solvent models. One exception is the work of Pliego and co-workers.27 They report that α ≥ 1.6 is necessary to compute physically reasonable solvation free energies of aqueous NH4+ and OH− with 32 or more explicit QM water molecules. Pliego and co-workers rationalized that smaller α values with large amounts of explicit solvent leads to “the formation of holes inside the cavity, which is occupied by the dielectric continuum.” These unphysical areas of dielectric continuum violate the assumptions of the apparent surface charge method invoked by PCM; increasing the α value removes these regions and recovers a physically reasonable solvation free energy. 10106


Article

The Journal of Physical Chemistry B the VDW surface, including the universal force field (UFF)53 radii and the intrinsic Coulomb radii of the SMD method,49 which we refer to as SMD radii. The SMD radii are parametrized for different solvents to reproduce experimental solvation free energies and are recommended for computing pKa values12 and vertical electronic excitation energies.4,7,54 Excited state calculations within the PCM framework can use a linear response (LR)55,56 or state-specific (SS)54,57−59 approach.60,61 Using time-dependent density functional theory (TDDFT),62,63 the LR PCM approach determines the electronic response of the dielectric continuum with respect to the ground-to-excited state transition density, whereas the SS PCM approach uses the difference between the ground and excited state electron densities to determine the response of the dielectric continuum. It has been shown that SS PCM methods are more accurate than LR PCM methods for modeling electronic transitions with significant charge-transfer character because they explicitly include the electronic response of the dielectric continuum to the excited state density, which is significantly different from the ground state density for chargetransfer transitions.64−66 In general, the accuracy and behavior of various SCRF methods for computing excited states when a large amount of explicit QM solvent is included in the calculation has not been thoroughly investigated. Advances in computational software and hardware have increased the size of the QM region that is computationally feasible, leading to studies examining the convergence of ground and excited state properties with respect to the size of the QM region.6,7,67−76 Our recent work compared the convergence of QM/MM and QM/PCM electronic excitation energies of aqueous solutes.77 We found that both QM/ classical methods require at least one solvation shell of explicit QM solvent to converge the excitation energy with respect to the size of the QM region. Since it is now possible to reach the size limits of QM computations that can yield converged molecular properties with the amount of QM solvent, these large-scale QM results can be compared to QM/PCM calculations to give insight into the best QM/PCM method for use with explicit QM solvent. For example, we found that the PCM cavity must be defined by an SES rather than a VDW surface (defined by UFF radii scaled by α = 1.1) to obtain physically reasonable QM/PCM excitation energies.77 With explicit QM solvent, we hypothesized that the VDW cavity artificially includes the response of the dielectric continuum for regions located between the solute and solvent molecules, whereas the SES “fills in” these regions with additional spheres, providing a more realistic model of solute−solvent polarization. We did not, however, test an SAS or scaled VDW surfaces, which we examine here. In this work, we investigate the dependence of QM/PCM excitation energies first on the definition of the PCM cavity and next on the SCRF method to identify the preferred QM/PCM approach for computing physically reasonable excitation energies of molecules surrounded with QM solvent. We directly compare QM/PCM and QM/MM excitation energies computed with a large number of QM solvent water molecules (1−2 solvation shells). For physically correct treatments of the PCM, we expect QM/MM and QM/PCM excitation energies to agree for large QM solvation regions, as they did in our previous work.77 Because water has a large dielectric constant, any excess polarization of the solute from double counting of solvation effects will be more substantial with aqueous solvation compared to a nonpolar solvent. We analyze TDDFT excitation

energies for three anionic solutes in water: trans-thiophenyl-pcoumarate (the chromophore of photoactive yellow protein, which we abbreviate as pCT−, 29 atoms) is characterized by a bright transition (oscillator strength f ∼ 1.2) with significant charge-transfer character from the phenolate oxygen to the π* system of the chromophore; phenolate (12 atoms) is a small, planar molecule that is characterized by a bright transition (f ∼ 0.2) with significant π → π* character; and the fluorescein anion (which we abbreviate as fluorescein−, 36 atoms) is a fluorophore with several low energy transitions; we choose to analyze a HOMO → LUMO + n transition of fluorescein− with moderate intensity (f ∼ 0.06). The TDDFT density difference plots for these transitions are shown in the Supporting Information, Figure S1.

■

COMPUTATIONAL DETAILS Classical Molecular Dynamics. For PCM with explicit QM solvent, solute−solvent configurations are typically optimized within the PCM cavity at zero temperature or extracted from a molecular dynamics simulation and inserted into a PCM cavity without further optimization. The latter approach accounts for thermal fluctuations within the solute− solvent system and provides more reliable solvatochromic shifts8,65,78,79 and pKa values;12 this is the approach used here. The computational details of the classical molecular dynamics simulation used to sample solute−solvent configurations may be found in ref 77. Each solute was solvated by a preequilibrated 32 Å sphere of TIP3P waters80 that included over 4000 water molecules. The Cartesian coordinates are provided in Table S1 for a single solute−solvent configuration of pCT− with 300 water molecules that is used throughout this work. QM/Classical Calculations. The Gaussian09 electronic structure program41 was used for all QM/classical calculations. All excitation energies were calculated using linear response TDDFT62,63 within the Tamm-Dancoff approximation.81 The long-range corrected hybrid functional ωB97X was used with the default amount of short-range exact exchange (about 16%) and the default range-separation parameter of 0.3 bohr−1.82 Although dispersion interactions between the QM solute and QM solvent may be important for computing accurate solutionphase properties, this is not the focus of the current work and no DFT corrections to the energy or density were included for dispersion interactions. The 6-31G(d) basis set was used for pCT− and the 6-31G basis set was used for phenolate and fluorescein−. Our tests suggest that the 6-31G and 6-31G(d) basis sets have similar agreement for the trends in the QM/MM and QM/PCM excitation energies computed with explicit QM solvent (Figure S2). The QM region was defined by the solute and the nearest water molecules. The remaining TIP3P water molecules80 were included as electrostatically embedded point charges in QM/MM calculations or removed in QM/PCM calculations. Unless otherwise noted, QM/PCM excitation energies were computed using nonequilibrium linear response PCM (LRQM/PCM),55,56 which is the default SCRF method for TDDFT calculations in Gaussian09. Nonequilibrium statespecific QM/PCM excitation energies (SS-QM/PCM) were computed by the external iteration method developed by Improta et al.58,59 using the Gaussian09 keywords TDA(root = N) and SCRF(ExternalIteration). This method requires two calculations. The solvent reaction field from a ground state calculation is saved using the Gaussian09 keyword SCRF(noneq = write) and subsequently read into an excited state 10107


Article

The Journal of Physical Chemistry B

moment, leading to increased polarization of the ground state electron density. This polarization would preferentially stabilize the ground state over the excited state, leading to an increase in the excitation energy, corresponding to a blue solvatochromic shift. Conversely, a positive Δμ0i (i.e., a smaller ground state dipole moment relative to the excited state dipole moment) would lead to greater polarization and stabilization of the excited state and a decreased excitation energy. Because the QM-PCM electrostatic interaction term in the energy expression is distance dependent, the extent of the stabilization provided by the PCM will be affected by the proximity of the PCM cavity surface points to the QM solute, which is determined by the PCM cavity definition. Thus, we hypothesize that the dependence of the QM/PCM excitation energy on the PCM cavity is related to the difference in the ground and excited state dipole moment (Δμ0i). We first compare the QM/classical excitation energies for different amounts of QM solvent for one configuration of three aqueous anionic solutes: pCT−, phenolate, and fluorescein−. The Δμ0i values for these structures are shown in Figure S3. The QM/MM and QM/PCM excitation energies for different amounts of QM solvent using the IEF-PCM method and various PCM cavity definitions are shown in Figure 2. We compare QM/PCM excitation energies using three common PCM cavity definitions: a VDW surface defined by UFF radii (scaled by α = 1.1, which is the default in Gaussian09),41 an SES, and an SAS. For the configurations of the three solutes examined in Figure 2, the negative Δμ0i values suggest that the excitation energy should be blue-shifted in going from vacuum to polar solvent. Experimentally, pCT− indeed has a large solvatochromic blue-shift of 0.44 eV.83 The present study does not aim to compute spectral shifts, which would require averaging over many solute and explicit solvent configurations, and TDDFT is likely not able to accurately model this shift for pCT− because of the challenging electronic structure of pCT− in vacuum.73,84−86 However, we report the values for a single configuration of each solute in order to examine the effects of explicit solvent combined with a PCM. The solvatochromic shifts going from vacuum to the PCM computed with the α = 1.1 VDW cavity and zero QM water molecules are 0.05 eV for pCT−, 0.23 eV for phenolate, and 0.36/0.65 eV for fluorescein− (two values are given since there were two states with nonzero oscillator strength in vacuum). A smaller solvatochromic shift is predicted for larger cavities, whereas a substantially larger solvatochromic shift is predicted when using MM point charges with zero QM water molecules. Upon inclusion of QM water molecules, all PCM excitation energies increase, leading to a larger solvatochromic shift than predicted with only a PCM. For all solutes, the PCM excitation energy with zero QM water molecules is underestimated compared to the converged result; the underestimation is substantial (0.3−0.7 eV) for phenolate and fluorescein−, but less so (∼0.1 eV) for pCT−. This underestimation of the PCM excitation energy is consistent with overestimation of the transition dipole moment with zero QM water molecules (Figure S4). If the nearest five water molecules are included in the QM region, the QM/PCM excitation energies are closer to the converged value. The underestimation of the excitation energy with zero QM water molecules shows that the ground state of these solutes is more stabilized by QM water molecules than by a PCM. All solutes with zero QM water molecules have a clear trend in the excitation energy with the size of the PCM cavity: smaller

calculation using the Gaussian09 keyword SCRF(noneq = read). The SS-QM/PCM excitation energy is then computed as a difference between the ground and excited state energies. Excited state dipole moments were also computed using this state-specific method. For the IEF-PCM40 and C-PCM43,44 SCRF methods, the VDW surface was defined using UFF53 radii scaled by α = 1.1 (the default in Gaussian09). The SMD method49 defines the VDW surface using unscaled (α = 1.0) SMD radii. Scaled VDW surfaces were defined using the Gaussian09 keyword Alpha = 1.2, 1.4, 1.5, 1.6, 1.8, 2.0, 2.5, or 3.0. The SES PCM cavity was defined using the Gaussian09 keyword Addsph with the default values of RMin = 0.200 and OFac = 0.890 and the default solvent radius of water (1.385 Å). The SAS PCM cavity was defined using the Gaussian09 keyword Surface = SAS and the default solvent radius of water (1.385 Å). All PCM calculations were computed using the default dielectric constants for water (ϵs = 78.355300, ϵ∞ = 1.777849). We tested the IPCM32 SCRF method for pCT− with up to 250 QM water molecules. The isodensity surface was computed using keywords Phi =10, 20, 40, 60, 80 and Theta = 5, 10, 20, 30, and 40 values per atom, respectively, to define the integration grid and an isovalue of 0.001. The ground and excited state calculations of the isodensity surface did not converge for QM regions with more than 15 QM water molecules. This problematic convergence is likely due to regions of low electron density between solvent molecules, which may cause large fluctuations in the cavity size and shape. An approximate IPCM excitation energy was computed using the isodensity surface defined by the first iteration of the excited state calculation, but gave large deviations from the converged QM/MM excitation energy (see Table S2); thus, further analysis of IPCM QM/PCM excitation energies is not presented here.

■

RESULTS AND DISCUSSION Comparison of PCM Cavities. In this section, we explore how the PCM cavity affects the excitation energy. The size of the cavity determines how close the surface representing the polarizable continuum is to the electron density, so the choice of cavity will affect the excitation energy regardless of whether or not there are explicit solvent molecules in the QM region. Our focus here is on how the accuracy of the excitation energy changes with choice of PCM cavity when including QM solvent. Accuracy is determined by comparison with the converged QM/MM result that includes 250−300 QM water molecules. Our previous work77 showed that the QM/MM excitation energy is not affected by the choice of nonpolarizable water model (i.e., TIP3P versus SPC/E); here, we use the TIP3P water model. To rationalize the accuracy of QM/PCM excitation energies computed using various PCM cavities, we use the change in the dipole moment upon electronic excitation. Similar analysis of the dipole moment and change in the dipole moment upon excitation was used to rationalize the prediction of aqueous solvatochromic shifts using a discrete/continuum solvation model.8 For most of the solute−solvent configurations studied here, the ground state dipole moment (μ0) is larger in magnitude than the excited state dipole moment (μi), yielding a negative value for the ground-excited state dipole moment difference (Δμ0i = μi − μ0). For these configuratons, we expect the dielectric medium to interact more strongly with the larger ground state dipole moment than the excited state dipole 10108


Article


computed with the largest QM region (300 QM water molecules) with ∼100 QM water molecules (Figure 2a). However, the α = 1.1 VDW cavity converges the QM/PCM excitation energy to a value 0.081 eV higher in energy than the converged QM/MM value. This overestimation of the excitation energy is due to the QM/PCM method including the unphysical response of the dielectric continuum at points between the QM solute and solvent molecules. We quantify the location of the points of the PCM dielectric by computing the minimum distance between the PCM cavity surface points and any QM solute nuclei in Figure 3. For pCT−, we see that smaller α values (α = 1.1, 1.2) give cavities with surface points that are closer to the solute than the explicit QM water molecules. These points overpolarize the QM electron density and preferentially stabilize the anionic ground state of pCT− relative to the excited state, which is consistent with a negative Δμ0i value (Figure S3).

Figure 2. QM/classical excitation energies of a single solute−solvent configuration of (a) pCT−, (b) phenolate, and (c) fluorescein− with increasing number of QM water molecules. Classical solvent is modeled using IEF-PCM (PCM) or fixed point charges (MM). For IEF-PCM calculations, the molecular cavity is defined by a VDW surface (using scaled UFF radii), an SES, or an SAS.

cavities (smaller α values) give higher excitation energies and larger cavities (larger α values) give smaller excitation energies. This result is consistent with the results of Improta and Barone that attribute a “smaller solvent effect” to the increased volume of cavities defined by larger α values.64 Similarly, we rationalize that smaller α values reduce the distance between PCM cavity surface points and the QM solute, which induces more polarization. This is reflected in larger magnitudes (by about 1−3 D) for the ground and excited state dipole moments with smaller α values. The increased excitation energies with smaller α values are consistent with negative Δμ0i values for zero QM water molecules (Figure S3), because the ground state is more stabilized by the enhanced QM-PCM interaction than the excited state. Focusing first on the QM/PCM excitation energies of pCT−, we see that as QM water molecules are added around the solute the QM/PCM excitation energies of pCT− computed with an α = 1.1 VDW, SES, and SAS cavities, as well as the QM/MM excitation energy, converge to within ±0.050 eV of the values

Figure 3. Distribution of the minimum distance between cavity surface points and any QM solute nucleus for various PCM cavity definitions for a single configuration of (a) pCT−, (b) phenolate, and (c) fluorescein− with 250 QM water molecules. The PCM cavity is defined by a VDW surface (using scaled UFF radii), an SES, or an SAS. The distribution of the minimum distance between QM solvent nuclei (Explicit) and any QM solute nucleus is shown using the right-hand yaxis. 10109


Article


water molecules to converge the solvation free energy.27 With 250 QM water molecules, the α = 2.5 and 3.0 VDW QM/PCM excitation energies deviate from the converged QM/MM value by −0.020 and −0.022 eV, respectively. Therefore, of all the PCM cavities examined here, the PCM cavities defined by an SES or an α = 1.5 VDW surface give the best agreement with converged QM/MM excitation energies. These cavities remove or reduce the unphysical response of the dielectric continuum evaluated at points within the explicit solvent shell, but do not push the cavity surface too far from the QM molecules as to create an artificial vacuum. Although both an SES and an α = 1.5 VDW surface give the correct physical result for pCT− with explicit QM solvent, an SES QM/PCM calculation is more computationally demanding. For example, the QM/PCM calculation for pCT− with 100 QM water molecules using an SES requires at least three times the number of CPU hours than a scaled VDW surface due to the computational cost of solving the PCM equations at the additional cavity surface points used to define an SES (31 914 points) compared to an α = 1.5 VDW surface (13 570 points). Thus, a scaled VDW cavity provides a computationally affordable alternative to an SES QM/PCM calculation. We next examine the convergence and accuracy of the QM/ classical excitation energies of phenolate (Figure 2b), which has a smaller Δμ0i value than pCT− (Figure S3). With 15 QM water molecules, the QM/PCM excitation energy of phenolate computed with an α = 1.1 VDW surface is converged to within ±0.050 eV of the value computed with 250 QM water molecules (which we take as the converged value). An SES and SAS require 50 QM water molecules to be converged, which again shows that larger PCM cavities require a larger QM region to yield converged excitation energies. With 250 QM water molecules, the SES and SAS QM/PCM excitation energies of phenolate agree with the QM/MM value with deviations less than ±0.005 eV, showing that the larger SES and SAS provide an accurate cavity for this system. The QM/PCM excitation energy computed with an α = 1.1 VDW surface overestimates the converged QM/MM value by 0.026 eV due to cavity surface points located within the QM solvent shell that are too close to QM solute and solvent molecules (Figure 3b). As the α value increases from 1.1 to 2.5, the QM/PCM excitation energy of phenolate with 250 QM water molecules systematically shifts down in energy (Figure 2b), similar to that shown in Figure 2a for pCT−. This shift in excitation energy is consistent with a negative Δμ0i value (Figure S3). Also, larger α values increase the distance between cavity surface points and the QM nuclei of phenolate and remove cavity surface points within the QM solvent shell (Figure 3b). Although scaled VDW QM/PCM excitation energies of phenolate depend on the α value, they deviate less from QM/ MM values than for pCT−. For example, at the point of convergence (∼50 QM water molecules), all PCM cavities tested here give QM/PCM excitation energies that deviate from the QM/MM excitation energy by less than ±0.011 eV compared to ±0.089 eV for pCT− (computed with ∼100 QM water molecules). This smaller sensitivity to the cavity is likely due to the smaller magnitude of Δμ0i for phenolate compared to pCT− (Figure S3). In fact, our results show that all PCM cavities tested here for phenolate with at least one solvation shell of explicit QM solvent give reasonable excitation energies, despite there being surface points located between QM solute and solvent nuclei for scaled VDW surfaces with small α values (Figure 3b).

An SES cavity provides the best agreement between QM/ PCM and QM/MM excitation energies when large amounts of QM solvent are included: with 300 QM water molecules, the deviation between the SES QM/PCM and QM/MM values is only +0.002 eV (Figure 2a). The good agreement is likely because an SES removes the unphysical surface points between the QM solute and the nearest QM solvent molecules and reduces the number of surface points located within the QM solvent shell (Figure 3a). An SAS cavity underestimates the converged QM/MM excitation energy by nearly 0.050 eV with 150 and 200 QM water molecules, which is reduced to about 0.010 eV with 250 and 300 QM water molecules (Figure 2a). This underestimation of the excitation energy using an SAS cavity is due to the cavity surface being too far from the QM solute and solvent molecules (Figure 3a), which reduces the electrostatic interaction between the QM density and dielectric medium and destabilizes the ground state relative to the excited state (consistent with a negative Δμ0i value). Overall, both SES and SAS PCM cavities give physically reasonable QM/PCM excitation energies with at least one solvation shell of explicit QM solvent, but an SES provides the best agreement with the converged QM/MM excitation energy. We next examine how scaling the VDW cavity, using UFF radii scaled by α = 1.1−3.0, changes the excitation energy. Increasing the α value increases the volume of the molecular cavity and, up to α ≈ 1.5, reduces the surface area by removing convex areas from the surface and “filling in” regions within the cavity between QM molecules (Figures S5 and S6). The increased volume of the cavity at larger α values increases the distance between cavity surface points and QM solute nuclei (Figure 3a), which decreases the ability of the PCM to stabilize the ground state. This change is reflected in a systematic shift to lower excitation energies for pCT− with increasing α value for most QM region sizes (Figure 2a). In fact, for a given QM region size, QM/PCM excitation energies computed with an increasing α value approximately approach the value computed without a PCM (Figure S7). Note that with 40 QM water molecules, the excited state of pCT− has a larger dipole moment than the ground state, giving a positive Δμ0i value (Figure S3); as the α value increases, the weaker QM-PCM interaction destabilizes the excited state relative to the ground state, increasing the excitation energy with larger α values (Figure 2a). The QM/PCM excitation energies computed with an α = 1.5 VDW surface agree with those computed with an SES for QM regions with at least 75 QM water molecules (Figure 2a). The good agreement is due to the location of the surface points being similar for an α = 1.5 VDW and an SES (Figure 3a). The α = 2.0 VDW QM/PCM excitation energies agree with those computed with an SAS for all QM region sizes (Figure 2a) because an α = 2.0 VDW surface pushes the bulk of the cavity surface points past the largest QM solute−solvent distance, consistent with an SAS (Figure 3a). Larger α values (α = 2.5 and 3.0) increase the molecular cavity volume and surface area above that of an SAS cavity (Figure S5), which increases the distance between cavity surface points and QM atoms (Figure 3a). With these large α values, an artificial vacuum is formed between the QM solvent molecules furthest from the solute and most of the cavity surface points. As a result, the α = 2.5 and 3.0 VDW surfaces require larger QM regions (∼250 QM water molecules) to converge the QM/PCM excitation energy to within ±0.050 eV of the values computed with 300 QM water molecules (Figure 2a). This result is consistent with de Lima et al. where a larger α value required more explicit QM 10110


Article

The Journal of Physical Chemistry B We next analyze the trends in the QM/classical excitation energies of fluorescein−, which has a small, negative Δμ0i value for most configurations (Figure S3). All QM/classical excitation energies of fluorescein− (Figure 2c) require ∼75 QM water molecules to be converged within ±0.050 eV of the values computed with the largest QM region (250 QM water molecules). A similar shift in QM/PCM excitation energy with different PCM cavities is observed for fluorescein− as seen with pCT− and phenolate. For example, with 100 QM water molecules, the QM/PCM excitation energy of fluorescein− shifts down in energy with increasing α value, which is consistent with larger distances between the QM nuclei and the PCM cavity surface points (Figure 3c) and a negative Δμ0i value (Figure S3). The magnitude of the shift in excitation energy is on par with that of phenolate. With 250 QM water molecules, all PCM cavities tested here give QM/PCM excitation energies of fluorescein− that deviate from the QM/ MM value by less than +0.015 eV, which is consistent with a small magnitude of Δμ0i (Figure S3). Overall, the dependence of the QM/PCM excitation energies of fluorescein− on the PCM cavity is similar to phenolate and less significant than pCT−. To test if our results vary with individual solute configurations, we compare results with various PCM cavity definitions for five solute−solvent configurations of pCT−, phenolate, and fluorescein− with one layer of QM solvent. Accuracy is determined by the deviation of the QM/PCM excitation energies from the converged QM/MM values (Figure 4). The Δμ0i values for each configuration are given in Table 1, which range from +1.1 to −5.8 D. For the results shown in Table 1 and Figure 4, the converged QM/MM excitation energies were computed using 250 QM water molecules and the QM/PCM excitation energies were computed using 50 QM water molecules for phenolate and 100 QM water molecules for pCT− and fluorescein−. The QM/ PCM excitation energies of Configuration 3 of pCT− are not converged with 100 QM water molecules, resulting in significant underestimation of the QM/PCM excitation energy relative to the converged QM/MM value for an SES, an SAS, and most α values of scaled VDW surfaces (Figure S8). With 250 QM water molecules, the QM/PCM excitation energies of Configuration 3 of pCT− are converged and have similar dependence on the PCM cavity as the other configurations of pCT− (Figure S9). The results for Configuration 3 of pCT− with 100 QM water molecules are not included in our analysis below. In Figure 4, panels a−c show the deviation of QM/PCM excitation energies computed with VDW cavities defined by UFF radii scaled by α = 1.1−3.0. For pCT− (Figure 4a), three of the four configurations give QM/PCM excitation energies that deviate from the converged QM/MM value by more than ±0.050 eV with a maximum deviation of ∼0.130 eV. The α value that minimizes the deviation is α = 1.5; however, α values of 1.4−1.8 give deviations less than ±0.050 eV. All solute− solvent configurations of pCT− have a negative difference dipole moment and show an inverse relationship between the QM/PCM excitation energy and α value, which is consistent with the trends observed in Figure 2a. Of the pCT − configurations, Configuration 2 has the smallest difference dipole moment and the least variation in excitation energy with different PCM cavities. Unlike for pCT−, the deviations of the QM/PCM excitation energies for scaled VDW cavities of all five configurations of

Figure 4. Deviation of QM/PCM excitation energies (ωPCM N ) from QM/MM excitation energies computed with 250 QM water molecules − (ωMM 250 ) for five solute−solvent configurations of (a) pCT (N = 100 QM water molecules), (b) phenolate (N = 50 QM water molecules), and (c) fluorescein− (N = 100 QM water molecules). For panels a−c, the PCM cavity is defined by a scaled VDW surface with increasing atomic radii scale factor (α). For panels d−f, the PCM cavity is defined by an SES or an SAS. See Figure S8 for Configuration 3 of pCT− results. Configuration 1 of each solute is analyzed in Figures 2 and 3.

Table 1. Difference between Ground (μ0) and Excited State (μi) Dipole Moments (Δμ0i = μi − μ0) for Three Solutes and the Surrounding QM Solvent in Debye (D)a Configuration Configuration Configuration Configuration Configuration

1 2 3 4 5

pCT− b

phenolatec

fluorescein− b

−4.14 −2.35 −5.17 −5.87 −5.82

+1.10 −1.43 −0.62 −0.54 −1.09

−0.59 −2.43 −0.06 −0.26 +0.37

a

Computed for a QM solute and N number of QM water molecules using nonequilibrium state-specific IEF-PCM and an SES cavity defined by UFF radii. bN = 100. cN = 50.

phenolate (Figure 4b) and fluorescein− (Figure 4c) are near or less than ±0.050 eV for all α values tested. Some configurations (e.g., Configuration 3 of fluorescein−) show negligible changes in QM/PCM excitation energy for different α values. For two configurations, Configuration 1 of phenolate and Configuration 5 of fluorescein−, the excited state dipole moment is larger than the ground state dipole moment, leading to a positive Δμ0i 10111


Article


calculation with large amounts of explicit QM solvent, but tests should be performed to determine the appropriate α value. Comparison of SCRF Methods. In the previous section, we focused solely on the effect of the PCM cavity on the QM/ PCM excitation energies. The formulation of the electrostatic interaction between the QM region and the solvent dielectric, i.e. the reaction field that is produced, will also affect the QM/ PCM excitation energies, which we focus on in this section. Comparison of IEF-PCM, SMD, and C-PCM Methods. The results presented in Figures 2−4 were computed using IEFPCM with cavities defined by UFF radii (the default SCRF method and atomic radii in Gaussian09). The SMD method includes a nonelectrostatic term in the total energy expression to account for cavitation, dispersion, and local solvent structure, in addition to the electrostatic term computed by IEF-PCM.49 The nonelectrostatic term is added to the total ground state energy and does not contribute to the excited state calculation. Therefore, excitation energies computed by the SMD method are equivalent to IEF-PCM excitation energies if the cavity is defined by the same surface (e.g., VDW, SES, SAS) and the same set of atomic radii (e.g., UFF, SMD). For example, the SMD and IEF-PCM QM/PCM excitation energies of pCT−with increasing amounts of QM solvent agree when the same cavity is used for both methods (Table S3). On the other hand, the C-PCM method is a conductor-like PCM method that computes the QM-PCM electrostatic interaction differently than IEF-PCM. In our analysis below, we compare first the effects of replacing UFF radii with SMD radii using the IEFPCM method and second the effects of using C-PCM versus IEF-PCM, where UFF radii are used for both SCRF methods. For each case, the computed QM/PCM excitation energies are shown in Figure 5 for the same solute−solvent configuration of pCT− as in Figure 2a.

value (Table 1). Only these two configurations show an increase in QM/PCM excitation energy with increasing α value, whereas the remaining configurations show a decrease in QM/ PCM excitation energy. These trends agree with our hypothesis regarding overpolarization by cavities with smaller α values and underpolarization by cavities with larger α values. Therefore, the dependence of the QM/PCM excitation energy on the PCM cavity may be predicted by the sign and magnitude of Δμ0i. In addition, solute−solvent configurations with larger Δμ0i values are more sensitive to the PCM cavity definition and care must be taken to ensure an appropriate cavity is used. For phenolate and fluorescein−, the best agreement with converged QM/MM values for all five configurations is given by α = 1.4 VDW cavities, whereas α = 1.5 minimizes the deviation for pCT−. These α values give QM/PCM excitation energies similar to those computed by an SES cavity (panels d− f in Figure 4), which may be attributed to the fact that α = 1.4 or 1.5 minimizes the surface area of the scaled VDW cavity of each solute surrounded by a large amount of explicit QM water molecules (see Figure S5). The VDW cavities scaled by α = 1.8 or 2.0 give similar QM/PCM excitation energies as SAS cavities. These results for five solute−solvent configurations of pCT−, phenolate, and fluorescein− are consistent with those presented in Figure 2. As a check on our interpretation for these anionic solutes in a mixed QM/PCM solvent environment, we also tested two solute−solvent configurations of the neutral phenol solute for how the radii scale factor for the cavity affected the excitation energy. The results given in Figure S10 are consistent with those presented for the anionic solutes in Figure 4. The changes in excitation energy with cavity size correlate with the sign and magnitude of the difference dipole moment, and for the two configurations of phenol solvated by 50 QM water molecules, an SES cavity and scaled VDW cavities with α = 1.1−1.4 give results in good agreement with the QM/MM excitation energy with 250 QM water molecules. Taken together, our results show that the addition of explicit QM solvent in a QM/PCM approach leads to an unphysical PCM cavity when a VDW surface scaled by a small α value is used. PCM cavities defined by an SES, SAS, or scaled VDW surface with a large α value (α ∼ 1.5) remove PCM cavity surface points near the QM solute nuclei and increase the distance of most points to form a physically reasonable surface that accounts for explicit QM solvent molecules. As the PCM cavity increases in size (e.g., SES to SAS or VDW surfaces scaled by larger α values), the QM/PCM excitation energy of all three solutes shifts to higher or lower energies depending on the relative magnitude of the ground and excite state dipole moments. This shift in energy is evidence that the QM-classical electrostatic interaction is changing with an expanding PCM cavity. The ideal α value, however, may depend on the solvent and set of atomic radii used to define the VDW surface. For the aqueous anionic solutes modeled by UFF radii studied here, an SES or an α ≈ 1.5 scaled VDW cavity is a reliable PCM approach for computing electronic excitation energies with explicit QM solvent. However, the QM/PCM calculation using an SES requires more CPU time than a scaled VDW surface. Additionally, in our other solution phase studies, we have found that using an SES cavity makes it difficult to optimize structures, which is likely due to different numbers of spheres being added when the geometry of the system changes slightly. Therefore, for improved computational efficiency, a scaled VDW cavity may be the preferred cavity of a QM/PCM

Figure 5. QM/classical excitation energies of a single solute−solvent configuration of pCT− with increasing number of QM water molecules. Classical solvent is modeled using IEF-PCM, C-PCM, or fixed point charges (MM). For IEF-PCM calculations, the molecular cavity is defined by a VDW surface using scaled (α = 1.1) UFF radii, an SES using UFF radii, a VDW surface using unscaled (α = 1.0) SMD radii, or an SES using SMD radii. For C-PCM calculations, the molecular cavity is defined by a VDW surface using scaled (α = 1.1) UFF radii or an SES using UFF radii.

Using a VDW surface, we compare unscaled (α = 1.0) SMD radii (default SMD cavity in Gaussian09) to UFF radii scaled by α = 1.1 (default IEF-PCM cavity in Gaussian09). The SMD VDW QM/PCM excitation energy converges to a value 0.022 eV higher in energy than the UFF VDW value (Figure 5). The converged SMD VDW QM/PCM excitation energy is also 10112


Article

The Journal of Physical Chemistry B 0.105 eV higher in energy than the QM/MM excitation energy computed with 250 QM water molecules. Similar to the UFF radii, the SMD radii give reasonable QM/PCM excitation energies when an SES is used to define the PCM cavity. The SMD VDW cavity may be scaled (α > 1.5) to reproduce the SMD SES QM/PCM or QM/MM excitation energies of pCT− with explicit QM solvent (see Figure S11). However, the ideal α value for SMD radii (α = 1.8) is larger than for UFF radii (α = 1.5) because the unscaled values of the SMD radii are smaller than the UFF radii. In general, QM/PCM excitation energies computed with an unscaled VDW cavity are not improved by using SMD radii instead of UFF radii, but both sets of radii give reasonable QM/PCM excitation energies when an SES or VDW surface scaled by a larger α value is used to define the PCM cavity with explicit QM solvent. Note that QM/PCM excitation energies computed using default cavity parameters of the SMD and IEF-PCM SCRF methods differ simply because the SMD and UFF radii differ. In Figure 5, we also compare the IEF-PCM and C-PCM SCRF methods, using both an α = 1.1 VDW and an SES cavity defined by UFF radii. We find that the QM/PCM excitation energy is dependent on the SCRF method when an α = 1.1 VDW surface is used. For example, the C-PCM QM/PCM excitation energy converges to a value 0.037 eV lower in energy than the IEF-PCM QM/PCM excitation energy. Using an SES to define the molecular cavity gives similar IEF-PCM and CPCM QM/PCM excitation energies for all QM region sizes. At the point of convergence (∼100 QM water molecules), the SES QM/PCM excitation energies with both SCRF methods agree to within +0.040 eV of the QM/MM value; this deviation is reduced to less than −0.005 eV with 250 QM water molecules. Thus, agreement between QM/PCM and QM/MM excitation energies with large amounts of QM solvent is more sensitive to the PCM cavity definition than the SCRF method. Comparison of Linear Response and State Specific PCM Methods. Within the PCM framework, electronic transitions can be modeled using a linear response (LR) or a state-specific (SS) approach.60,61 Using TDDFT to model vertical electronic transitions, the SS approach explicitly includes the electronic relaxation of the solvent with respect to the excited state density, whereas the LR approach uses the transition density to determine the response of the solvent. For this reason, the SS method is recommended for modeling charge-transfer transitions in polar solvents,64−66 such as the bright state of aqueous pCT− analyzed here. Using an external iteration method,58 SS PCM includes the mutual polarization of the solvent and both the ground and excited state electron densities, which may overpolarize the ground state and lead to unphysical SS-QM/PCM excitation energies.66 The results shown in Figures 2, 4, and 5 were computed using a LR-QM/PCM method. In Figure 6 we compare LRand SS-QM/PCM excitation energies for the same solute− solvent configuration of pCT− as in Figure 2a using the IEFPCM SCRF method and various PCM cavities defined by UFF radii. For clarity, the results shown in Figure 6 are within 3.35 and 3.85 eV; additional SS-QM/PCM excitation energies of pCT− computed using an SES cavity outside this energy range are shown in Figure S12a. The SS-QM/PCM excitation energies of phenolate and fluorescein− (same configurations as Figures 2b and 2c, respectively) are shown in Figures S12b and S12c, respectively. Because of the small change in dipole moment upon excitation (i.e., small Δμ0i values) for these solutes (Figure S3), the LR- and SS-QM/PCM excitation

Figure 6. LR- and SS-QM/PCM excitation energies and QM/MM excitation energies of pCT− with increasing number of QM water molecules. Classical solvent is modeled using LR or SS IEF-PCM or fixed point charges (MM). For IEF-PCM calculations, the molecular cavity is defined by a VDW surface using UFF radii scaled by α = 1.1, 1.5, or 2.0, an SES, or an SAS.

energies converge to similar values with ∼250 QM water molecules for each solute and are relatively independent of the PCM cavity definition compared to those of pCT−. Thus, additional analysis of the SS-QM/PCM excitation energies of phenolate and fluorescein− is not presented here. To investigate the discrepancy between the SS- and LR-QM/ PCM excitation energies of pCT−, we focus first on the SES cavity, which predicts physically reasonable LR-QM/PCM excitation energies with large amounts of explicit QM solvent.77 With zero QM water molecules, the SES SS-QM/PCM excitation energy is ∼0.30 eV higher in energy than the LRQM/PCM value (Figure 6). This difference is reduced to less than 0.15 eV with 100 QM water molecules and to ∼0.05 eV with 250 QM water molecules. Similarly, the deviation between the SS-QM/PCM excitation energy and the QM/MM excitation computed with 300 QM water molecules (which we take as the converged value) is significantly reduced with 250 QM water molecules (∼0.04 eV) compared to 0 QM water molecules (∼0.25 eV). This agreement is expected to improve with larger QM regions. Therefore, the SS- and LR-QM/PCM excitation energies converge to similar values that agree with QM/MM excitation energies for large amounts of explicit solvent. We next examine the effects of the PCM cavity size on the SS-QM/PCM excitation energy, starting with a continuum-only solvation model. With zero QM water molecules, an α = 1.1 VDW surface overestimates the SS-QM/PCM excitation energy relative to the converged QM/MM value by 0.23 eV. This deviation is reduced to +0.10 eV with an α = 1.5 VDW surface and +0.03 eV with an α = 2.0 VDW surface. In comparison, the LR-QM/PCM excitation energy computed by various PCM cavities is underestimated relative to the converged QM/MM value by at least 0.04 eV. Therefore, the SS PCM approach with a large enough α value improves the predicted QM/PCM excitation energy of the charge-transfer transition of aqueous pCT− studied here relative to the LR PCM approach when a continuum-only solvent model (i.e., no explicit QM solvent) is used. 10113


Article


explicit QM solvent with a PCM. We show that unphysical double counting of solvation effects can occur for some cavity definitions due to including both the QM solvent and the bulk solvent dielectric between QM molecules. In this work we have focused on how this double counting affects the TDDFT excitation energies of molecules in solution. However, the excess polarization provided when including both explicit QM solvent and the PCM dielectric near the solute is not TDDFT specific; it will be present for any electronic structure method, for both ground and excited state properties. Our results show that an SES cavity is the most reliable PCM cavity for computing electronic excitation energies with explicit QM solvent, especially if there is a large difference in the ground and excited state dipole moment. To reduce computational cost, a scaled VDW cavity may be used, but the ideal α scaling value depends on the solute−solvent configuration, the amount of explicit QM solvent, the SCRF method, and the set of radii used to define the VDW surface. Although not explored here, the ideal α scaling value will also likely depend on the identity of the solvent. For example, the discrete-continuum boundary definition is important for accurate simulations of proteins.87−89 For pCT− in water, which has a larger dipole moment in the ground state than the excited state, the QM/PCM excitation energy is physically reasonable when an SES or VDW surface scaled by α = 1.5 is used. Scaling the VDW cavity by a smaller α value leads to an overestimation of the excitation energies due to the unphysical presence of the dielectric between QM molecules. Scaling the VDW cavity by a larger α value or using an SAS cavity leads to an underestimation of the excitation energies due to the dielectric being too far away from the QM molecules. For phenolate and fluorescein−, negligible deviation between the converged QM/PCM and QM/MM excitation energies are observed for all PCM cavities studied here. These results suggest that the QM/PCM excitation energy of pCT− is more dependent on the PCM cavity than phenolate or fluorescein−, which we rationalize by the magnitude of the difference between the ground and excited state dipole moment. By comparing three SCRF methods (IEF-PCM, SMD, and C-PCM), we demonstrate that the choice of PCM cavity definition is more important than the choice of SCRF method for obtaining agreement with the converged QM/MM results; that is, using an SES to define the cavity gives similar results for all SCRF methods considered, whereas a VDW cavity gives large deviations among SCRF methods. By directly comparing SMD and IEF-PCM excitation energies computed using the same PCM cavity, we show that the deviations between the SMD and IEF-PCM methods are simply due to the different atomic radii invoked by the default cavity parameters of each method. Using IEF-PCM, we show that SMD radii require larger α values than UFF radii to reproduce SES QM/PCM excitation energies. We also show that a state-specific (SS) PCM method requires a larger PCM cavity than a linear response PCM method to obtain physically reasonable excitation energies. This result may be due to the external iteration method used to compute the SS-QM/PCM excitation energies and may be resolved by using an SS PCM method that does not reoptimize the ground state molecular orbitals during the iterative SCRF procedure, such as cLR or VEM. On the basis of these results, we suggest using an SES to compute reliable solution phase properties with mixed QM and classical solvent models. A scaled VDW surface may be used to

As QM water molecules are included in the PCM cavity, the SS-QM/PCM method using an α = 1.1 VDW surface gives unreliable excitation energies that converge to a value 0.24 eV higher in energy than the converged QM/MM excitation energy (Figure 6). Convergence improves when the VDW surface is scaled by α = 1.5, but the SS-QM/PCM value computed with 250 QM water molecules is still 0.06 eV higher in energy than the converged QM/MM excitation energy. Increasing the atomic radii scale factor to α = 2.0 reduces the SS-QM/PCM excitation energy computed with 250 QM water molecules to within 0.02 eV of the converged QM/MM excitation energy. Thus, the SS-QM/PCM excitation energy decreases with increasing α values, similar to the LR-QM/PCM excitation energies of pCT− (Figure 2a). This trend is consistent with a negative Δμ0i value (Figure S3) and the removal of unphysical PCM cavity surface points within the explicit QM solvent shell with larger α values (Figure 3a). Similarly, an SAS pushes most surface points beyond the QM solvent nuclei to give a larger PCM cavity than an SES. With 100 QM water molecules, the SS-QM/PCM excitation energy computed with an SAS cavity is only 0.03 eV higher in energy than the converged QM/MM excitation energy, compared to 0.16 eV with an SES cavity (Figure 6). Therefore, a rather large PCM cavity defined by an α = 2.0 VDW surface or an SAS is required to achieve physically reasonable SS-QM/PCM excitation energies. The origin of the discrepancy between the SS- and LR-PCM methods with smaller α values may be due to the external iteration method58 used to compute the SS-QM/PCM excitation energies. This method self-consistently determines the response of the dielectric continuum with respect to the total electron density (i.e., the ground state density and the difference between the ground and excited state densities), whereas other SS-QM/PCM methods, such as corrected linear response (cLR)57 PCM and the vertical excitation model (VEM),54 use only the difference between the ground and excited state electron densities. It has been suggested that the external iteration method overpolarizes the electron density by reoptimizing the ground state molecular orbitals during the excited state SCRF procedure.66 Here, this overpolarization results in an overestimated SS-QM/PCM excitation energy relative to the LR-QM/PCM value, which is significant for VDW cavities scaled by small α values and an SES cavity with 100 or fewer QM water molecules (Figure 6). These PCM cavities include the response of the dielectric continuum at points near the QM solute (Figure 3a and Figure S13), which may enhance the effects of the external iteration SCRF procedure on the SS-QM/PCM excitation energy. Larger PCM cavities, such as an α = 2.0 VDW or SAS cavity, remove cavity surface points near the QM solute nuclei and increase the distance between the QM electron density and the PCM cavity surface points, which may reduce the effects of the external iteration SCRF procedure on the SS-QM/PCM excitation energy. As a result, larger PCM cavities improve the agreement between SS- and LR-QM/PCM excitation energies and bring the SS-QM/PCM excitation energy in better agreement with the converged QM/MM value.

■

CONCLUSIONS Combining explicit QM solvent with a polarizable continuum model should yield accurate short and long-range solvation effects. However, the type of surface used to define the PCM molecular cavity must be chosen with care when combining 10114


Article


Molecule in Solution: A Quantum Mechanical Study. Phys. Chem. Chem. Phys. 2010, 12, 1000−1006. (6) Ma, H.; Ma, Y. Solvent Effect on Electronic Absorption, Fluorescence, and Phosphorescence of Acetone in Water: Revisited by Quantum Mechanics/Molecular Mechanics (QM/MM) Simulations. J. Chem. Phys. 2013, 138, 224505. (7) Eilmes, A. Solvatochromic Probe in Molecular Solvents: Implicit versus Explicit Solvent Model. Theor. Chem. Acc. 2014, 133, 1538. (8) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Electronic Absorption Spectra and Solvatochromic Shifts by the Vertical Excitation Model: Solvated Clusters and Molecular Dynamics Sampling. J. Phys. Chem. B 2015, 119, 958−967. (9) Zuehlsdorff, T. J.; Haynes, P. D.; Payne, M. C.; Hine, N. D. M. Predicting Solvatochromic Shifts and Colours of a Solvated Organic Dye: The Example of Nile Red. J. Chem. Phys. 2017, 146, 124504. (10) Pliego, J. R.; Riveros, J. M. The Cluster−Continuum Model for the Calculation of the Solvation Free Energy of Ionic Species. J. Phys. Chem. A 2001, 105, 7241−7247. (11) Pliego, J. R.; Riveros, J. M. Theoretical Calculation of pKa Using the Cluster−Continuum Model. J. Phys. Chem. A 2002, 106, 7434− 7439. (12) Marenich, A. V.; Ding, W.; Cramer, C. J.; Truhlar, D. G. Resolution of a Challenge for Solvation Modeling: Calculation of Dicarboxylic Acid Dissociation Constants Using Mixed Discrete− Continuum Solvation Models. J. Phys. Chem. Lett. 2012, 3, 1437− 1442. (13) Thapa, B.; Schlegel, H. B. Calculations of pKa’s and Redox Potentials of Nucleobases with Explicit Waters and Polarizable Continuum Solvation. J. Phys. Chem. A 2015, 119, 5134−5144. (14) Thapa, B.; Schlegel, H. B. Density Functional Theory Calculation of pKa’s of Thiols in Aqueous Solution Using Explicit Water Molecules and the Polarizable Continuum Model. J. Phys. Chem. A 2016, 120, 5726−5735. (15) Ho, J.; Ertem, M. Z. Calculating Free Energy Changes in Continuum Solvation Models. J. Phys. Chem. B 2016, 120, 1319−1329. (16) Korsun, O. M.; Kalugin, O. N.; Fritsky, I. O.; Prezhdo, O. V. Ion Association in Aprotic Solvents for Lithium Ion Batteries Requires Discrete−Continuum Approach: Lithium Bis(oxalato)borate in Ethylene Carbonate Based Mixtures. J. Phys. Chem. C 2016, 120, 16545− 16552. (17) da Silva, C. O.; Mennucci, B.; Vreven, T. Combining Microsolvation and Polarizable Continuum Studies: New Insights in the Rotation Mechanism of Amides in Water. J. Phys. Chem. A 2003, 107, 6630−6637. (18) Floris, F. M.; Filippi, C.; Amovilli, C. Electronic Excitations in a Dielectric Continuum Solvent with Quantum Monte Carlo: Acrolein in Water. J. Chem. Phys. 2014, 140, 034109. (19) Floris, F. M.; Filippi, C.; Amovilli, C. A Density Functional and Quantum Monte Carlo Study of Glutamic Acid in Vacuo and in a Dielectric Continuum Medium. J. Chem. Phys. 2012, 137, 075102. (20) Verdolino, V.; Cammi, R.; Munk, B. H.; Schlegel, H. B. Calculation of pKa Values of Nucleobases and the Guanine Oxidation Products Guanidinohydantoin and Spiroiminodihydantoin Using Density Functional Theory and a Polarizable Continuum Model. J. Phys. Chem. B 2008, 112, 16860−16873. (21) Jackson, V. E.; Felmy, A. R.; Dixon, D. A. Prediction of the pKa’s of Aqueous Metal Ion + 2 Complexes. J. Phys. Chem. A 2015, 119, 2926−2939. (22) Nazarski, R. B.; Justyna, K.; Leśniak, S.; Chrostowska, A. A Benefit of Using the IDSCRF- over UFF-Radii Cavities and Why Joint Correlations of NMR Chemical Shifts Can Be Advantageous: Condensed Pyridines as an IEF-PCM/GIAO/DFT Case Study. J. Phys. Chem. A 2016, 120, 9519−9528. (23) Martínez, J. M.; Pappalardo, R. R.; Sánchez Marcos, E.; Mennucci, B.; Tomasi, J. Analysis of the Opposite Solvent Effects Caused by Different Solute Cavities on the Metal−Water Distance of Monoatomic Cation Hydrates. J. Phys. Chem. B 2002, 106, 1118− 1123.

reduce computational cost, but testing is required to determine an appropriate atomic radii scale factor for a given system. Furthermore, the magnitude of the ground and excited state dipole moment difference may be used to predict the sensitivity of the QM/PCM excitation energy on the PCM cavity.

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcb.7b06693. Ground and excited state electron density difference plots, Cartesian coordinates, basis set comparison, approximate IPCM excitation energies, ground-to-excited state dipole moment difference, transition dipole moments, molecular cavity parameters and surface plots, excitation energies without a PCM, deviation in QM/ PCM excitation energies of pCT− (including configuration 3), and phenol, comparison of SMD and IEFPCM excitation energies, QM/classical excitation energies of pCT− using SMD radii, state-specific QM/PCM excitation energies, and PCM cavity surface point distances (PDF)

■

AUTHOR INFORMATION

Corresponding Author

*(C.M.I.) E-mail: [email protected]. ORCID

Makenzie R. Provorse Long: 0000-0003-4034-0248 Christine M. Isborn: 0000-0002-4905-9113 Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Acknowledgement is made to the donors of Petroleum Research Fund, administered by the American Chemical Society, for supporting the early stages of this research examining convergence, under Award Number 53674-DNI6. The later stages of this research, including comparison of various PCM cavities, were supported by the Department of Energy, Office of Basic Energy Sciences CTC and CPIMS programs, under Award Number DE-SC0014437. M.R.P.L. is grateful to the University of California Merced Chancellor’s Postdoctoral Fellowship for financial support. We acknowledge computing time on the Multi-Environmental Computer for Exploration and Discovery (MERCED) cluster which is supported by the National Science Foundation, Grant No. ACI-1429783.

■

REFERENCES

(1) Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem., Int. Ed. 2009, 48, 1198−1229. (2) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3093. (3) Mennucci, B. Modeling Environment Effects on Spectroscopies through QM/classical Models. Phys. Chem. Chem. Phys. 2013, 15, 6583. (4) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Sorting Out the Relative Contributions of Electrostatic Polarization, Dispersion, and Hydrogen Bonding to Solvatochromic Shifts on Vertical Electronic Excitation Energies. J. Chem. Theory Comput. 2010, 6, 2829−2844. (5) Pedone, A.; Bloino, J.; Monti, S.; Prampolini, G.; Barone, V. Absorption and Emission UV-Vis Spectra of the TRITC Fluorophore 10115


Article


(45) Klamt, A.; Schüürmann, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and Its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 799−805. (46) 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−260. (47) Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem. 1995, 99, 2224−2235. (48) Pye, C. C.; Ziegler, T. An Implementation of the Conductor-like Screening Model of Solvation within the Amsterdam Density Functional Package. Theor. Chem. Acc. 1999, 101, 396−408. (49) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B 2009, 113, 6378−6396. (50) Weijo, V.; Mennucci, B.; Frediani, L. Toward a General Formulation of Dispersion Effects for Solvation Continuum Models. J. Chem. Theory Comput. 2010, 6, 3358−3364. (51) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Uniform Treatment of Solute-Solvent Dispersion in the Ground and Excited Electronic States of the Solute Based on a Solvation Model with StateSpecific Polarizability. J. Chem. Theory Comput. 2013, 9, 3649−3659. (52) Cupellini, L.; Amovilli, C.; Mennucci, B. Electronic Excitations in Nonpolar Solvents: Can the Polarizable Continuum Model Accurately Reproduce Solvent Effects? J. Phys. Chem. B 2015, 119, 8984−8991. (53) Rappe, A. K.; Casewit, C. J.; Colwell, K. S.; Goddard, W. A.; Skiff, W. M. UFF, a Full Periodic Table Force Field for Molecular Mechanics and Molecular Dynamics Simulations. J. Am. Chem. Soc. 1992, 114, 10024−10035. (54) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G.; Guido, C. A.; Mennucci, B.; Scalmani, G.; Frisch, M. J. Practical Computation of Electronic Excitation in Solution: Vertical Excitation Model. Chem. Sci. 2011, 2, 2143−2161. (55) Cammi, R.; Mennucci, B.; Tomasi, J. Fast Evaluation of Geometries and Properties of Excited Molecules in Solution: A Tamm-Dancoff Model with Application to 4-Dimethylaminobenzonitrile. J. Phys. Chem. A 2000, 104, 5631−5637. (56) Cossi, M.; Barone, V. Time-Dependent Density Functional Theory for Molecules in Liquid Solutions. J. Chem. Phys. 2001, 115, 4708. (57) 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. (58) 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. (59) Improta, R.; Scalmani, G.; Frisch, M. J.; Barone, V. Toward Effective and Reliable Fluorescence Energies in Solution by a New State Specific Polarizable Continuum Model Time Dependent Density Functional Theory Approach. J. Chem. Phys. 2007, 127, 074504. (60) 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. (61) Corni, S.; Cammi, R.; Mennucci, B.; Tomasi, J. Electronic Excitation Energies of Molecules in Solution within Continuum Solvation Models: Investigating the Discrepancy between StateSpecific and Linear-Response Methods. J. Chem. Phys. 2005, 123, 134512. (62) Runge, E.; Gross, E. K. U. Density-Functional Theory for TimeDependent Systems. Phys. Rev. Lett. 1984, 52, 997−1000.

(24) Cossi, M.; Barone, V.; Cammi, R.; Tomasi, J. Ab Initio Study of Solvated Molecules: A New Implementation of the Polarizable Continuum Model. Chem. Phys. Lett. 1996, 255, 327−335. (25) Preat, J.; Loos, P.-F.; Assfeld, X.; Jacquemin, D.; Perpète, E. A. DFT and TD-DFT Investigation of IR and UV Spectra of Solvated Molecules: Comparison of Two SCRF Continuum Models. Int. J. Quantum Chem. 2007, 107, 574−585. (26) Psciuk, B. T.; Lord, R. L.; Munk, B. H.; Schlegel, H. B. Theoretical Determination of One-Electron Oxidation Potentials for Nucleic Acid Bases. J. Chem. Theory Comput. 2012, 8, 5107−5123. (27) de Lima, G. F.; Duarte, H. A.; Pliego, J. R. Dynamical Discrete/ Continuum Linear Response Shells Theory of Solvation: Convergence Test for NH4+ and OH− Ions in Water Solution Using DFT and DFTB Methods. J. Phys. Chem. B 2010, 114, 15941−15947. (28) Connolly, M. L. Analytical Molecular Surface Calculation. J. Appl. Crystallogr. 1983, 16, 548−558. (29) Pascual-ahuir, J. L.; Silla, E.; Tunon, I. GEPOL: An Improved Description of Molecular Surfaces. III. A New Algorithm for the Computation of a Solvent-Excluding Surface. J. Comput. Chem. 1994, 15, 1127−1138. (30) Cossi, M.; Mennucci, B.; Cammi, R. Analytical First Derivatives of Molecular Surfaces with Respect to Nuclear Coordinates. J. Comput. Chem. 1996, 17, 57−73. (31) Connolly, M. Solvent-Accessible Surfaces of Proteins and Nucleic Acids. Science 1983, 221, 709−713. (32) Foresman, J. B.; Keith, T. A.; Wiberg, K. B.; Snoonian, J.; Frisch, M. J. Solvent Effects. 5. Influence of Cavity Shape, Truncation of Electrostatics, and Electron Correlation on Ab Initio Reaction Field Calculations. J. Phys. Chem. 1996, 100, 16098−16104. (33) Guareschi, R.; Floris, F. M.; Amovilli, C.; Filippi, C. Solvent Effects on Excited-State Structures: A Quantum Monte Carlo and Density Functional Study. J. Chem. Theory Comput. 2014, 10, 5528− 5537. (34) Barone, V.; Cossi, M.; Tomasi, J. A New Definition of Cavities for the Computation of Solvation Free Energies by the Polarizable Continuum Model. J. Chem. Phys. 1997, 107, 3210. (35) Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161−2200. (36) Fattebert, J.-L.; Gygi, F. Density Functional Theory for Efficient Ab Initio Molecular Dynamics Simulations in Solution. J. Comput. Chem. 2002, 23, 662−666. (37) Fattebert, J.-L.; Gygi, F. First-Principles Molecular Dynamics Simulations in a Continuum Solvent. Int. J. Quantum Chem. 2003, 93, 139−147. (38) Scherlis, D. A.; Fattebert, J.-L.; Gygi, F.; Cococcioni, M.; Marzari, N. A Unified Electrostatic and Cavitation Model for FirstPrinciples Molecular Dynamics in Solution. J. Chem. Phys. 2006, 124, 074103. (39) Dziedzic, J.; Helal, H. H.; Skylaris, C.-K.; Mostofi, A. A.; Payne, M. C. Minimal Parameter Implicit Solvent Model for Ab Initio Electronic-Structure Calculations. Europhys. Lett. 2011, 95, 43001. (40) Tomasi, J.; Mennucci, B.; Cancès, E. The IEF Version of the PCM Solvation Method: An Overview of a New Method Addressed to Study Molecular Solutes at the QM Ab Initio Level. J. Mol. Struct.: THEOCHEM 1999, 464, 211−226. (41) 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.; et al. Gaussian 09, revision D.01; Gaussian, Inc.: Wallingford, CT, 2013. (42) Scalmani, G.; Frisch, M. J. Continuous Surface Charge Polarizable Continuum Models of Solvation. I. General Formalism. J. Chem. Phys. 2010, 132, 114110. (43) 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−2001. (44) 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−681. 10116


Article

The Journal of Physical Chemistry B (63) Casida, M. E. Time-Dependent Density Functional Response Theory for Molecules. In Recent Advances in Density Functional Methods; Chong, D. P., Ed.; Recent Advances in Computational Chemistry; World Scientific: 1995; Vol. 1, pp 155−192. (64) Improta, R.; Barone, V. PCM/TD-DFT Study of the Two Lowest Excited States of Uracil Derivatives in Solution: The Effect of the Functional and of the Cavity Model. J. Mol. Struct.: THEOCHEM 2009, 914, 87−93. (65) Cerezo, J.; Petrone, A.; Ferrer, F. J. A.; Donati, G.; Santoro, F.; Improta, R.; Rega, N. Electronic Spectroscopy of a Solvatochromic Dye in Water: Comparison of Static Cluster/Implicit and Dynamical/ Explicit Solvent Models on Structures and Energies. Theor. Chem. Acc. 2016, 135, 263. (66) 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−5790. (67) Hu, L.; Söderhjelm, P.; Ryde, U. On the Convergence of QM/ MM Energies. J. Chem. Theory Comput. 2011, 7, 761−777. (68) Liao, R.-Z.; Thiel, W. Convergence in the QM-Only and QM/ MM Modeling of Enzymatic Reactions: A Case Study for Acetylene Hydratase. J. Comput. Chem. 2013, 34, 2389−2397. (69) Sumowski, C. V.; Ochsenfeld, C. A Convergence Study of QM/ MM Isomerization Energies with the Selected Size of the QM Region for Peptidic Systems. J. Phys. Chem. A 2009, 113, 11734−11741. (70) Meier, K.; Thiel, W.; van Gunsteren, W. F. On the Effect of a Variation of the Force Field, Spatial Boundary Condition and Size of the QM Region in QM/MM MD Simulations. J. Comput. Chem. 2012, 33, 363−378. (71) Flaig, D.; Beer, M.; Ochsenfeld, C. Convergence of Electronic Structure with the Size of the QM Region: Example of QM/MM NMR Shieldings. J. Chem. Theory Comput. 2012, 8, 2260−2271. (72) Retegan, M.; Neese, F.; Pantazis, D. A. Convergence of QM/ MM and Cluster Models for the Spectroscopic Properties of the Oxygen-Evolving Complex in Photosystem II. J. Chem. Theory Comput. 2013, 9, 3832−3842. (73) Isborn, C. M.; Götz, A. W.; Clark, M. A.; Walker, R. C.; Martínez, T. J. Electronic Absorption Spectra from MM and Ab Initio QM/MM Molecular Dynamics: Environmental Effects on the Absorption Spectrum of Photoactive Yellow Protein. J. Chem. Theory Comput. 2012, 8, 5092−5106. (74) Murugan, N. A.; Jha, P. C.; Rinkevicius, Z.; Ruud, K.; Ågren, H. Solvatochromic Shift of Phenol Blue in Water from a Combined Car− Parrinello Molecular Dynamics Hybrid Quantum Mechanics-Molecular Mechanics and ZINDO Approach. J. Chem. Phys. 2010, 132, 234508. (75) Zuehlsdorff, T. J.; Haynes, P. D.; Hanke, F.; Payne, M. C.; Hine, N. D. M. Solvent Effects on Electronic Excitations of an Organic Chromophore. J. Chem. Theory Comput. 2016, 12, 1853−1861. (76) Manzoni, V.; Lyra, M. L.; Coutinho, K.; Canuto, S. Comparison of Polarizable Continuum Model and Quantum Mechanics/Molecular Mechanics Solute Electronic Polarization: Study of the Optical and Magnetic Properties of Diazines in Water. J. Chem. Phys. 2011, 135, 144103. (77) Provorse, M. R.; Peev, T.; Xiong, C.; Isborn, C. M. Convergence of Excitation Energies in Mixed Quantum and Classical Solvent: Comparison of Continuum and Point Charge Models. J. Phys. Chem. B 2016, 120, 12148−12159. (78) Brancato, G.; Barone, V.; Rega, N. Theoretical Modeling of Spectroscopic Properties of Molecules in Solution: Toward an Effective Dynamical Discrete/Continuum Approach. Theor. Chem. Acc. 2007, 117, 1001−1015. (79) Martínez-Fernández, L.; Pepino, A. J.; Segarra-Martí, J.; Banyasz, A.; Garavelli, M.; Improta, R. Computing the Absorption and Emission Spectra of 5-Methylcytidine in Different Solvents: A Test-Case for Different Solvation Models. J. Chem. Theory Comput. 2016, 12, 4430− 4439.

(80) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926. (81) Hirata, S.; Head-Gordon, M. Time-Dependent Density Functional Theory within the Tamm−Dancoff Approximation. Chem. Phys. Lett. 1999, 314, 291−299. (82) Chai, J.-D.; Head-Gordon, M. Systematic Optimization of LongRange Corrected Hybrid Density Functionals. J. Chem. Phys. 2008, 128, 084106. (83) Nielsen, I. B.; Boyé-Péronne, S.; El Ghazaly, M. O. A.; Kristensen, M. B.; Brøndsted Nielsen, S.; Andersen, L. H. Absorption Spectra of Photoactive Yellow Protein Chromophores in Vacuum. Biophys. J. 2005, 89 (4), 2597−2604. (84) He, Z.; Martin, C. H.; Birge, R.; Freed, K. F. Theoretical Studies on Excited States of a Phenolate Anion in the Environment of Photoactive Yellow Protein. J. Phys. Chem. A 2000, 104, 2939−2952. (85) Wang, Y.; Li, H. Excited State Geometry of Photoactive Yellow Protein Chromophore: A Combined Conductorlike Polarizable Continuum Model and Time-Dependent Density Functional Study. J. Chem. Phys. 2010, 133, 034108. (86) Zuev, D.; Bravaya, K. B.; Crawford, T. D.; Lindh, R.; Krylov, A. I. Electronic Structure of the Two Isomers of the Anionic Form of PCoumaric Acid Chromophore. J. Chem. Phys. 2011, 134, 034310. (87) Welch, W. R. W.; Kubelka, J. DFT-Based Simulations of Amide I′ IR Spectra of a Small Protein in Solution Using Empirical Electrostatic Map with a Continuum Solvent Model. J. Phys. Chem. B 2012, 116, 10739−10747. (88) Lipparini, F.; Lagardère, L.; Raynaud, C.; Stamm, B.; Cancès, E.; Mennucci, B.; Schnieders, M.; Ren, P.; Maday, Y.; Piquemal, J.-P. Polarizable Molecular Dynamics in a Polarizable Continuum Solvent. J. Chem. Theory Comput. 2015, 11, 623−634. (89) Setny, P.; Dudek, A. Explicit Solvent Hydration Benchmark for Proteins with Application to the PBSA Method. J. Chem. Theory Comput. 2017, 13, 2762−2776.

10117


Combining Explicit Quantum Solvent with a Polarizable Continuum

Recommend Documents