Implicit Electrostatic Solvent Model with Continuous Dielectric

Jan 29, 2010 - of a variable ε function inherent to the solvent region. The obtained .... Consider now the case where ε(r) is a continuous function...
0 downloads 0 Views 166KB Size
Download PDF
J. Phys. Chem. B 2010, 114, 2457–2466

2457

Implicit Electrostatic Solvent Model with Continuous Dielectric Permittivity Function Mikhail V. Basilevsky,† Fedor V. Grigoriev,‡ Ekaterina A. Nikitina,*,§ and Jerzy Leszczynski| Photochemistry Center and Institute of Applied Mechanics, Russian Academy of Sciences, Moscow, Russia, Research Computing Center, Moscow State UniVersity, Moscow, Russia, Dimonta LLC, Moscow, Russia, and Interdisciplinary Nanotoxicity Center, Jackson State UniVersity, Jackson, Mississippi 39217 ReceiVed: July 8, 2009; ReVised Manuscript ReceiVed: October 22, 2009

The modification of the electrostatic continuum solvent model considered in the present work is based on the exact solution of the Poisson equation, which can be constructed provided that the dielectric permittivity ε of the total solute and solvent system is an isotropic and continuous spatial function. This assumption allows one to formulate a numerically efficient and universal computational scheme that covers the important case of a variable ε function inherent to the solvent region. The obtained type of solution is unavailable for conventional dielectric continuum models such as the Onsager and Kirkwood models for spherical cavities and the polarizable continuum model (PCM) for solute cavities of general shape, which imply that ε is discontinuous on the boundary confining the excluded volume cavity of the solute particle. Test computations based on the present algorithm are performed for water and several nonaqueous solvents. They illustrate specific features of this approach, called the “smooth boundary continuum model” (SBCM), as compared to the PCM procedure, and suggest primary tentative results of its parametrization for different solvents. The calculation for the case of a binary solvent mixture with variable ε in the solvent space region demonstrates the applicability of this approach to a novel application field covered by the SBCM. 1. Introduction Continuum solvent models are widely used in chemistry, biochemistry, biophysics, and related fields.1-5 In such models, the solvation effects are usually separated into electrostatic and nonelectrostatic components and combined as additive contributions coming from these components. The present work is devoted to studying the electrostatic continuum theory of solvation based on the solution of the Poisson equation. Currently, a number of well-developed computational methodologies utilize the treatment of electrostatic solvation effects, in practice. In the conventional formulation, the solvent is considered as a structureless medium with static dielectric permittivity ε ) ε0. It surrounds the cavity of excluded volume assembled around the solute particle; the condition ε ) 1 is assumed inside the cavity. According to this model, the dielectric permittivity changes in a stepwise manner on the surface representing the cavity boundary. This is the well-defined electrostatic problem consisting of the computation of the electric response field created by the continuum solvent, which is polarized by the electric charge distribution inherent to the solute. The popular PCM (polarizable continuum model) procedure1–3 is confined to the most common case of constant ε0, which allows one to reduce the three-dimensional (3D) Poisson equation in real space to the two-dimentional (2D, either iterative or direct linear) scheme of finding the polarization charge density on the cavity surface. The alternative methodology suggests the straightforward numerical 3D procedure for solving the Poisson equation with properly defined boundary conditions.4,5 A number * Corresponding author. Phone: +7 095 9350115. E-mail: [email protected]. † Photochemistry Center, Russian Academy of Sciences. ‡ Moscow State University and Dimonta LLC. § Institute of Applied Mechanics, Russian Academy of Sciences. | Jackson State University.

of more simplistic models also exist [for instance, the generalized Born approximation (GBA)3,6 and the conductor-like screening model (COSMO)7], suggesting various more or less rigorous approximations for easily solving the same basic electrostatic problem. The condition of the stepwise rupture of ε on the cavity surface is the common feature for all of these approaches. The original idea was first applied to the simplest model solute cases by Onsager8 and Kirkwood.9 The procedure of matching electric fields on the boundary surface is the most essential feature of this methodology. Owing to such a strategy, the solution of the Poisson problem cannot be avoided. On the other hand, there exists a special approach that seriously simplifies the desired solution. It can be appropriate for application provided that the dielectric displacement field D ) εE (where E represents the total electric field strength) exactly coincides with the vacuum electric field G0, created by the same solute charge in empty space (ε ) 1 everywhere). This idea was formulated as the D-G theorem,10 that is

D ) G0

(1)

the relation that provides the desired solution with no extra labor. It is therefore very important to establish the conditions under which eq 1 is exact. This issue was thoroughly analyzed by Kharkats, Kornyshev, and Vorotyntsev,10,11 the earlier literature on this point can be found in the references cited in these works. The basic conclusion was that, for the electrostatic problems formulated in terms of the stepwise ε functions, this relation is hardly typical; it remains exactly true only in exceptional cases (for instance, in the case of a Born spherical ion). On the other hand, eq 1 is, as a rule, obeyed, provided that ε is considered as a continuous and isotropic function. In the present work, we accept the idea of continuous and isotropic ε in order to

10.1021/jp9064399  2010 American Chemical Society Published on Web 01/29/2010

2458

J. Phys. Chem. B, Vol. 114, No. 7, 2010

formulate an exact computational algorithm for solving the Poisson equation for applications in the solution chemistry. Its important additional benefit is that ε can be variable within this approach. Few dielectric continuum computations of molecularshaped solutes with variable ε have been reported in the recent literature.12–14 The attempt using a smooth ε model12 was provided with no further elaboration. The important study based on variable ε smoothly changing across the planar interface13 retained the basic PCM postulate of an abrupt ε function on the solute boundary. In the present work, the idea of variable and purely continuous ε is developed and systematically elaborated. The new algorithm designed for treating a nonuniform continuum with variable ε is formulated in section 2. In section 3, it is supplemented by an explicit prescription for constructing the continuous ε function on the boundary of the solute cavity, including the case of arbitrarily shaped solutes. (The earlier approaches that considered variable ε models were confined to the few simplest solution systems: either spherical Born-like ions14 or two-component solvent layers with a smooth planar interface13). In the present work, the combination of the two suggested algorithms gives rise to a universal and internally consistent computational methodology. Its applications and efficiency are illustrated in sections 4 and 5. An additional remark is required regarding the interrelation between the present approach, henceforth called SBCM (smooth boundary continuum model), and the conventional and widespread model with stepwise ε, called PCM.1–3 In the SBCM scheme, the ε function depends on several input parameters whose changing deliberately could ultimately result in a stepwise ε. It would be natural to expect that the PCM and the corresponding limit of the SBCM are equivalent. This issue is not so simple, however: the identity of the two stepwise approaches can be gained only in a few exceptional cases when eq 1 is true. Usually, the two schemes (the PCM and the stepwise limit of the SBCM) provide results that are numerically different, although both remain physically relevant and applicable for treating real solution systems under proper parametrization. Nevertheless, the resulting unusual and confusing situation should be clarified. Its analysis requires careful mathematical elaboration, which is deferred to Appendix A and summarized in the concluding section 6. 2. Theoretical Background 2.1. Basic Energy Relation. Consider electrostatic fields created by a charge distribution F(r) immersed in a continuous medium with dielectric permittivity ε(r), where r denotes the space point. The free energy corresponding to the total electric potential field ψ(r) and the corresponding field strength (E ) -∇ψ) is

