Combining Explicit Quantum Solvent with a Polarizable Continuum

Oct 9, 2017 - A promising approach for accurately modeling both short-range and long-range solvation effects is to combine explicit quantum mechanical...
0 downloads 7 Views 2MB Size
Subscriber access provided by University of Florida | Smathers Libraries

Article

Combining Explicit Quantum Solvent with a Polarizable Continuum Model Makenzie Rae Provorse, and Christine Marie Isborn J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b06693 • Publication Date (Web): 09 Oct 2017 Downloaded from http://pubs.acs.org on October 18, 2017

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

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

Page 1 of 51

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

The Journal of Physical Chemistry

Combining Explicit Quantum Solvent with a Polarizable Continuum Model Makenzie R. Provorse Longa and Christine M. Isbornb,* a

Department of Chemistry, University of Central Arkansas, Conway, Arkansas 72035, United

States b

Chemistry and Chemical Biology, University of California Merced, Merced, California 95343,

United States Corresponding Author * E-mail: [email protected]

ACS Paragon Plus Environment

1

The Journal of Physical Chemistry

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

Page 2 of 51

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 selfconsistent 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. 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.

ACS Paragon Plus Environment

2

Page 3 of 51

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

The Journal of Physical Chemistry

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 solute and solvent 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 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,

ACS Paragon Plus Environment

3

The Journal of Physical Chemistry

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

Page 4 of 51

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 water,4–8 but also for non-polar solvents that interact with the solute via π-stacking.9 A discrete-continuum solvation model is also recommended to accurately predict solvation free energies 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 extremely 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,24,25 Many molecular-based cavities have been proposed and a schematic of common PCM molecular cavities is shown in Figure 1. The most common 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

ACS Paragon Plus Environment

4

Page 5 of 51

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

The Journal of Physical Chemistry

(i.e. VDW surface).2 The atomic radii used to define the VDW surface can be scaled by a parameter (a) 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,26–28 but is not typically used with mixed discrete-continuum solvent models. One exception is the work of Pliego and co-workers.29 They report that a ≥ 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 a 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 a value removes these regions and recovers a physically reasonable solvation free energy.

a)

van der Waals (VDW)

b)

Scaled VDW

c)

Solvent Excluded Surface (SES)

d)

Solvent Accessible Surface (SAS)

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.

ACS Paragon Plus Environment

5

The Journal of Physical Chemistry

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

Page 6 of 51

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.30 Using the GEPOL algorithm,31 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.32 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).33 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,32 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 self-consistent IPCM methods use an isosurface based on the electron density to define the PCM cavity.34 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,35 but isodensity methods can often be plagued by convergence issues.36,37 The formulation of Fattebert and Gygi38,39 is the basis for other isodensity-based continuum solvation models.40,41

ACS Paragon Plus Environment

6

Page 7 of 51

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

The Journal of Physical Chemistry

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 QMPCM 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 (IEFPCM),42 which, in its implementation in Gaussian09,43 uses a continuous surface charge formalism44 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,45,46 COSMO,47 and others.48–50 In addition to electrostatic interactions, non-electrostatic contributions such as dispersion and cavitation energies may be included in continuum solvation models.4,51–54 The SMD SCRF method uses semiempirical functions to model the non-electrostatic interactions in addition to the electrostatic interactions modeled by the IEF-PCM.51 Other SCRF methods include nonelectrostatic terms in the effective Hamiltonian to explicitly account for non-electrostatic effects during the SCRF procedure, thus affecting the final QM density.52,54 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 the VDW surface, including the universal force field (UFF)55 radii and the intrinsic Coulomb radii of the SMD method,51 which we refer to as SMD radii. The SMD radii are parameterized for different solvents to reproduce experimental solvation free energies and are recommended for computing pKa values12 and vertical electronic excitation energies.4,7,56

ACS Paragon Plus Environment

7

The Journal of Physical Chemistry

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

Page 8 of 51

Excited state calculations within the PCM framework can use a linear response (LR)57,58 or state-specific (SS)56,59–61 approach.62,63 Using time-dependent density functional theory (TDDFT),64,65 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 chargetransfer 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 charge-transfer transitions.66–68 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,69–78 Our recent work compared the convergence of QM/MM and QM/PCM electronic excitation energies of aqueous solutes.79 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 a = 1.1) to obtain physically reasonable QM/PCM

