9056
J. Phys. Chem. 1995,99, 9056-9061
Simplifications in the Theoretical Treatment of Solvation Effects Based on the Quadratic Expansion of the Energy Alessandro Fortunelli Istituto di Chimica Quantistica ed Energetica Molecolare del C.N.R., Via Risorgimento 35, I-56126 Pisa, Italy Received: January 12, 1995; In Final Form: March 23, 1995@
An analysis of the electrostatic free energies of a set of neutral test molecules in water is presented, as calculated within the polarizable continuum model (PCM). The analysis is based on the formalism of the quadratic expansions of the energy. A general expression is derived for the solvation energy and a sequence of approximations to it, ranging from the rigorous linear response, to the linear rigime, to an approximate linear response, to the zeroth-order response. It is shown that the linear response of the solvated molecule to the solvent reaction field potential is of sufficient accuracy for most purposes, whereas the zeroth-order response may not be. The merits of each approximate expression are also compared, showing that the use of an approximate linear response allows one to drastically reduce the computational effort associated with the treatment of solvation effects to a negligible portion of a standard quantum chemical calculation, without degrading the overall accuracy.
1. Introduction The theoretical simulation of chemical processes in solution is difficult due to the large number of solvent molecules to be considered. This hinders a pure quantum mechanical approach to the study of solvated systems and makes the use of simplified methods necessary. Among them, the polarized continuum model (PCM), first proposed by MiertuS, Scrocco, and Tomasi,' has proven to be a reliable tool for the description of electrostatic solute-solvent interactions (see ref 2 for a comprehensive review of the continuum models for solvation). In this approach, a quantum mechanical description of the solute and a quasicontinuum representation of the solvent are used. The PCM was originally implemented at the Hartree-Fock (HF) level of quantum mechanical description of the solute' and successively generalized to include most of the methods of quantum chemistry: MC-SCF,3 CI,3,4MBPT,5 density functional (DF),6 semiempirical,' etc. methods have been explicitly considered within the approach. In the present article an analysis of the PCM electrostatic free energies of a set of neutral test molecules in water is performed, based on the formalism of the quadratic expansions of the energy. A general expression will be derived for the solvation energy and a sequence of approximations to it. In particular, it will be shown that the linear response of the solvated molecule to the reaction field potential due to the solvent is sufficient to give accurate values of interaction free energies. Moreover, further simplifications will be derived in the linear rigime, which drastically reduce the computational effort associated with the treatment of solvation effects to a negligible portion of a standard quantum chemical calculation, without degrading the overall accuracy. These simplifications may be helpful to study the behavior of large molecules in solution. 2. Method
The Polarizable Continuum Model. The PCM is based on the theory of electrostatic interactions in fluids. It assumes that the solvent is a continuum dielectric, which reacts against the @Abstractpublished in Advance ACS Absfracts, May 1, 1995.
solute charge distribution (both electronic and nuclear), generating a reaction field. The effect of this reaction field is introduced into the solute Hamiltonian (&) by means of a perturbation operator (V,) according to
&Y = (&(O)
+ qR)Y = E Y
(1)
where &(') stands for the unperturbed gas-phase Hamiltonian. A molecular-shaped algorithm is used to define the cavity in which the solute is placed and which is embedded in the infinite polarizable dielectric medium. The polarization of the solvent by the solute charge distribution induces a surface charge distribution, a(s), on the solutelsolvent interface. The reaction field generated by a(s) is included into the Hamiltonian .@cO) by means of the perturbation operator V, as in (1). It the solutelsolvent interface is divided into M elements, S,, small enough to consider the charge distribution a(s) constant inside them, the perturbation operator can be further approximated as
where qi are the polarization charges on the cavity surface. Evaluation of the charge density at the different surface elements, ~ ( s , ) ,is performed by solving the Laplace equation
(3) where VTis the total electrostatic potential, which includes the contribution of both solvent (apparent) and solute charge distributions. A dielectric constant of 78.39 was used in the present calculations to mimic the properties of pure water at 298 K. Since the potential due to a(s) can be extracted from VT and taken over to the lhs of (3), a(s) itself comes out to be a linear (though intricate) function of the solute potential, and therefore of the solute charge distribution (for fixed e). The algorithm utilized to build up the cavity defines each atomic sphere in terms of a pentakisdodecahedron, whose faces are progressively divided in triangles, defining in this way a sequence of levels of approximation (see ref 8 for more details).
0022-365419512099-9056$09.0010 0 1995 American Chemical Society
Theoretical Treatment of Solvation Effects The present calculations were carried out at the standard surface cavity definition (level 1 in ref 8), except for the resserae at the intersection of two sphere, for which the largest surface cavity definition was used (level 5 in ref 8). Alternative definitions or refinements are discussed in the literature.* The cavity was determined from the atomic van der Waals radii multiplied by a scaling factor of 1.2. Standard van der Waals radii were used: H, 1.2; C, 1.7; N, 1.6; 0, 1.5 (see refs 9 and 10 for more refined definitions of the atomic radii). Following linear free response theory, the electrostatic contribution to the free energy (AGele)was computed according to the formula (see ref 2 for details)
J. Phys. Chem., Vol. 99, No. 22, 1995 9057
Theory of the Quadratic Expansions of the Energy. Let us recall some well-known equations of the quadratic expansion theory. Let the energy E of the system be considered as a function of some variational parameters {A}:E(A). Let then E be written as a sum of an unperturbed term, EO, and a small perturbative correction, U. Defining {Ao} as the set of parameters for which EO has a minimum and therefore the gradient of EOwith respect to A vanishes, v~Eo(A0)= 0, it is possible to derive an expression for the energy difference, AE, between Eo(A0) and the new energy minimum E(& obtained after introducing U into the energy expression:
AE = E(A) - Eo(Ao)
(5)
where A are the values of the parameters {A} which minimize Eo U. By first expanding E ( A ) up to quadratic terms in AA = & . Ao, one has
+
where pnucis the nuclear charge density, Vgis the elec-trostatic potential generated by the solvation charges qi, and K a n d Y (respectively, 5%") and W0))are the solute Hamiltonian and wave function in the solvent (respectively, in vacuo). According to this definition, the third term on the rhs of (4) is included to account for the work involved in the polarization of the solvent (see, e.g., ref 11, pp 1,41-142). As noted above, F4 is, a liqear ope;ator in the space of one-electron functions: % $1 [e]. % is a symmetric operator (see ref 11, p 143), however-due to numerical approxjmations-in the present calculations small deviations of Vo from hermiticity were found (see below): Jdr @A s"[@B] f Jdr @E V k ' A l .
Apart from the electrostatic component, the total Gibbs free energy of hydration AGy&is evaluated in the PCM approach2 by adding a "steric" (or, better, nonelectrostatic) contribution AG,,,, accounting for the van der Waals interactions and the work needed to disrupt the solvent structure and build the solute cavity (the cavitation free energy). Other thermal contributions to AGydrare not usually considered. In the present work, eqs 1-4 were solved within a KohnSham density functional (KS-DF) scheme,l2-I4 since DF calculations have proven to be able to furnish accurate values of AGydrfor neutral molecules in water.6 A standard convergence procedure was applied: first the in vacuo calculation was performed, then the electrostatic potential at the surface was evaluated, and the corresponding polarization charges were calculated by using the iterative solution of the self-consistent equations (see ref 2): they were then introduced into the standard Kohn-Sham Hamiltonian according to (1 and 2) and the whole process was iterated until convergence. Three to five iterations were necessary (as usual) to achieve convergence up to The sum of all the self-consistent cycles (including the in vacuo part) was 2-5 times the number necessary for a simple in vacuo calculation. The evaluation of the surface charges always implied a negligible CPU time with respect to the calculation of the electronic energy. The same applies to the evaluation of the solute potential and of the additional one-electron integrals necessary for the description of the field generated by the surface charges. Altogether, they represented a few percent of the total CPU time (see below for a more detailed discussion). In the following, e(o) and V(O)= V[e(O)] will denote the gasphase charge density (both electronic and nuclear) and the corresponding reaction field potential, respectively. Accordingly, $') and @') = V[Q(')] will denote the same quantities after the first self-consistent iteration etc. until @(fland Vcfl], which are the same quantities at convergence.
E(A) = E(&)
+ VAE(Ao).AA + '/2AAt*VV,E(Ao).AA
(6)
where VVAE(AO)is the (symmetric) Hessian matrix at Ao. By imposing moreover that VAE(A)= 0, one obtains
VAE(&)
+ AA~~VV,E(A,)= o
(7)
and thus
E(.&) = E@,) 4-'/,VAE(A0).AA
(8)
Now, since VAEO(AO) = 0, one finally obtains
AE = E(A) - Eo(Ao)= U(&)
+ '/,VAU(A0).AA
(9)
It is to be noted that EOand its derivatives do not formally appear in the expression (9) for AE,which depends explicitly only on U. Moreover, defining AEpol= EO(& - Eo(A0) as the solute polarization energy (i.e., the difference in the intemal energy of the solute molecule due to solvation effects), one finds
(an alternative definition of polarization energy, which should not be confused with the present one, can be found in ref 15). In passing, we note that if A does not differ too much from A or Ao, one also has
E(A) = E(A)
+ 1 / 2 (-~ A)+.VV,E(A,>(A
- A) (11)
Let us specify (5)-( 11) to the PCM case. The parameters {A} are left unspecified, but it is assumed that the electron density e is linear in A so that VA,.AA = A@,where A@is the change due to the solvent effects: this is tantamount to neglecting the Hessian of e with respect to the parameters {A}. Within the PCM approach one has
u = '12Jdr e@>V[el(r)
(12)
(9) thus gives
+ + '/,Sdr
Actf,, = 'I2Jdr e"'(r) V[e'o'](r)
1/4j d r e"'(r) V[Ae"](r)
Ae"(r) %'[~(~)](r) (13)
(if one does not neglect the Hessian of @, one should add a term -l/4 Jdr [AA+.VVA~(Ao).AA](r).V[~(o)](r) to the rhs of this equation).
9058 J. Phys. Chem., Vol. 99, No. 22, 1995
Fortunelli
If one considers, instead, U as a modification of the oneelectron operator, U = Jdr @(r)VL'(r), where V~'(r) is considered as an external potential (note that the factor disappears as U now represents an energy, not a free energy contribution), one has that the difference in the solute energy due to the solvent, AE, is given by
+ l/,Sdr Ae"(r)
@(r)V[@](r)), and, by also exploiting the fact that V is a symmetric operator, one obtains
A