Gtot )

1 8π

1 ∫ -∇ψ(-ε∇ψ) d3r ) 8π ∫ d3r ε(∇ψ)2

where -ε∇ψ denotes the dielectric displacement D. The vacuum free energy equals

Gvac )

1 8π

∫ d3r (∇ψ0)2

where ψ0 is the vacuum potential of the charge F(r). The electrostatic solvation energy is then defined as

Basilevsky et al.

1 8π

Gel ) Gtot - Gvac )

∫ d3r [ε(∇ψ)2 - (∇ψ0)2] (2)

This expression is known in the literature.10,11 Consider now the case where ε(r) is a continuous function. In particular, if S denotes the boundary of the solute cavity of excluded volume (where ε ) 1 by the definition), then normal derivatives (∂ψ/∂n)i and (∂ψ/∂n)e, counted on S from its internal (subscript i) and external (subscript e) regions, are also continuous provided ε(r) is continuous, that is

εi(r) ) εe(r),

∂ψ )( ) ( ∂ψ ∂n ) ∂n i

e

(3)

This is a consequence of the standard matching condition formulated for the dielectric displacement field D, namely, εi(∂ψ/ ∂n)i ) εe(∂ψ/∂n)e. For this special case, the following basic relation is valid [i.e., the D-G theorem ( eq 1)10]: D ) -∇ψ0 or

1 ∇ψ ) ∇ψ0 ε

(4)

This is derived below. By inserting eq 4 into eq 2, we obtain

Gel ) -

1 8π

∫ d3r (1 - 1ε )(∇ψ0)2

(5)

Because ε(r) and ∇ψ0 are supposed to be known, eq 5 provides the solvation energy without solving the Poisson equation. To certify eq 4, we consider the pertaining Poisson equation

∇(ε∇ψ) ) -4πF

(6)

Under the conditions in eq 3, it is obeyed identically. This is seen when eq 4 is substituted in its left-hand part, because the conventional Poisson equation ∇ψ0 ) -4πF is satisfied for the vacuum field with ψ0 given by

ψ0(r) )

∫ d3r |rF(r') - r'|

(7)

Thereby, eq 4 is indeed the exact solution to 6, provided that ε(r) is continuous. The continuity condition in eq 3 is essential. Discontinuity of ε(r) (on the boundary of the surface S, for example) generates singularities in eq 6 that make an additional contribution to its solution. As mentioned in section 1, the typical case of the discontinuity of ε(r) on S almost always discards eqs 4 and 5, with the exclusions of this rule being exceptional. The detailed consideration of this issue is deferred to Appendix A. 2.2. Total Electrostatic Potential and Response Potential Field. We accept now the new condition that the charge distribution F(r) is completely confined in the solute cavity [i.e., F(r) ) 0 when r is outside S]. This condition can be formulated as

(ε - 1)F ) 0

(8)

Implicit Electrostatic Solvent Model

J. Phys. Chem. B, Vol. 114, No. 7, 2010 2459

Equation 8 allows one to derive the explicit expression for ε(r) starting from eq 4. We introduce the medium response field Φ(r) ) ψ(r) - ψ0(r). This relation is substituted into Poisson eq 6 to give

1 ∇ε · ∇ψ0 - 4π(ε - 1)F ) -ε∇2Φ ε Equations 4 and 6 were used here. If we invoke now the boundary condition in eq 8, the second term on the left-hand side vanishes. Thereby, we obtain the final result

∇2Φ ) -4πg 1 g) (∇ε · ∇ψ0) 4πε2

∫ d3r' |rg(r') - r'|

(9)

(10)

1 2

∫ d3r F(r) Φ(r) ) 21 ∫ ∫ d3r d3r' F(r) |rg(r') - r'| (12)

This expression has the transparent interpretation as the interaction energy of solute and solvent charges. We can also define the energy corresponding to the total charge in eq 10 as

Gtot )

1 2

+ g(r')] ∫ ∫ d3r d3r' [F(r) + g(r)][F(r') |r - r'|

and the self-energies of charges F and g as

GF )

1 2

F(r')] , ∫ ∫ d3r d3r' [F(r) |r - r'| Gg )

1 2

g(r')] ∫ ∫ d3r d3r' [g(r) |r - r'|

In this way, we find that

Gel ) Gtot - GF - Gg

(14)

where qi and ri denote the charge and location, respectively, of site i of the solute. Function ε(r) must be formulated so as to capture its steep variation on the solute boundary surface S from ε ) 1 inside the cavity to the bulk static value ε ) 〈ε〉 ) ε0. This represents the simplest case when no other changes of ε are introduced; other models can also be considered within the present algorithm. The desired formulation of the space dependence of the dielectric permittivity will be made in terms of variable number density n(r) (the number of solvent particles per unit volume). We introduce the bulk density 〈n〉 and the dimensionless number density z ) n/〈n〉. The appropriate expression arises based on dielectric susceptibilities χ∞ and χin, namely, the inertialess electronic and inertial nuclear components, respectively, of the total dielectric susceptibility of the solvent15

(11)

The solvation energy in eq 2 is then obtained in the alternative and more usual form

Gel )

q

∑ |r -i ri| i

Finally, Φ(r) is expressed in terms of g(r) as

Φ(r) )

3. Parameterization of the Permittivity Function 3.1. Dielectric Permittivity as a Function of the Number Density. Equations 5, 9, and 11 provide the new explicit algorithm for performing a complete solvation calculation. Quantities ψ0(r) and ε(r) are required only as input data. For the vacuum potential, the obvious site representation can be suggested

ψ0(r) )

This is a formal Poisson equation for Φ(r). Here, g(r) serves as a charge distribution, which can be interpreted as the polarization charge induced in the medium by the solute charge F(r). Note that g(r) ) 0 inside the cavity because ∇ε ) 0 there. Therefore, distributions F(r) and g(r) do not overlap. We can also define the total charge distribution, due to both the solute and the solvent, as

Ftot ) F + g

fields or, alternatively, in terms of electrostatic charges corresponding to these fields.

(13)

So eqs 2 and 5 on one hand and eqs 12 and 13 on the other represent the solvation energy Gel either in terms of electrostatic