ACS Paragon Plus Environment

8

Page 9 of 51

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

The Journal of Physical Chemistry

excitation energies.79 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.79 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-p-coumarate (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.

ACS Paragon Plus Environment

9

The Journal of Physical Chemistry

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

Page 10 of 51

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,67,80,81 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 of pCT– and phenolate in aqueous solvent may be found in Ref. 79. Each solute was solvated by a pre-equilibrated 32 Å sphere of TIP3P waters82 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 program43 was used for all QM/classical calculations. All excitation energies were calculated using linear response TDDFT64,65 within the Tamm-Dancoff approximation.83 The long-range corrected hybrid functional w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.84 Although dispersion interactions between the QM solute and QM solvent may be important for computing accurate solution-phase 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

ACS Paragon Plus Environment

10

Page 11 of 51

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

The Journal of Physical Chemistry

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 waters82 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 non-equilibrium linear response PCM (LR-QM/PCM),57,58 which is the default SCRF method for TDDFT calculations in Gaussian09. Non-equilibrium state-specific QM/PCM excitation energies (SSQM/PCM) were computed by the external iteration method developed by Improta et al.60,61 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 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-PCM42 and C-PCM45,46 SCRF methods, the VDW surface was defined using UFF55 radii scaled by a = 1.1 (the default in Gaussian09). The SMD method51 defines the VDW surface using unscaled (a = 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

ACS Paragon Plus Environment

11

The Journal of Physical Chemistry

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

Page 12 of 51

calculations were computed using the default dielectric constants for water (𝜖% = 78.355300, 𝜖& = 1.777849). We tested the IPCM34 SCRF method for pCT– with up to 250 QM water molecules. The isodensity surface was computed using 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 would could 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 this 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 work79 showed that the QM/MM excitation energy is not affected by the choice of non-polarizable water model (i.e. TIP3P versus SPC/E); here, we use the TIP3P water model.

ACS Paragon Plus Environment

12

Page 13 of 51

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

The Journal of Physical Chemistry

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 𝜇( is larger in magnitude than the excited state dipole moment 𝜇) , yielding a negative value for the ground-excited state dipole moment difference ∆𝜇() = 𝜇) − 𝜇( . We expect the dielectric medium to interact more strongly with the larger ground state dipole moment than the excited state dipole 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 ∆𝜇() (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 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 ∆𝜇() . 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 ∆𝜇() 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 of using

ACS Paragon Plus Environment

13

The Journal of Physical Chemistry

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

Page 14 of 51

three common PCM cavity definitions: a VDW surface defined by UFF radii (scaling by a = 1.1 is the default in Gaussian09),43 an SES, and an SAS. For the three configurations examined in Figure 2, the negative ∆𝜇() 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.85 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,75,86–88 but we report these values for a single configuration in order to examine the effects of explicit solvent combined with a PCM. The solvatochromic shifts going from vacuum to the PCM with the a = 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 non-zero 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

ACS Paragon Plus Environment

14

Page 15 of 51

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

The Journal of Physical Chemistry

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 cavities (smaller a values) give higher excitation energies and larger cavities (larger a 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 a values.66 Similarly, we rationalize that smaller a 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 a values. The increased excitation energies with smaller a values are consistent with negative ∆𝜇() 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 a = 1.1 VDW, SES, and SAS cavities, as well as the QM/MM excitation energy, converge to within ±0.050 eV of the values computed with the largest QM region (300 QM water molecules) with ~100 QM water molecules (Figure 2a). However, the a = 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.

ACS Paragon Plus Environment

15

The Journal of Physical Chemistry

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

Page 16 of 51

For pCT–, we see that smaller a values (a = 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 ∆𝜇() value (Figure S3). 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 ∆𝜇() 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.

ACS Paragon Plus Environment

16

Page 17 of 51

PCM VDW 1.1 PCM VDW 1.2 PCM VDW 1.4 PCM VDW 1.5 PCM VDW 1.6 PCM VDW 1.8 PCM VDW 2.0 PCM VDW 2.5 PCM VDW 3.0 PCM SES PCM SAS MM

4.0 3.9

1 3.8 2

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

Excitation Energy [eV]

3.7

3.8

3.7 3.6 100 150 200 250 Number3.6 of QM Water Molecules 3.5

PCM VDW 1.1 PCM VDW 1.2 PCM VDW 1.4 PCM VDW 1.5 PCM VDW 1.6 PCM VDW 1.8 PCM VDW 2.0 PCM VDW 2.5 PCM VDW 3.0 PCM SES 300 PCM SAS MM

(a)

3.5 3.4 3.4 0 5.2

0

50 50 100 100 150 150 200 200 250 250 300 300 (b) Number of QM Molecules Number of Water QM Water Molecules

5.0 Excitation Energy [eV]

Excitation Energy [eV]

4.8 4.6 4.4

5.05 5.00 4.95 15 50 100 150 200 250 Number of QM Water Molecules

4.2

5.0

0

50

100 150 200 250 Number of QM Water Molecules

300

(c)

4.8 Excitation Energy [eV]

50

3.8

3.9

Excitation Energy [eV]

0

Excitation Energy [eV]

4.0

The Journal of Physical Chemistry

4.6 4.4 4.2

4.95 4.90 4.85 4.80

4.0 0

50

15 50 100 150 200 250 Number of QM Water Molecules 100 150 200 250 300

Number of QM Water Molecules

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.

ACS Paragon Plus Environment

17

0.07 0.06 0.05

VDW 1.1 0.04 VDW 1.2 0.03 VDW 1.4 VDW 1.5 0.02 VDW 1.6 VDW 1.8 0.01 VDW 2.0 0.00 VDW 2.5 VDW 3.0 0 SES SAS Explicit

0.08 0.07 0.02 0.06 0.05 0.04 0.01 0.03 0.02 0.01 0.00 0.00 00

1

22

0.03

0.50 0.08 0.48 0.46 0.44 0.42 0.40 0.38 0.06 0.36 0.34 0.32 0.30 0.28 0.26 1 2 3 4 5 6 0.24 7 0.04 0.22 Distance [Angstroms] 0.20 0.18 0.16 0.14 0.12 0.02 0.10 0.08 0.06 0.04 0.02 0.00 0.00 3 4 5 6 7 8 3 4 5 6 7 8 9 10 Distance Distance[Angstroms] [Angstroms] 0.08 (b)

(a)

0.06 0.02 0.04 0.01 0.02

0.00 0

1

2

0.03

3 4 5 6 7 Distance [Angstroms]

8

9

10

(c)

0.00

0.08

0.06 0.02 0.04 0.01 0.02

0.00 0

1

2

3 4 5 6 7 Distance [Angstroms] Distance [Å]

8

9

10

0.00

8

0.48 0.46 0.44 0.42 0.40 0.38 0.36 0.34 0.32 0.30 0.28 0.26 0.24 0.22 0.20 0.18 0.16 0.14 0.12 0.10 0.08 0.06 0.04 0.02 0.00

Page 18 of 51

Normalized Count of Solvent Atoms

0.08

Normalized Count of Solvent Atoms Normalized Count of Solvent Atoms

0.09

The Journal of Physical Chemistry 0.50

Normalized Count of QM Solvent Atoms Nuclei Normalized Count of Solvent

0.10 0.03 0.09

VDW 1.1 VDW 1.2 VDW 1.4 VDW 1.5 VDW 1.6 VDW 1.8 VDW 2.0 VDW 2.5 VDW 3.0 SES SAS Explicit

Normalized Count of Solvent Atoms

Normalized Count of Surface Points

Normalized Count of Surface Points Normalized Count ofCount PCM Cavity Surface Points Normalized Normalized Count of Surface Points Normalized of Surface Points Count of Surface Points

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

0.10

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 y-axis.

ACS Paragon Plus Environment

18

Page 19 of 51

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

The Journal of Physical Chemistry

We next examine how scaling the VDW cavity, using UFF radii scaled by a = 1.1 – 3.0, changes the excitation energy. Increasing the a value increases the volume of the molecular cavity and, up to a ≈ 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 a 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 a value for most QM region sizes (Figure 2a). In fact, for a given QM region size, QM/PCM excitation energies computed with an increasing a 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 ∆𝜇() value (Figure S3); as the a value increases, the weaker QM-PCM interaction destabilizes the excited state relative to the ground state, increasing the excitation energy with larger a values (Figure 2a). The QM/PCM excitation energies computed with an a = 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 a = 1.5 VDW and an SES (Figure 3a). The a = 2.0 VDW QM/PCM excitation energies agree with those computed with an SAS for all QM region sizes (Figure 2a) because an a = 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 a values (a = 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 a values, an

ACS Paragon Plus Environment

19

The Journal of Physical Chemistry

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

Page 20 of 51

artificial vacuum is formed between the QM solvent molecules furthest from the solute and most of the cavity surface points. As a result, the a = 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 a value required more explicit QM water molecules to converge the solvation free energy.29 With 250 QM water molecules, the a = 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 a = 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 a = 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 a = 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 ∆𝜇() value than pCT– (Figure S3). With 15 QM water molecules, the QM/PCM excitation energy of phenolate computed with an a = 1.1 VDW surface

ACS Paragon Plus Environment

20

Page 21 of 51

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

The Journal of Physical Chemistry

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 a = 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 a 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 ∆𝜇() value (Figure S3). Also, larger a 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 a 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 ∆𝜇() 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

ACS Paragon Plus Environment

21

The Journal of Physical Chemistry

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

Page 22 of 51

located between QM solute and solvent nuclei for scaled VDW surfaces with small a values (Figure 3b). We next analyze the trends in the QM/classical excitation energies of fluorescein–, which has a small, negative ∆𝜇() 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 a value, which is consistent with larger distances between the QM nuclei and the PCM cavity surface points (Figure 3c) and a negative ∆𝜇() 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 ∆𝜇() (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 ∆𝜇() 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

ACS Paragon Plus Environment

22

Page 23 of 51

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

The Journal of Physical Chemistry

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 both an SES and SAS and most a 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 a = 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 a value that minimizes the deviation is a = 1.5; however, a 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 a 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.

ACS Paragon Plus Environment

23

0.15

0.15

0.10

0.10

0.05

0.05

0.00

0.00

Configuration 1 Configuration 2 Configuration 3 Configuration 4 Configuration 5

'()

-0.05 -0.05 -0.10 -0.10

Deviation in Excitation Energy (eV)

%% [eV] in Excitation Energy (eV) #$% Deviation in Excitation Energy (eV) of QM/PCM Deviation Excitation in Excitation Energy!(eV) Deviation Energy − !Deviation

(a)0.15

Page 24 of 51

(d)

0.10 0.05 0.00 -0.05 -0.10

0.10 0.05 0.00 -0.05 -0.10 1.0 0.15

1.5 2.0 2.5 VDW Radii Scale Factor (α)

0.10 0.05 0.00 -0.05 -0.10 1.0

1.5 2.0 2.5 Atomic radiiScale scaleFactor factor (*) VDW Radii (α)

Deviation in Excitation Energy (eV)

Deviation in Excitation Energy (eV)

1.0 1.5 1.5 2.0 2.0 2.5 2.5 3.0 3.0 1.0 0.151.0 1.5 (b)0.15 VDW Radii Scale (e) Factor (α) VDWScale RadiiFactor Scale (α) Factor (α) VDW Radii

"

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

Deviation in Excitation Energy (eV)

The Journal of Physical Chemistry

0.10 0.05 0.00 -0.05 -0.10 1.0 1.5 3.0 (c)0.15 (f) VDW Radii Scale Factor (α)

0.10 0.05 0.00 -0.05 -0.10

SAS 3.0 SES 1.0 1.5 VDW Radii Scale Factor (α)

/01 Figure 4. Deviation of QM/PCM excitation energies 𝜔. from QM/MM excitation energies 11 computed with 250 QM water molecules 𝜔23( for five solute-solvent configurations of (a)

pCT– (N = 100 QM water molecules), (b) phenolate (N = 50 QM water molecules), and (c) fluorescein anion (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 (a). 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.

ACS Paragon Plus Environment

24

Page 25 of 51

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

The Journal of Physical Chemistry

Table 1. Difference between ground (𝜇( ) and excited state (𝜇) ) dipole moments (𝜇) − 𝜇( ) for three solutes and surrounding QM solvent in debye (D).a pCT– b

phenolatec

fluorescein– b

Configuration 1

-4.14

+1.10

-0.59

Configuration 2

-2.35

-1.43

-2.43

Configuration 3

-5.17

-0.62

-0.06

Configuration 4

-5.87

-0.54

-0.26

Configuration 5

-5.82

-1.09

+0.37

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

Unlike for pCT–, the deviations of the QM/PCM excitation energies for scaled VDW cavities of all five configurations of phenolate (Figure 4b) and fluorescein– (Figure 4c) are near or less than ±0.050 eV for all a values tested. Some configurations (e.g. Configuration 3 of fluorescein–) show negligible changes in QM/PCM excitation energy for different a 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 ∆𝜇() value (Table 1). Only these two configurations show an increase in QM/PCM excitation energy with increasing a 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 a values and underpolarization by cavities with larger a values. Therefore, the dependence of the QM/PCM excitation energy on the PCM cavity may be predicted by the sign and magnitude of ∆𝜇() . In addition, solute-solvent configurations with larger ∆𝜇() values are

ACS Paragon Plus Environment

25

The Journal of Physical Chemistry

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

Page 26 of 51

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 a = 1.4 VDW cavities, whereas a = 1.5 minimizes the deviation for pCT–. These a 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 a = 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 a = 1.8 or 2.0 give similar QM/PCM excitation energies as SAS cavities. These results for five solutesolvent 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 50 QM water molecules, an SES cavity and scaled VDW cavities with a = 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 a value is used. PCM cavities defined by an SES, SAS, or scaled VDW surface with a large a value (a ~

ACS Paragon Plus Environment

26

Page 27 of 51

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

The Journal of Physical Chemistry

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 a 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 a 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 a ≈ 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 calculation with large amounts of explicit QM solvent, but tests should be performed to determine the appropriate a 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.

ACS Paragon Plus Environment

27

The Journal of Physical Chemistry

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

Page 28 of 51

Comparison of IEF-PCM, SMD, and C-PCM methods. The results presented in Figures 2 - 4 were computed using IEF-PCM with cavities defined by UFF radii (the default SCRF method and atomic radii in Gaussian09). The SMD method includes a non-electrostatic 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.51 The non-electrostatic term is added to the total energy determined by the ground state calculation and does not contribute to the excited state calculation. Therefore, excitation energies computed by the SMD method are equivalent to IEFPCM 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 IEF-PCM 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. Using a VDW surface, we compare unscaled (a = 1.0) SMD radii (default SMD cavity in Gaussian09) to UFF radii scaled by a = 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 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 (a > 1.5) to

ACS Paragon Plus Environment

28

Page 29 of 51

reproduce the SMD SES QM/PCM or QM/MM excitation energies of pCT– with explicit QM solvent (see Figure S11). However, the ideal a value for SMD radii (a = 1.8) is larger than for UFF radii (a = 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 a 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.

IEF-PCM VDW UFF IEF-PCM SES UFF IEF-PCM VDW SMD IEF-PCM SES SMD C-PCM VDW UFF C-PCM SES UFF MM

3.8

Excitation Energy [eV]

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

The Journal of Physical Chemistry

3.7 3.6 3.5 3.4 0

50 100 150 200 Number of QM Water Molecules

250

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, CPCM, or fixed point charges (MM). For IEF-PCM calculations, the molecular cavity is defined by a VDW surface using scaled (a = 1.1) UFF radii, an SES using UFF radii, a VDW surface using unscaled (a = 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 (a = 1.1) UFF radii or an SES using UFF radii.

ACS Paragon Plus Environment

29

The Journal of Physical Chemistry

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

Page 30 of 51

In Figure 5, we also compare the IEF-PCM and C-PCM SCRF methods, using both an a = 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 a = 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 C-PCM 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.62,63 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 electron density to determine the response of the solvent. For this reason, the SS method is recommended for modeling charge-transfer transitions in polar solvents,66–68 such as the bright state of aqueous pCT– analyzed here. Using an external iteration method,60 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.68

ACS Paragon Plus Environment

30

Page 31 of 51

The results shown in Figures 2, 4, and 5 were computed using a LR-QM/PCM method. In Figure 6 we compare LR- and SS-QM/PCM excitation energies for the same solute-solvent configuration of pCT– as in Figure 2a using the IEF-PCM 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 S11b and S11c, respectively. Due to the small change in dipole moment upon excitation (i.e. small ∆𝜇() values) for these solutes (Figure S3), the LR- and SS-QM/PCM excitation 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.

LR VDW 1.1 LR VDW 1.5 1.1 1.1 LR LRVDW VDW2.0 1.1 1.5 LR VDW 1.5 LR VDW 1.1 SS 1.5 LR VDW 1.1 2.0 SS LR VDW 2.0 SS LRVDW VDW 1.5 2.0 1.1 LR 1.5 1.1 SS SES SS 1.1 LRVDW VDW 2.0 SS 1.1 1.5 LR 2.0 SS 1.5 LR SES 1.5 SS 1.1 1.5 LRVDW VDWMM 2.0 SS 1.1 SSVDW VDW1.5 1.5 1.1 4.4 SS SS VDW 1.5 4.2

LRVDW VDW2.0 1.1 SS LRVDW VDW 1.5 SS SES SS 2.0 SS 2.0 LRVDW VDW 2.0 LR SES SS SS SS SES SSVDW VDW 1.1 2.0 MM SSSES SES SS 2.0 LR LR SS SS VDW 1.5 SS LRSES SES 2.0 SES MM MM LRSES SES MM SS LR MM LR MM SES SS SAS MM

SS VDW 2.0 SS SES LR SES MM

Excitation Energy [eV]

LR VDW 1.1 LR VDW 1.5 LR VDW 2.0 SS VDW 1.1 SS VDW 1.5 4.4 4.4 4.4 4.2 3.8 4.4 4.4 4.2 4.4 3.8 4.0 4.2 4.0 4.2 4.4 4.2 4.0 4.2 3.7 3.8 4.0 3.8 4.0 4.2 4.0 3.8 4.0 3.6 3.8 3.6 3.8 3.6 4.0 3.6 3.8 3.6 3.8 3.4 3.4 3.6 3.6 3.8 3.6 0 3.4 3.6 0 50 250 50 100 150 200 100 250 150 300 200 3.4 3.4 3.6 3.5 0 3.4 50 100 150 200 250 300 Number of QM Water Molecules Number of QM Water Molecules 3.4 0 3.4 50 100 0 50 100 150 150 200 200 250 250 300 300 3.4 00 Number of QM150 Water Molecules 50 100 150 200 250 250 300 300 100 200 0 50 100 150 50 200 250 300 Number of QM Water Molecules Number of QM Water Molecules 50 100 150 200 250 300 3.4 0 Number QMWater WaterMolecules Molecules Number ofofQM Number of QM Water Molecules Number of QM Water Molecules 0 50 100 150 200 250 300 Number of QM Water Molecules Excitation Excitation Excitation Excitation Energy Energy Energy Energy [eV] [eV] [eV] [eV] Excitation Energy [eV] Excitation Excitation Energy Energy [eV] [eV] Excitation Energy [eV] Excitation Energy [eV]

Excitation Energy [eV]

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

The Journal of Physical Chemistry

300

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 IEF-PCM

ACS Paragon Plus Environment

31

The Journal of Physical Chemistry

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

Page 32 of 51

(PCM) or fixed point charges (MM). For IEF-PCM calculations, the molecular cavity is defined by a VDW surface using UFF radii scaled by a = 1.1, 1.5, or 2.0, an SES, or an SAS.

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.79 With zero QM water molecules, the SES SS-QM/PCM excitation energy is ~0.30 eV higher in energy than the LR-QM/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 a = 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 a = 1.5 VDW surface and +0.03 eV with an a = 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 improves the predicted QM/PCM

ACS Paragon Plus Environment

32

Page 33 of 51

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

The Journal of Physical Chemistry

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. As QM water molecules are included in the PCM cavity, the SS-QM/PCM method using an a = 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 a = 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 a = 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 a values, similar to the LR-QM/PCM excitation energies of pCT– (Figure 2a). This trend is consistent with a negative ∆𝜇() value (Figure S3) and the removal of unphysical PCM cavity surface points within the explicit QM solvent shell with larger a 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 a = 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 a values may be due to the external iteration method60 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

ACS Paragon Plus Environment

33

The Journal of Physical Chemistry

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

Page 34 of 51

ground and excited state densities), whereas other SS-QM/PCM methods, such as corrected linear response (cLR)59 PCM and the vertical excitation model (VEM),56 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 re-optimizing the ground state molecular orbitals during the excited state SCRF procedure.68 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 a 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 a = 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 LRQM/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 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

ACS Paragon Plus Environment

34

Page 35 of 51

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

The Journal of Physical Chemistry

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 a 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 a scaling value will also likely depend on the identity of the solvent. The discrete-continuum boundary definition is also important for accurate simulations of proteins.89–91 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

a = 1.5 is used. Scaling the VDW cavity by a smaller a 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 a 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

ACS Paragon Plus Environment

35

The Journal of Physical Chemistry

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

Page 36 of 51

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 IEFPCM 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 a values than UFF radii to reproduce SES QM/PCM excitation energies. We also show that a statespecific (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 re-optimize the ground state molecular orbitals during the iterative SCRF procedure, such as cLR or VEM. Based on 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 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.

ACS Paragon Plus Environment

36

Page 37 of 51

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

The Journal of Physical Chemistry

ACKNOWLEDGMENT Acknowledgement is made to the donors of the American Chemical Society Petroleum Research Fund 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. acknowledges the University of California Merced Chancellor’s Postdoctoral Fellowship Function for financial support. We acknowledge computing time on the Multi-Environmental Computer for Exploration and Discovery (MERCED) cluster which is supported by National Science Foundation Grant No. ACI-1429783.

ASSOCIATED CONTENT Supporting Information. The Supporting Information is available free of charge on the ACS Publications website at DOI: XX. Ground and excited state electron density difference plots; Cartesian coordinates; basis set comparison; approximate IPCM excitation energies; ground-toexcited 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); comparison of SMD and IEF-PCM excitation energies; QM/classical excitation energies of pCT– using SMD radii; state-specific QM/PCM excitation energies; PCM cavity surface point distances.

REFERENCES

ACS Paragon Plus Environment

37

The Journal of Physical Chemistry

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

Page 38 of 51

(1)

Senn, H. M.; Thiel, W. QM/MM Methods for Biomolecular Systems. Angew. Chem. Int. Ed. Engl. 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 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.

ACS Paragon Plus Environment

38

Page 39 of 51

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

The Journal of Physical Chemistry

(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, 34109.

(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, 75102.

(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.; Leś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.; Marcos, E. S.; 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, 118–1123.

(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)

Martínez, J. M.; Pappalardo, R. R.; Marcos, E. S.; 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.

(26)

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.

(27)

Preat, J.; Loos, P.-F.; Assfeld, X.; Jacquemin, D.; Perpète, E. A. DFT and TD-DFT Investigation of IR and UV Spectra of Solvated Molecules: Comparison of Two SCRF

ACS Paragon Plus Environment

39

The Journal of Physical Chemistry

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

Page 40 of 51

Continuum Models. Int. J. Quantum Chem. 2007, 107, 574–585. (28)

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.

(29)

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.

(30)

Connolly, M. L. Analytical Molecular Surface Calculation. J. Appl. Crystallogr. 1983, 16, 548–558.

(31)

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.

(32)

Cossi, M.; Mennucci, B.; Cammi, R. Analytical First Derivatives of Molecular Surfaces with Respect to Nuclear Coordinates. J. Comput. Chem. 1996, 17, 57–73.

(33)

Connolly, M. Solvent-Accessible Surfaces of Proteins and Nucleic Acids. Science 1983, 221, 709–713.

(34)

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.

(35)

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.

(36)

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.

(37)

Cramer, C. J.; Truhlar, D. G. Implicit Solvation Models: Equilibria, Structure, Spectra, and Dynamics. Chem. Rev. 1999, 99, 2161–2200.

(38)

Fattebert, J.-L.; Gygi, F. Density Functional Theory for Efficient Ab Initio Molecular Dynamics Simulations in Solution. J. Comput. Chem. 2002, 23, 662–666.

(39)

Fattebert, J.-L.; Gygi, F. First-Principles Molecular Dynamics Simulations in a Continuum Solvent. Int. J. Quantum Chem. 2003, 93, 139–147.

(40)

Scherlis, D. A.; Fattebert, J.-L.; Gygi, F.; Cococcioni, M.; Marzari, N. A Unified Electrostatic and Cavitation Model for First-Principles Molecular Dynamics in Solution. J. Chem. Phys. 2006, 124, 74103.

(41)

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.

ACS Paragon Plus Environment

40

Page 41 of 51

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

The Journal of Physical Chemistry

Europhys. Lett. 2011, 95, 43001. (42)

Tomasi, J.; Mennucci, B.; Cancè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.

(43)

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.

(44)

Scalmani, G.; Frisch, M. J. Continuous Surface Charge Polarizable Continuum Models of Solvation. I. General Formalism. J. Chem. Phys. 2010, 132, 114110.

(45)

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.

(46)

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.

(47)

Klamt, A.; Schüü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.

(48)

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.

(49)

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.

(50)

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.

(51)

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.

(52)

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.

(53)

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 State-Specific Polarizability. J. Chem. Theory Comput. 2013, 9, 3649–3659.

(54)

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.

ACS Paragon Plus Environment

41

The Journal of Physical Chemistry

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

Page 42 of 51

(55)

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.

(56)

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.

(57)

Cammi, R.; Mennucci, B.; Tomasi, J. Fast Evaluation of Geometries and Properties of Excited Molecules in Solution: A Tamm-Dancoff Model with Application to 4Dimethylaminobenzonitrile. J. Phys. Chem. A 2000, 104, 5631–5637.

(58)

Cossi, M.; Barone, V. Time-Dependent Density Functional Theory for Molecules in Liquid Solutions. J. Chem. Phys. 2001, 115, 4708.

(59)

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.

(60)

Improta, R.; Barone, V.; Scalmani, G.; Frisch, M. J. A State-Specific Polarizable Continuum Model Time Dependent Density Functional Theory Method for Excited State Calculations in Solution. J. Chem. Phys. 2006, 125, 54103.

(61)

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, 74504.

(62)

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.

(63)

Corni, S.; Cammi, R.; Mennucci, B.; Tomasi, J. Electronic Excitation Energies of Molecules in Solution within Continuum Solvation Models: Investigating the Discrepancy between State-Specific and Linear-Response Methods. J. Chem. Phys. 2005, 123, 134512.

(64)

Runge, E.; Gross, E. K. U. Density-Functional Theory for Time-Dependent Systems. Phys. Rev. Lett. 1984, 52, 997–1000.

(65)

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.

(66)

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.

(67)

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.

ACS Paragon Plus Environment

42

Page 43 of 51

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

The Journal of Physical Chemistry

Theor. Chem. Acc. 2016, 135, 263. (68)

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.

(69)

Hu, L.; Söderhjelm, P.; Ryde, U. On the Convergence of QM/MM Energies. J. Chem. Theory Comput. 2011, 7, 761–777.

(70)

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.

(71)

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.

(72)

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.

(73)

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.

(74)

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.

(75)

Isborn, C. M.; Gö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.

(76)

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.

(77)

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.

(78)

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.

(79)

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

ACS Paragon Plus Environment

43

The Journal of Physical Chemistry

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

Page 44 of 51

Models. J. Phys. Chem. B 2016, 120, 12148–12159. (80)

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.

(81)

Martínez-Ferná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.

(82)

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.

(83)

Hirata, S.; Head-Gordon, M. Time-Dependent Density Functional Theory within the Tamm–Dancoff Approximation. Chem. Phys. Lett. 1999, 314, 291–299.

(84)

Chai, J.-D.; Head-Gordon, M. Systematic Optimization of Long-Range Corrected Hybrid Density Functionals. J. Chem. Phys. 2008, 128, 84106.

(85)

Nielsen, I. B.; Boyé-Pé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.

(86)

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.

(87)

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, 34108.

(88)

Zuev, D.; Bravaya, K. B.; Crawford, T. D.; Lindh, R.; Krylov, A. I. Electronic Structure of the Two Isomers of the Anionic Form of P-Coumaric Acid Chromophore. J. Chem. Phys. 2011, 134, 34310.

(89)

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.

(90)

Lipparini, F.; Lagardère, L.; Raynaud, C.; Stamm, B.; Cancè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.

(91)

Setny, P.; Dudek, A. Explicit Solvent Hydration Benchmark for Proteins with Application to the PBSA Method. J. Chem. Theory Comput. 2017, 13, 2762–2776.

ACS Paragon Plus Environment

44

Page 45 of 51

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

The Journal of Physical Chemistry

TOC Graphic

ACS Paragon Plus Environment

45

The Journal of Physical Chemistry

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

82x67mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 46 of 51

Page 47 of 51

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

The Journal of Physical Chemistry

82x155mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

82x137mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 48 of 51

Page 49 of 51

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

The Journal of Physical Chemistry

82x121mm (300 x 300 DPI)

ACS Paragon Plus Environment

The Journal of Physical Chemistry

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

82x60mm (300 x 300 DPI)

ACS Paragon Plus Environment

Page 50 of 51

Page 51 of 51

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

The Journal of Physical Chemistry

82x72mm (300 x 300 DPI)

ACS Paragon Plus Environment