ε(r) ) 1 + 4π{χ∞z(r) + χinz(r) φ[z(r)]}

(15)

Here

φ(z) ) z

1 + z0 z + z0

[

φ(z) ) exp

(z - 1) z0

(16a)

]

(16b)

The inertialess susceptibility in eq 15 corresponds to the fast electronic component of solvent polarization variable, whereas the inertial one represents the sluggish, slowly relaxing polarization component; it is mainly related to nuclear (orientational and translational) solvent motions. With z ) 1, quantities χ∞ and χin represent the uniform bulk solvent model. For this case, 1 + 4πχ∞ ) ε∞ is the optical dielectric permittivity, and the total permittivity of eq 15 is the static quantity ε0 ) 1 + 4π(χ∞ + χin). According to eq 15, the inertialess component of ε is proportional to the density (which is a conventional assumption16,17), whereas the inertial component displays a stronger dependence on z;16–18 its actual trend is controlled by the parameter z0.15 Two versions for scaling function φ(z) are suggested in eqs 16a and 16b; they are normalized such that φ(z) ) 1 when z ) 1. 3.2. Space Dependence of the Number Density. We suggest the following model for introducing the continuous density change occurring on the cavity boundary. It is based on the PCM representation of a cavity surface S as a collection of overlapping atomic-site spheres. Figure 1 shows the piece of S containing two neighboring spheres i and j with radii Ri and Rj,

2460

J. Phys. Chem. B, Vol. 114, No. 7, 2010

Basilevsky et al. TABLE 1: Electrostatic Components of the Hydration Free Energies (kcal/mol) for Ammonium-Type Ions cation +

Figure 1. Construction of the density function z(r), eq 17. Curve AB divides the regions of spheres i and j.

Figure 2. Morse-like interpolation of the solvent density function z(r) on the solute boundary. See Figure 1 for the definition of li.

respectively. For every space point r, we draw vectors b L k, connecting this point with the centers of spheres k (k ) 1, 2, bk| and the ...), and obtain the corresponding distances Lk ) |L separation lengths lk ) Lk - Rk (i.e., the distance between r and the point closest to it on each sphere k). The following procedure is then performed (see Figure 1): We consider all distances lk (k ) 1, 2, 3, ...) and find the one (k ) i, Figure 1) that is the smallest: li < lk (k * i). Then, we acknowledge that the point r belongs to the region of sphere i, that is, r ∈{ri}. With this specification, the prescription for calculating z within the region of sphere i reads as

z(r) ) zi(r) ) 1 + exp(-2Rli) - 2 exp(-Rli) li > 0, r ∈ {ri} li < 0, r ∈ {ri} 0 (17)

{

Equation 17 represents a fragment of a Morse-like function (Figure 2). It vanishes on the sphere surface (li ) 0), together with its derivative, and asymptotically reaches the value z ) 1, when Rli . 1. Parameter R regulates the steepness of this z function. Finally, consider the boundary surface dividing the regions of the neighboring spheres i and j. It is schematically drawn as curve AB in Figure 1 and counts the points equidistantly remote from the surfaces of spheres i and j (i.e., li ) lj at points belonging to this dividing surface). One can infer from eqs 15 and 17 that ε(r) is continuous on this boundary surface. This property is assured by the present algorithm automatically, and no special effort is required to determine this surface explicitly. By this means, the suggested procedure yields, indeed, such a function ε(r) that has no rupture either on the solute boundary or in the bulk solvent region. Its continuity is assured both on the sphere surfaces (owing to function 17) and on the surfaces dividing neighboring regions {ri} and {rj}. This is the necessary property underlying the methodology elaborated in section 2.19 The normal derivatives are continuous on the boundary but are not matched on the dividing surface AB. The latter fact creates no essential problems. 4. Illustrative Computations 4.1. Hydration Free Energies of Organic Ions. Test computations were carried out for a series of ammonium-like

CH3NH3 CH3CH2NH3+ CH3(CH2)2NH3+ (CH3)2CHNH3+ CH3(CH2)3NH3+ (CH3)3CNH3+ (CH3)2NH2+ (CH3CH2)2NH2+ (CH3)3NH+ pyrrolidinium R-methylpyrrolidinium pyridinium piperidinium anilinium imidazolium K+ e Cl- e I- e Na+ e

hydrocarbona

∆Gs_Expb

PCMc

SBCM-2d

CH3CH3 CH3CH2CH3 CH3(CH2)2CH3 (CH3)2CHCH3 CH3(CH2)3CH3 (CH3)3CCH3 (CH3)2CH2 (CH3CH2)2CH2 (CH3)3CH cyclopentane methylcyclopentane benzene cyclohexane toluene cyclopentadiene Ar Ar Xe Ne

-73.1 -70.4 -68.8 -67.8 -68.5 -65.6 -65.9 -61.2 -58.9 -62.8 -56.2 -55.2 -61.2 -65.11 -59.7 -83.1 -79.1 -64.5 -100.8

-72.4 -69.6 -68.8 -67.2 -68.3 -65.0 -64.4 -59.7 -57.6 -61.0 -59.0 -55.7 -59.4 -67.7 -60.3

-72.9 -71.8 -67.0 -65.8 -68.3 -62.3 -68.0 -58.1 -58.1 -62.5 -58.1 -57.1 -58.9 -66.6 -61.81 -81.75 -793 -63.3 -100.1

a Hydrocarbon structure simulating the pertaining cation. Electrostatic part of experimental free energies obtained by eliminating their nonelectrostatic contributions as explained in the text. The background experimental data, collected in refs 20 and 21, are based on original measurements cited in these references. c κ ) 1.06 (see eq 18). d See parameter values in section 4.1. e Experimental estimates based on refs 21 and 24. b

cations for which experimental solvation free energies in water are available. The electrostatic fraction of this energy was extracted by subtracting the solvation free energies of the related hydrocarbons with the same structures (as illustrated by pairs NH4+/CH4, pyridinium/benzene, etc.). Experimental data derived by this procedure are collected in Table 1. The basic experimental material is published in refs 20 and 21. The free energy values cited in these articles are supplementary; when they overlap, the discrepancies do not exceed 1-2 kcal/mol. The electrostatic computations were based on atomic-site partial charges obtained by quantum-chemical Hartree-Fock computations [6-31G(d,p) basis]22,23 according to the electrostatic potential (ESP) charge localization scheme.22 For later convenience the present computational procedure is abbreviated as SBCM (smooth boundary continuum model). We consider two versions: SBCM-1 (based on eq 5) and SBCM-2 (based on eqs 9, 11, and 12). As discussed in section 5, algorithm SBCM-2 is more complicated and becomes computationally more expensive at the stage of preparation of the integrand entering the final energy expression. However, the energy integral converges significantly more rapidly with SBCM-2. Therefore, we preferred this method in computations of ions, because the longrange tails of the integrands significantly hampered the convergence of the integral in eq 5. The SBCM-2 results are compared with conventional PCM computations in Table 1. Both procedures were applied within the simplified scheme neglecting the polarization of solutes by the solvent response field, that is, beyond the SCRF (self-consistent reaction field) nonlinear methodology. The electrostatic atomic radii, Rel, were calibrated by scaling the van der Waals radii, RvdW, as

Rel ) κRvdW

(18)

The studies in this and the following sections are entirely based on the list of RvdW values given in Table 5 of ref 21. The scaling factor κ ) 1.06 was found in our PCM computations. A similar quality of fitting the experimental free energies with the present

Implicit Electrostatic Solvent Model

J. Phys. Chem. B, Vol. 114, No. 7, 2010 2461 components. The SBCM-2 computations for the same electrostatic contributions were also performed. The two arrays of data are compared in Table 2. Fitting of the computational data to their quasi-experimental counterparts was gained by the variation of the SBCM parameters κ, R, and z0, yielding the following values (with van der Waals radii from ref 21):

κ ) 0.9, Figure 3. SBCM electrostatic interaction profiles in a model solvent with variable R and ε0 ) 78 for ion pair K+ Cl- and their convergence to the PCM.

SBCM method was accomplished by applying κ ) 0.92, R ) 3.1 Å-1, and z0 ) 0.5 [eq 16b was applied for φ(z)]. Simple monatomic cations and anions are also included in Table 1. They are involved in the forthcoming presentation. 4.2. Parameters of the SBCM Scheme. The meaning of parameter κ in eq 18 is transparent. It scales electrostatic radii controlling the solute shape in both the PCM and SBCM approaches. The new parameters R (eq 17) and z0 (eqs 16a and 16b) are specific for the SBCM. They determine the steepness of the dielectric permittivity profile on the boundary of the solute cavities. According to eq 17, R serves for tuning the width of the smoothed-over stepwise density change. It is expected that the large R limit shifts the SBCM methodology asymptotically toward the PCM. This conjecture is illustrated by Figure 3, where the electrostatic R-dependent curves, representing the solvation contribution to the ion-pair interaction in K+ Cl-, are displayed as functions of the interionic distance. For consistency of comparison, the parameter κ ) 1.06 was assumed to be common for both the PCM and SBCM computations; we employed z0 ) 1.0 for the SBCM, whereas R was varied in a wide range. It should be remembered that the complete coincidence of the PCM and SBCM results is impossible even in the limit of infinite R, because of the different matching conditions on the boundary surface (see Appendix A). According to eqs 15, 16a, and 16b, parameter z0 discriminates the density dependences inherent to inertial and inertialess solvent polarizations. Combined with R, this provides finer tuning of the electrostatic behavior of solution systems when they change on the cavity boundary. We expect also that 1/R must correlate with the characteristic size of the solvent particles. This idea is supported by computations of free solvation energies for several nonwater solvents. The experimental data required for such calibration of parameters is notably worse than that available in hydration systems. Systematic measurements of solvation free energies for organic ions in nonaqueous media are usually lacking. Therefore, we considered the experimental data for electrically neutral organic compounds, whose absolute values are almost by an order of magnitude smaller than those representing solvation energies of ions in Table 1. The relative reliability of these data is certainly less convincing than that observed in the ionic case. Nonetheless, general trends found in these studies do not look unreasonable. Therefore, we examined series of organic nitrogencontaining compounds related to the ions included in Table 1. The pertinent information was extracted from the extensive collection of data in ref 6. An additional problem arises at the stage of separation of the electrostatic effect from the total experimental data. For this purpose, we applied the nonelectrostatic contributions calculated in the same ref 6 in terms of the SM5 version25 of the SASA (solvent-accessible surface area) scheme. They were subtracted from the true (total) experimental solvation energies to give the “experimental” electrostatic

R ) 1.5 Å-1

(benzene, ε0 ) ε∞ ) 2.2)

κ ) 0.5,

R ) 0.5 Å-1,

z0 ) 2.5 (chloroform, ε0 ) 4.7, ε∞ ) 2.2)

κ ) 0.9,

R ) 0.8 Å-1,

z0 ) 3.0 (isobutanol, ε0 ) 17.7,ε∞ ) 2.2)

In the data for benzene, introduced following eqs 15-17, we assumed χin ) 0, which implies that the parameter z0 should be withdrawn. For simplicity we also assumed ε∞ ) 2.2 for all three solvents, the experimental values being 2.25 (benzene), 2.1 (chloroform), 2.0 (isobutanol). Invoking the SM5 results6 was the most easily available way to estimate tentatively the quasi-experimental electrostatic free energy Gel for nonaqueous solvents. Otherwise, we would be involved in the problem of a consistency with the SBCM parametrization of the nonelectrostatic free energy component, which is beyond the scope of the present work. Note that the alternative simple route that we applied for water solutions (section 4.1 and Table 1) is now unavailable because of the lack of the experimental data for solvation energies of the hydrocarbon systems similar to the compounds considered in Table 2. It is clearly seen that decreasing R, as compared to the water case, correlated with the fact that the molecules of the nonwater solvents considered are significantly bulkier than water particles. For nonwater solvents, the so-determined parameter values provide preliminary tentative estimates to be refined in a more extended investigation based on more diverse and representative experimental background. (For instance, the value κ ) 0.5, found above for chloroform, seems doubtful.) It is also worth commenting that data for aromatic solutes were excluded from Table 2. Their presence strongly decreases the quality of fitting, which can be attributed to the solute polarization effect associated with aromatic π-electron systems. It is neglected according to the present linearized SBCM scheme but can notably contribute to the solvation energies. On the other hand, this effect is less visible in Table 1, where several aromatic cations are included. Seemingly, this is so because the cationic polarizabilities are reduced, as compared to those of neutral molecules. We have also excluded from Table 1 the cation NH4+, for which the solvation free energy proved to be much too high (85 kcal/mol as found with the present parametrization versus 81 kcal/mol according to the experimental estimate). This discrepancy could be associated with the usually encountered variety of problems appearing for the continuum solvent models when the solvation of small ions is considered19 (see also the discussion in refs 15, 26, and 27 and references therein). 4.3. Potential of Mean Force for an Ion Pair in a Binary Solvent Mixture. We considered finally the case of nontrivial space-dependent dielectric permittivity. This effect is specific for a binary mixture of polar and nonpolar solvents. When the solute in such a system is polar or charged, the polar solvent ingredient accumulates in its vicinity. This means that its

2462

J. Phys. Chem. B, Vol. 114, No. 7, 2010

Basilevsky et al.

TABLE 2: SBCM-2 Computations of the Electrostatic Free Energy Component (kcal/mol) for Nonaqueous Solutions benzene molecule

experiment

methylamine ethylamine dimethylamine propylamine trimethylamine butylamine diethylamine piperidine dipropylamine ammonia

-1.8 -1.0 -1.3 -1.1 -0.5 -1.0 -0.7 -0.3 -0.2 -1.1

a

a

chloroform

computation

b

experiment

-1.2 -1.1 -0.8 -1.1 -0.4 -1.0 -0.7 -0.7 -0.6 -1.4

a

isobutanol

computation

-2.1 -2.1 -1.6 -2.1 -1.1 -1.9 -1.7 -1.4

-2.1 -2.1 -1.5 -2.0 -0.8 -2.2 -1.5 -1.3

-2.4

-2.3

b

experiment

a

computationb

-3.0

-2.4

-2.0

-1.7

-1.2 -2.0 -1.3 -1.2 -1.0

-0.8 -2.3 -1.4 -1.5 -1.5

b

Reference 6; see text (section 4.2). See parameter values in section 4.2.

concentration, as well as the overall composition of the system, becomes essentially space-dependent. Provided that such situation (called “dielectric enrichment”) is treated within a continuum solvent model, the dielectric permittivity of the solvent must be considered as a space-dependent function.15 We denote as c1 and c2 local concentrations of nonpolar and polar solvent ingredients, respectively, with c ) c1 + c2 being the total local concentration. The average concentrations are denoted as 〈c1〉, 〈c2〉, and 〈c〉 ) 〈c1〉 + 〈c2〉, respectively. The most essential dimensionless quantities are y ) c2/c (the local fraction of the polar ingredient), x ) 〈c2〉/〈c〉 (the input quantity, a measure of the average solvent composition) and z ) c/〈c〉 (which is identified with the dimensionless number density in eqs 15 and 17). The extension of dielectric permittivity expression (eq 15) for a binary solvent reads

ε(r) ) 1 + 4π{χ∞z(r) + χinz(r) y(r) φ[z(r)]}

(19)

We assumed here that the optical dielectric susceptibility χ∞ represents the same constant for both polar and nonpolar components. The local equilibrium condition suggests the following equation for the function y(r)15

ln

y x F(r) φ(z) ) ln + 1-y 1-x kBT〈c〉

(20)

where F(r) ) (χin/2ε2)(∇ψ0)2, ε ) ε(r), y, and z are r-dependent as well; φ is defined in eqs 16a and 16b; and ∇ψ0 is the gradient of vacuum field ψ0(r) (eq 14) for the ion pair. After solving simultaneously eqs 19 and 20, with z(r) given by eq 17, both functions ε(r) and y(r) become available. This strategy was applied earlier15 to treat the simplest solute case of a spherical Born-like ion. It is numerically realized now for a more sophisticated problem, a pair of spherical ions. We calculated the potential of mean force (PMF) for the K+ Cl- pair in a benzene (nonpolar)/dimethylsulfoxide (DMSO, polar) mixture. The basic permittivities are ε∞ ) 2.2 for the two solvent ingredients, ε0 ) 46.5 for DMSO, and ε0 ) ε∞ for benzene; other details can be found in ref 15. Figure 4 represents the PMF U(d) as a function of interionic distance d for several sample cases of solvent composition x (i.e., the average DMSO fraction). The following approximation was implemented

U(d) ) -

Q2 + UvdW + Uel d

(21)

Figure 4. Potentials of mean force for ion pair K+ Cl- in the binary solvent mixture benzene/DMSO for different compositions x ) 〈c2〉/ 〈c〉 represented as the percentage of DMSO.

where Q ) (1 represents the ion charges, UvdW is the conventional van der Waals (Lennard-Jones) pair potential with LJ parameters21 σ ) 4.935 and 4.417 Å for K+ and Cl-, respectively, and ε ) 0.00033 and 0.118 kcal/mol for K+ and Cl-, respectively. The electrostatic component of solvation energy, Uel, was calculated by eq 5, where ε was found by solving eqs 19 and 20, with parameters κ ) 0.9, R ) 1.5 Å-1, and z0 ) 2 being applied [with the exponential model in eq 16b for φ(z)]. Figure 4 displays the typical case of a contact ion pair (for x ) 0, i.e., the pure benzene solvent) and also a sequence of steadily transformed solvent-separated pairs representing the gradual increase of the DMSO fraction. Note that the important free energy component is missing in eq 21. This is the solvent reorganization energy, promoted by the excluded-volume effect, which is distance-dependent. If taken into account, it would modify the depth of the PMF well and its shape at short separations. It could be evaluated in a molecular-level simulation as an additive correction to eq 21 under the condition that the electrostatic part of the solute/solvent interaction is eliminated. We have disregarded this nonelectrostatic “cavitation” correction in the present computation, essentially devoted to the electrostatic solvation effects. Figure 4 serves here to the objective of illustrating the application of the new dielectric continuum technique developed in the first part of the present work. A complete and detailed investigation and discussion of the behavior of ion pairs in a binary solvent mixture will be addressed in a separate publication. 5. Computational Experiments According to computations reported in section 4, the SBCM proved to be a workable method naturally adapted for treating systems in which the dielectric permittivity ε is an essentially variable quantity (such as that considered in section 4.3). For standard solvation problems, which can be alternatively treated in terms of the conventional PCM, the only principal advantage

Implicit Electrostatic Solvent Model of the SBCM would be the increased flexibility, gained owing to the increase in the number of parameters (which has an obvious physical interpretation; see section 4.2). From a purely computational viewpoint, the efficiency of the SBCM is determined by two factors. Being based on the exact solution of the Poisson equation (eq 4), the SBCM requires no extra effort for solving it numerically, which is an obligatory step both in the PCM version1-3 and within the PoissonBoltzmann (PB) algorithm.4,5 On the other hand, the threedimensional (3D) integration, needed at the stage of the free energy evaluation (eq 5) (SBCM-1), is expected to be numerically disadvantageous compared to the PCM, which is a 2D procedure. The situation is improved, however, provided that the alternative algorithm in eqs 9-12 (SBCM-2) is implemented for a free energy computation. Here, the 3D integration scheme is reduced to treat a relatively thin solvent layer around the solute cavity, in which ∇ε reaches values of significant magnitude (the thickness of this layer being on the order of 1/R). However, a technical complication still exists. It is generated by large gradients of ε, arising within this layer as well as on the boundary surfaces that divide the solvent regions belonging to neighboring atomic spheres of the solute (see section 3 and Figure 1). Treatment of the pertaining solvent regions requires finer numerical grids, and this problem becomes noticeable for R> 3 because higher accuracy in calculations of ∇ε is then required. The variable-ε case (in the bulk solvent region), which cannot be treated by the PCM, is, in principle, available in terms of the PB technique.4,5 The necessary amounts of 3D computations would be approximately the same in both the SBCM and the PB approach, with a reservation that the job of solving the 3D Poisson equation (in the PB approach) is substituted by a less laborious 3D integration (in the SBCM). Our practical experience showed approximately the same duration of PCM and SBCM-1 computations for systems included in Table 1. For the largest system with the solvent chloroform in Table 2, the acceleration of SBCM-1 relative to PCM is a factor of 3. Such a significant efficiency difference for different solvents arises because the computations in Table 1 (ionic solutes) and those in Table 2 (neutral solutes) require different 3D integration limits in the SBCM-1: the integration boxes used for cations were significantly larger than those used for neutral solutes. Both algorithms SBSM-1 and 2 were applied in these calculations, but only SBSM-2 results are listed in Tables 1 and 2. Finally, we considered a large solute, taking structure A (3,16androstan) as a sample case.

A cationic form of this compound in polar solvents in its ground electronic state has a pair of minima A1 and A2, in which the positive charge is located on the biphenyl and naphthalene fragments, respectively. The charge transfer is represented by the A1 f A2 process.28 The following computations for structure A2 were performed: (1) standard PCM computations in water (gas-phase geometry optimization and subsequent computation of partial atomic charges by SCRF/PCM22,29) and (2) one PCM and two SBCM (SBCM-1 and SBCM-2) computations for the given geometry and charge distribution

J. Phys. Chem. B, Vol. 114, No. 7, 2010 2463 TABLE 3: Efficiency of SBCM Computations (with Structure A2 in Water as a Test System) method

integration box edge (nm)

τ/τ(PCM)a

Gel/Gel(PCM)b

SBCM-1

10 15 20 25 30 4 6 10

0.13 0.30 0.55 0.92 1.41 0.21 0.37 0.88

0.927 0.944 0.952 0.962 0.965 0.986 0.986 0.986

SBCM-2

b

a Duration of the computation relative to that for PCM. Solvation free energy relative to that for PCM.

found in the first step (solvent water) with the same parameters as used in Table 1. The SBCM-1 run corresponded to algorithm 5 for the free energy. The SBCM-2 run consisted of eqs 9, 11, and 12 with the gradient of ε. The ratios of computational times and electrostatic solvation free energies are listed in Table 3. The PCM hydration free energy is -43.6 kcal/mol. Inspection of Table 3 demonstrates the advantage of the SBCM-2 algorithm. For boxes of the same size, it works more slowly than SBCM1; however, it rapidly converges with much smaller boxes. Therefore, the convergence of SBCM-1 was not fully completed even with a box edge of 30 nm, whereas an edge of 4 nm was quite sufficient for the SBCM-2 algorithm. 6. Concluding Remarks Combining the estimates of numerical efficiency given in section 5 with the results of the test SBCM applications in section 4, we come out with a rather optimistic conclusion on the perspectives of this methodology. For water solutions, the results, their quality, and the computational efficiency of the algorithm proved to be more or less similar between the PCM and SBCM techniques. For nonaqueous solvents, the present SBCM results look satisfactory, and from general reasoning, we expect them to be preferable to the PCM results. Although the present realization of the SBCM does not incorporate the SCRF procedure (required to treat the effect of the solute polarization), this refinement can be easily made in terms of the SBCM-2 algorithm. Meanwhile, the ability to effortlessly treat problems with variable ε values remains its essential advantage. The methodological background of the present computational scheme is less transparent. At the fundamental level, the interrelation between the SBCM and PCM schemes presents a rather sophisticated issue. For its understanding, a careful examination of the stepwise limit transition of the variable ε function in the SBCM framework is required. This topic is specifically addressed in Appendix A. We briefly summarize here the essential outcomes of the analysis. In the limit of an abrupt stepwise change in ε, the SBCM and PCM techniques produce values that can be significantly different. The two schemes, which formally represent the exact solution for the same Poisson equation, should be considered as two alternative models of the cavity with sharp boundary walls. A priori, the two models are equally legitimate as purely mathematical constructions. The discrepancy arises because of the different ways of performing the stepwise limit transition. The main SBCM problem in performing this limit is the ability to correctly reproduce the continuity of the total electric field potential ψ and its tangent derivative, which constitutes the fundamental condition of electrostatic theory (namely, curl E

2464

J. Phys. Chem. B, Vol. 114, No. 7, 2010

Basilevsky et al.

) 0 with E ) -∇ψ). Within the SBCM approach, the ε function always remains continuous provided that the parameters controlling its steepness remain finite. The resulting electric fields, remaining continuous as well, are forced to adjust to the discrete (or almost discrete) boundary change of ε, which becomes essentially stepwise under the PCM formulation. This is why the singular behavior of the SBCM fields in the close vicinity of the cavity boundary contrasts markedly with what appears under smooth boundary conditions, alternatively postulated in the PCM scheme. Related numerical problems are expected to arise ultimately for the SBCM computations in the deep stepwise extreme, owing to the large gradient of ε. In real practice, this situation is associated with large R values introduced in eq 17. As demonstrated in Appendix A, the case of two closely positioned opposite charges within the body of a solute particle proves to be especially difficult. On the other hand, in actual computations we encountered no numerical problems in the treatment of ion pairs with R increasing up to 30 Å-1. Inspection of Figure 3 shows that, with increasing R, the SBCM energy profiles approach their PCM counterparts quite closely at large and moderate interionic separations. However, the discrepancy indeed becomes clearly visible at shorter distances. It is important to note that, for R < 3 Å-1 (i.e., for typical values that appear when commonly used solvents are addressed), no reflective difficulties, such as those discussed above, are expected to exist. The caution is only required when R is very large or the cavity shape is unusually tangled. For such exceptionally complicated cases, the nearly singular behavior of the SBCM solutions deserves further, more careful investigation. Acknowledgment. M.V.B., F.V.G., and E.A.N. are grateful to the Russian Fund for Fundamental Research (RFBR Grant 07-03-00479-a) for the financial support of this work. The support from the State Contract of Russian Federation No. 02.523.11.3014-DM/08 is also acknowledged. Participation of Dr. M. Devashis at the initial stage of the study is greatly appreciated. Appendix A: Formal Justification of the SBCM A1. D-G0 Theorem j the SBCM solution for the electric field We denote as ψ potential (instead of ψ as employed in the main text); the corrected exact solution (provided such a correction is required) is now denoted as ψ. The two Poisson equations read

∇(ε∇ψ) ) ∇2ψ0 ) -4πF

(A1)

where ψ0 is the vacuum potential. The basic assumption in eq 4 now appears in the form

1 ∇ψ ¯ ) ∇ψ0 ε

(A2)

The response field Φ is defined as

ψ ) ψ0 + Φ

(A3)

j , the field Φ is determined by relations 9 Provided that ψ ) ψ and 11. The possible correction to the SBCM is introduced as

ψ)ψ ¯ + ξ1 1 ∇ψ ) (∇ψ0 + ∇ξ) ε

(A4)

where the connection between ξ1 and ξ is given by

1 ∇ξ1 ) ∇ξ ε

(A5)

With such notation, ξ obeys the Laplace equation

∇2ξ ) 0

(A6)

The solutions for eq A6 (with the proper asymptotic behavior) are singular functions, called the single-layer potential and double-layer potential (SLP and DLP), respectively.30 They are discontinuous, displaying, on a given boundary surface S, ruptures of the normal derivative ∂/∂n (the case of the SLP) or of the function itself together with its tangent derivatives ∂/∂t (the case of the DLP). Here, we consider as S the solute cavity boundary and also the cuts that divide the space regions belonging to the different spherical fragments of this cavity (see Figure 1). It is easy to demonstrate now that both ξ and ξ1 vanish as long as ε is continuous. The thus-formulated conjecture follows directly from eqs 9 and 11, as applied to the response field Φ j . Indeed, both Φ arising in the zeroth-order case, when ψ ) ψ and ∇Φ prove to be continuous provided that ε is continuous, j and ∇ψ j (we always imply that and the same is true for ψ ψ0 and ∇ψ0 are continuous). The addition of discontinuous corrections, forming ξ and ξ1, would bring about ruptures in ψ or ∇ψ, which are prohibited. We denote the continuity condition as ε+ ) ε-, with the subscripts + and - labeling the values of any function on the outer (+) and inner (-) sides of the surface j if ε+ ) ε-, and this constitutes the present S. Thereby, ψ ) ψ version of the D-G0 theorem,10,31 which underlies the SBCM. Singular corrections are necessary only when the original solution itself creates singularities. The deliberately introduced singularities are invoked in order to cancel the singularities inherent in the original approximation, and thus they provide a smooth total solution. A2. Step Limit Transition Consider the step limit of the dielectric permittivity ε(r) on the boundary surface S. Two different algorithms to perform such a transition exist. The first one (the SBCM algorithm) defines the process as a limit of a continuously changing ε(r) function that appears during the course of the specially prescribed continuous change of its parameters. The parameter values always remain finite, but they change smoothly in the desired direction. The limit R f ∞ in eq 17 provides a typical example. Within such a prescription, the continuity of ε (ε+ ) ε-) is maintained. In the case of eq 17 for any finite R value, the continuity of ε(l) is formally assured in the range 0 < l < 1/R, although the slope of this function becomes extremely steep when R is large. The second method (the PCM algorithm) introduces ε+ * ε- as the original boundary condition. Singularities obviously appear in the solution, and they have to be canceled by invoking the extra singular corrections as explained in section A1. It is shown in section A3 how this gives rise to the standard PCM result.

Implicit Electrostatic Solvent Model

J. Phys. Chem. B, Vol. 114, No. 7, 2010 2465

The results of the two limit algorithms can be the same only in the very special case that10

∂ψ0 ∂ξ 4π ∂p + )∂t ε0 - 1 ∂t ∂t

∂ψ0 )0 ∂t

This is the integral equation defining the dipole density p(r′), r′ ∈ S. Now we return to the Onsager dipole considered above. Inside the cavity, ε ) 1, and then the combination of eqs A4 and A5 provides the simple result

(A7)

Otherwise, the SBCM solution, remaining formally continuous, cannot satisfy one of the boundary conditions (namely, the continuity of ∂ψ/∂t), and asymptotically, it does not exactly coincide with its PCM counterpart. A3. Illustrations A transparent example is given by a point dipole placed at the center of a spherical cavity. This is the simplest version of the construction introduced in Figure 1. We use the new notation for Ri ) a (the cavity radius) and for r ) (R, θ, φ) (the spherical coordinates). The dipole moment is denoted as m. The notation here partly overlaps with that appearing in the main text, but this produces no confusion below. In the present case, ε ) 1 (R < a) and ε ) ε0 (R> a), whereas ψ0 ) (m cos θ)/R2. The scenario suggested by the SBCM algorithm assumes R f ∞ in eq 17 as described in section A2. The response field Φ (eq 11) inside the cavity, evaluated in this limit, is

Φ)-

( )

(R < a)

(A8.1)

)

(R < a)

(A8.2)

2 mR cos θ ε0 - 1 3 ε0 a3

The standard (PCM) result is8

Φ)-

(

mR cos θ ε0 - 1 ε0 + 1/2 a3

1 ∂ ∫S ds' p(r') ∂n' |r - r'|

∂ ∂p (ξ - ξ-) ) 4π , ∂t + ∂t

(A11)

The straightforward solution of eq A10 gives the standard expression (eq A8.2) for the response field Φ ) ξ (R < a). Given p, the external field ξ for R> a is extracted from A9. It is added j ) ψ0/ε0, according to eqs A4 as a correction ξ1 ) ξ/ε0 to ψ and A5. A similar calculation performed for the normal derivative in order to ensure the continuity of ε(∂ξ/∂n) creates another component of ξ, the single-layer potential (SLP). This correction proves to be absent in the present case. The corresponding SLP is implicitly included in j , expressed through eq A8.1. It happens the original SBCM solution ψ to be so because of the condition ∂ε/∂t ) 0, which is always true for the constant ε ) ε0 outside the cavity. Therefore, only the normal contribution survives in the scalar product (∇ε∇ψ0) representing the polarization charge of eq 9. As shown in the next section, it generates j. the singular piece of the response field, which is the ingredient of ψ The tangent component of ∆ψ0 is projected off the solution in eq 11. It appears only as the extra DLP term A9. This exercise reveals the origin and importance of the condition in eq A7. It is seen how the PCM result can be regenerated starting from the SBCM methodology when eq A7 is violated. The procedure applying eq A10 is entirely equivalent to the standard integral equation technique that realizes the boundary conditions of the PCM procedure.2,32

The comments listed below briefly outline several specific features of the illustrative calculations reported in section A3. (a) The integral representation of the response field Φ, as suggested by eqs 9 and 11, for the particular case of an Onsager dipole has the structure

2 Φ(R, Ω) ) - mR cos θ 3

dε ∫a∞ dR' ( R'1 ) ε12 dR' 3

(A12) (A9)

where ds′ means the surface differential element expressed in terms of r′. It represents the electric field created by the singular dipole density p distributed over the surface S; ∂/∂n′ denotes the normal derivative for r′. The following boundary conditions are valid30

ξ+ - ξ- ) 4πp,

(R < a)

A4. Technical Commentary

The arising discrepancy is a consequence of violating the condition in eq A7. In the case of a point ionic solute, the similar computation gives the standard Born result in terms of both SBCM and PCM algorithms, because the condition in eq A7 continues to be valid in that case. In the general abrupt stepwise case (ε+ * ε-), when eq A7 is violated, the correction terms in eq A4 are expressed in terms of the double-layer potential as

ξ)

ψ ) ψ0 + ξ

(A10)

∂ (ξ ∂n + ξ-) ) 0

The matching condition providing the smooth tangent derivative of ψ, namely, (∂/∂t)(ψ+ - ψ-) ) 0 for the total function A4, produces the final expression

where ε ) ε(R′). This is gained after applying the expansion over spherical harmonics for 1/|r - r′|.33 Here, the position variables are denoted as r and r′, with R ) |r| and R′ ) |r′|; Ω and θ stand for angular variables. When ε(R′) tends to be a step function at R′ ) a, then dε/dR′ represents a narrow delta peak, and 1/R′3 can be substituted by 1/a3. Thereby, the integral in eq A12 reduces to

1 a3

dε 1 ε dR' ) 3 ∫1 ∫a∞ ε12 dR' a

0

(

dε 1 1 ) 3 12 ε0 ε a

)

(A15)

The change of variables [R′ f ε(R′)] is made before the step limit transition. The monotonic variation of ε(R′) is implied, according to eqs 15 and 17.

2466

J. Phys. Chem. B, Vol. 114, No. 7, 2010

(b) The solution of eq A10 again invokes the spherical expansion for 1/|r - r′|. In the derivation of this basic equation, special properties of the DLP (listed after eq A9) were used. For the similar equation derived from the continuity conditions for the normal derivative, namely, ε+(∂ψ+/∂n) - ε-(∂ψ-/∂n) ) 0, the solution vanishes, that is, the singular SLP component is absent. This note complements the pertinent remark appearing after eq A11. References and Notes (1) Continuum SolVation Models in Chemical Physics; Mennucci, B., Cammi, R. Eds.; Wiley: Chichester, U.K., 2007. (2) Tomasi, J.; Mennucci, B.; Cammi, R. Chem. ReV. 2005, 105, 2999. (3) Cramer, C. J.; Truhlar, D. G. Chem. ReV. 1999, 99, 2161. (4) Baker, N. A. In ReViews in Computational Chemistry; Lipkovitz, K. B., Larter, R., Cundari, T. R., Eds.; Wiley: New York, 2005; Vol. 21, p 349. (5) Lu, B. Z.; Zhou, Y. C.; Holst, M. J.; McCammon, J. A. Commun. Comput. Phys. 2008, 3, 5–973. (6) Li, J.; Zhu, T.; Hawkins, G. D.; Winget, P.; Liotard, D. A.; Cramer, C. J.; Truhlar, D. G. Theor. Chem. Acc. 1999, 103, 9. (7) Klamt, A.; Schuurman, G. J. Chem. Soc., Perkin Trans. 2 1993, 799. Klamt, A. J. Phys. Chem. 1995, 99, 2224. (8) Onzager, L. J. Am. Chem. Soc. 1936, 58, 1486. (9) Kirkwood, J. C. J. Chem. Phys. 1934, 2, 357. (10) Kharkats, Yu. I.; Kornyshev, A. A.; Vorotyntsev, M. A. J. Chem. Soc., Faraday Trans 2 1976, 72, 361. (11) Kornyshev, A. A.; Vorotyntsev, M. A. The Electrostatics of the Media with Spatial Dispertion; Nauka: Moscow, Russia, 1993. (12) Cossi, M.; Mennucci, B.; Tomasi, J. Chem. Phys. Lett. 1994, 228, 165. (13) Frediani, L.; Pomelli, C. S.; Tomasi, J. Phys. Chem. Chem. Phys. 2000, 2, 4876. Frediani, L.; Cammi, R.; Corni, S.; Tomasi, J. J. Chem. Phys. 2004, 120, 3893. 8. (14) Discussion of earlier solvation models with variable ε, as applied to spherically symmetrical solutes (ions), can be found in refs 2, 12, and 15. (15) Basilevsky, M.; Odinokov, A.; Nikitina, E.; Grigoriev, F.; Petrov, N.; Alfimov, M. J. Chem. Phys. 2009, 130, 024504–024505. (16) Kirkwood, J. C. J. Chem. Phys. 1939, 7, 911. (17) Frohlich, H. Theory of Dielectrics: Oxford University Press: New York, 1949. (18) Booth, F. J. Chem. Phys. 1951, 19, 391. (19) More complicated density functions were constructed in ref 15. They are required for reproducing properly the complicated solvation effects

Basilevsky et al. inherent to small spherical ions. The present approach addresses larger organic molecules for which such fine effects can be neglected. (20) Chudinov, G. E.; Napolov, D. V.; Basilevsky, M. V. Chem. Phys. 1992, 160, 41. (21) Jorgensen, W. L.; Ulmschneider, J. P.; Tirado-Rivers, J. J. Phys. Chem. B 2004, 108, 16264. (22) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J. W.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J. J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, M. C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision B.03; Gaussian Inc.: Pittsburgh, PA, 2004. (23) Gas-phase optimized geometries were used. The solute charges were computed within the PCM option for the pertaining solvent (water in section and other solvents in section 4.2). (24) Ben-Naim, A. SolVation Thermodynamics; Plenum Press: New York, 1987. (25) Chambers, C. C.; Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. 1996, 100, 16385. Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 1997, 101, 7147. 1998, 102, 3257. (26) Hummer, G.; Pratt, L. R.; Garcia, A. E. J. Phys. Chem. 1996, 100, 1206. Hummer, G.; Pratt, L. R.; Garcia, A. E. J. Am. Chem. Soc. 1997, 119, 8523. (27) Fedorov, M. V.; Kornyshev, A. A. Mol. Phys. 2007, 105, 1. (28) Johnson, M. D.; Miller, J. R.; Green, N. S.; Closs, G. L. J. Phys. Chem. 1989, 93, 1173. (29) Leontyev, I. V.; Tovmash, A. V.; Vener, M. V.; Rostov, I. V.; Basilevsky, M. V. Chem. Phys. 2005, 319, 4. (30) Korn, G. A.; Korn, T. M. Mathematical Handbook for Scientists and Engineers; McGraw-Hill: New York, 1968; Section 15.6.5. Courant, R. Partial Differential Equations; Wiley-Interscience: New York, 1962. (31) The formulation of the D-G0 theorem in ref addressed only the discontinuous case ε+ * ε- in terms of the restrictive condition in eq A7. (32) Basilevsky, M. V.; Chudinov, G. E. Chem. Phys. 1991, 157, 327. Cances, E.; Mennucci, B.; Tomasi, J. J. Chem. Phys. 1997, 107, 3032. (33) Jackson, J. D. Classical Electrodynamics, 3rd ed.; Wiley: New York, 1999.

JP9064399