- Email: [email protected]
The implementation of density functional theory within the polarizable continuum model for solvation
16 December 1994
CHEMICAL PHYSICS LETTERS Chemical Physics Letters 231 ( 1994) 34-39
The implementation of density functional theory within the polarizable continuum model for solvation Alessandro Fortunelli, Jacopo Tomasi Istituto di Chimica Quantistica
ed Energetica Molecolare
de1 CNR, Via Risorgimento
35, 56126 Piss, Italy
Received 1.5 September 1994; in final form 10 October 1994
Abstract The implementation of density functional theory within the polarizable continuum model for solvation is reported. Test calculations are performed on a set of representative compounds and the resulting free energies of hydration, dGhydrr are compared with experimental data and Hartree-Fock calculations. Two gradient-corrected functionals are considered and are found to assure an improved description of the solute-solvent electrostatic interactions with respect to the Hartree-Fock results. The mean square root deviation of the AGhydr values with respect to experimental values is found to be only 0.78 kcal/mol (respectively 0.77 kcal/mol) utilizing the Becke functional for exchange and the Perdew functional (respectively the Lee-Yang-Parr functional) for correlation, whereas it is 1.97 kcal/mol at the Hartree-Fock level.
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 [ 11, has proven to be a reliable tool for the description of electrostatic solute-solvent interactions (see Ref. [ 21 for a comprehensive review of the continuum models for solvation). In this approach, a quantum mechanical description of the solute and a quasi-continuum representation of the solvent are used. Since its proposal, the PCM has undergone a number of modifications and improvements. Among the most recent ones, we quote the possibility of treating non-homogeneous and/or non-isotropic solvents [ 31, which is of interest, e.g., in the biological field: 0009-2614/94/$07.00
the description of proteins, counter-ion effects, etc. The PCM was originally implemented at the Hartree-Fock (HF) level of quantum mechanical description of the solute [ 11, and the introduction of correlation effects has always followed the lines of standard ab initio quantum chemistry: MC-SCF [4], CI [ 51, MBPT [ 61, etc. methods have been explicitly considered within the approach. On the other hand, the merits of density functional (DF) theory, by far the most useful non-empirical alternative to conventional post-HF methods for studying physico-chemical properties of molecules, are nowadays well recognized [7-lo]. Gradientcorrected functionals have significantly increased the reliability of DF methods, and recent work has shown that small modifications of the functional form [lo] or partial inclusion of the HF exchange [ 111 can provide even better results. It is therefore of theoretical and practical interest to check the performances of DF-based approaches within the context of the polar-
@ 1994 Elsevier Science B.V. All rights reserved
.SSDIOOO9-2614(94)01253-9
A. Fortunelli, J. Tomb/Chemical
izable continuum model for solvation, to see whether they can represent also in this context a viable method for describing in an approximate way correlation effects within the solute. This is the aim of this Letter. The use of DF theory within continuum models for solvation has already been considered in the literature (see, e.g., Refs. [12,13]). In Ref. [12], in particular, DF theory is implemented within a self-consistent reaction field model which has several bearings with the PCM. The main differences lie in the fact that in Ref. [ 121 the electrostatic quantities of interest are obtained through a multipole expansion, whereas the PCM methodology is based on the boundary element method within the apparent surface charge approach (see Ref. [2] for details) and should therefore be more versatile in terms of the description of general cavities and an accurate representation of the interaction potential (we refer to Ref. [ 21 for a comparative discussion of these and other continuum approaches to solvation) .
2. The method 2.1. 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 solute charge distribution, generating a reaction field. The effect of this reaction field is introduced into the solute Hamiltonian by means of a perturbation operator (i& ) according to (“ria+Va)w.=EV,
(1)
where fia 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 solute/solvent interface. The reaction field generated by c+(s) is included into the Hamiltonian 7%~by means of the perturbation operator & as in (1). If the solute/solvent interface is divided into M elements, S,, small enough to consider the charge distribution (T(S) constant inside them, the perturbation operator can be further approximated as
Physics Letters 231 (1994) 34-39
M =gg+
(2)
where qi are the polarization charges on the cavity surface. Evaluation of the charge density at the different surface elements, I, is performed by solving the Laplace equation (3) where VT is the total electrostatic potential, which includes the contribution of both solvent 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. 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. [ 141 for more details). The present calculations were carried out at the largest surface cavity definition (level 5 in Ref. [ 141). Alternative definitions or refinements are discussed in the literature [ 21. 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. Following linear free response theory, the electrostatic contribution to the free energy ( AG,,) was computed according to the formula (see Ref. [2] for details)
where pnuci is the nuclear charge density, V, is the electrostatic potential generated by the solvation charges qi, fi and * (respectively fia and Wa) are the solute Hamiltonian and wavefunction in the solvent (respectively, in vacua) . According to this definition, the third term in (4) is included to account for the work involved in the polarization of the solvent. Apart from the electrostatic component calculated through the PCM, the total Gibbs free energy of hy-
36
A. Fortunelli. J. Tomasi/Chemical Physics Letters 231 (1994) 34-39
dration hChydr was evaluated by adding a steric contribution AGsk, accounting for the van der Waals interaction (distinguished into a repulsive AG,, and an attractive AGdisp term) and the work needed to disrupt the solvent structure and build the solute cavity (the cavitation free energy AGCaY). Other thermal contributions to Achydr were not considered explicitly. The van der Waals contributions were determined with the method proposed in Ref. [ 151, which gives the interaction as a sum of an atom-atom exponential term for the repulsive part and an atom-atom rP6 term for the dispersive part. The cavitation free energy was determined from the scaled particle theory according to Pierotti’s method, as modified by Langlet et al. [ 161. The only parameters utilized are the van der Waals radii, whose value is given above, the following numerical coefficients (which determine the c,~ coefficients, see Ref. [ 151 for more details) : H 1.O; C 1.5; N 1.lS; 0 1.36, and some data characteristic of the solvent (water) : number density, etc., which were taken from Ref. [ 171. 2.2. The density functional
approach
All DF computations were performed within the Kohn-Sham (KS) formalism [ 18,191 using the GAUSSIAN 92/DFI system of programs [20]. Among the characteristics of this code we mention the use of Gaussian basis functions, the avoidance of auxiliary functions, the implementation of pruned grids, and the availability of analytical first and second derivatives. The numerical integration procedure applied in the calculations was developed by Becke [ 211. On the basis of some test computations with large grids, the standard SG-1 grid [22] was selected for the present study. As is well known, the main problem with DF approaches is that the exact functional dependence of the exchange and correlation energies, E, and EC, on the density is not known, and one must resort to approximate expressions for them in practical calculations. The reliability of KS calculations thus strictly depends on the approximation for E, + EC. In the present work, an option of the program was selected which utilizes the functional of Becke [ 231 for exchange and either the Perdew [24] or the Lee-Yang-Parr [25] functional for correlation. These belong to the set of socalled non-local or (better) gradient-corrected func-
tionals, which have been shown [7] to provide improved results of formation energies, etc. with respect to the simplest approximation for E, + EC, namely the local-spin-density approximation [ 261. It must be noted that, since this is an exploratory work, only wellestablished functional forms have been considered (no hybrid forms as, for example, in Ref. [ 111) . All the computations were performed using the Huzinaga (5/3) (for hydrogen), (9, 5/5, 2) (for second row atoms) basis set [27], augmented by a single set of polarization functions, The contraction coefficients were taken from atomic calculations: only the innermost primitives were contracted and a scale factor of 1.2 was then introduced for the hydrogen s-functions [ 281. The GAUSSIAN 92/DFT set of programs was purposely modified and linked to the PCM routines developed by the Pisa group ’ (the program can also work with all the ab initio methods programmed in the GAUSSIAN 92/DFT code: restricted and unrestricted Hartree-Fock, Moller-Plesset perturbation, coupled cluster, quadratic configuration interaction, etc. calculations). A standard convergence procedure was applied: first the in vacua calculation was performed, then the electrostatic potential at the surface was evaluated and the corresponding polarization charges were calculated using the iterative solution of the self-consistent equations (see Ref. [ 21) : they were then introduced into the standard KS 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 lop6 Eb_ The sum of all the self-consistent cycles (including the in vacua part) was two to five times the number necessary for a simple in vacua calculation. The ratio between DF and HF CPU times was heavily dependent on the molecule under consideration; however, the DF CPU time never exceeded that of HF by a factor larger than 2-2.5. 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. 1 Courtesy
of Dr. Maurizio
Cossi.
A. Fortunelli, J. Tomasi/Chemical
We expect DF theory to assure an improved description of the electrostatic interactions in solution over the HF approach. As is well known [7], for molecules in vacua DF approaches are able to produce both dipole (in general, multipole) moments and static polarizabilities in better agreement with experiment than those of HE and of the same quality of the most refined post-HF methods. It is hardly necessary to be reminded that quantities such as multipole moments and polarizabilities are crucial in determining the extent of solute-solvent interactions.
3. The results In Table 1, the results of PCM-DF calculations are reported for a set of test molecules, together with those of HF calculations, the experimental data and the values of cavitation and van der Waals energies computed as explained above. The notation is as follows. AGr$!“d is the total Gibbs free energy of hydration according to the given methodology or experiment (method) : when from theoretical calculations, AGf;dtOd is the sum of an electrostatic component A Gsztiod and a steric component AGster, which, in turn, is the sum of the cavitation AGcav, dispersive AGdisr and repulsive AGrep van der Waals interaction free energies. As for the notation ‘method’, exp stands for experiment, HF stands for Hartree-Fock, BP86 (respectively, BLYP) stands for the DF approach using the functional of Becke [23] for exchange and the functional of Perdew [ 241 (respectively, Lee-YangParr [ 2.51) for correlation. In Table 1 the values of the in vacua dipole moments, dipmethod, are also reported as obtained from the various methodologies or experiment. The test molecules were chosen among simple compounds for which experimental data for AGt,,,dr were accessible, and each of them is representative of a class of substances. The geometries at which the computations were performed are the experimental ones (for the molecules in vacua) and were taken from Ref. [ 291 (or from Ref. [ 171 if not present in Ref. [ 291) . It may be questioned as to whether the choice of gasphase geometries is correct when considering solvation processes. However, previous experience (see, e.g., Ref. [ 3 1 ] ) suggests that AGhydr for the molecules considered should be insensitive to geometrical mod-
Physics Letters 231 (1994) 34-39
31
ifications in solutions. From the analysis of the results of Table 1 the following observations and conclusions can be drawn: (a) There is a substantial agreement between BP86 and BLYP methodologies, the discrepancies between the two sets of results always being smaller than the theoretical accuracy. (b) The DF values of AChydr are close to the experimental ones: the mean square root deviation of the theoretical versus the experimental values is 0.78 kcal/mol for the BP86 and 0.77 kcal/mol for the BLYP methodologies, while the mean absolute error is 0.68 and 0.66 kcal/mol, respectively. With the exception of CHsCN ( 1.16 kcal/mol for BPS6 and 1.05 kcal/mol for BLYP) and CHsNH2 ( 1.57 kcal/mol for BP86 and 1.73 kcal/mol for BLYP) molecules, all DF errors are smaller than 1 kcal/mol. The largest errors may be attributed to the presence of hydrogen atoms bound to electronegative elements, which give hydrogen bonding and thus probably require a smaller van der Waals radius for describing their electrostatic interactions. (c) The contribution of AG,,,, is crucial in bringing the theoretical results of AGhydr into close agreement with the experimental data. (d) As previously known (see, e.g., Ref. [7] ), the DF gas-phase dipole moments are in much better agreement with the experimental ones than those of HF: the mean square root deviation is 0.15 D for both BP86 and BLYP while it is 0.35 D for HF. (e) The DF results of AGt,,,dr are generally superior to the HF ones, for which the mean square root deviation with respect to the experimental values is 1.97 kcal/mol, while the mean absolute error is 1.73 kcal/mol. This is due to the improved representation of the electronic charge distribution and of its polarizability (as discussed above). The only exceptions are, again, CHjNH2 and (in a less significant way) CH30H. As a final conclusion, it can be asserted that, through the DF approach, AGhydr values are obtained which are close to the experimental values at low computational cost (minimum amount of CPU and/or disk-space facilities). It can finally be noted that here we have considered only calculations of a ‘standard’ type: we have computed just energies, at a fixed geometry, with an iterative convergence procedure and in a homogeneous,
38
A. Fortune&
J. Tomasi/Chemical
Physics Letters 231 (1994) 34-39
Table 1 Hydration free enthalpies (AC, in kcal/mol) and gas-phase dipole moments (dip, in D) for a set of test compounds. The free enthalpies are distinguished according to the various steric (ster) - further distinguished into van der Waals dispersive (disp), cavitation (cav) and van der Waals repulsive (rep) -, electrostatic (ele) and total (hydr) contributions and to the method used for their evaluation: experiment (exp), Hartree-Fock (HF), DF with the Becke functional [23] for exchange and the Perdew [24] (BP86) or the Lee-Yang-Parr [25] (BLYP) functional for correlation (see text for further details). The dipole moments are denoted accordingly. Experimental data of AGhydr were taken from Refs. [ 29,301, of the dipole moments from Ref. [ 171
CH3CH3 CH3CH0 CH3CN CH3COCH3 CHxCONH2 CH3COOCH.j CH3COOH CHjNH2 CH30CH3 CH3OH H20
HCN NH3
-8.63 -9.08 -8.53 -11.31 -10.95 -12.82 -10.48 -8.25 -10.15 -7.82 -4.91 -5.60 -5.60 ~@P86 ele
CHjCH3 CHjCHO CH3CN CH3COCHj CH3 CONH2 CHjCOOCHj CHjCOOH CHjNH2 CH30CH3 CHjOH H20
HCN NH3
-0.53 -5.96 -7.39 -6.25 -11.60 -6.24 -8.53 -5.28 -4.02 -6.09 -8.15 -6.00 -6.08
9.80
1.34
9.83 9.40 12.85 11.84 14.02 11.25 9.19 11.16 8.18 4.99 6.36 5.85
1.49 1.46 1.73 1.69 1.93 1.67 1.31 1.55 1.31 0.96 1.09 1.01
&BPS6 hydr
+ 1.98 -3.72 -5.06 -2.99 -9.02 -3.11 -6.08 -3.03 -1.47 -4.44 -7.11 -4.16 -4.8 1
2.51 2.24 2.33 3.27 2.58 3.13 2.44 2.25 2.56 1.67 1.04 1.85 1.26 AC?hydr +1.8 -3.5 -3.9 -3.8 -9.7 -3.3 -6.7 -4.6 -1.9 -5.1 -6.3 -3.2 -4.3
solvent. It is not difficult, however, to introduce in the code all the recent improvements of the PCM approach: simultaneous convergence of wavefunction and surface charges [ 321, analytical derivatives of the energy with respect to geometrical parameters or external fields [ 331, non-homogeneous and/or non-isotropic solvents [ 31, etc. (such extensions are in progress). isotropic
-0.44 -8.97 -8.68 -9.53 -14.66 -8.59 -11.36 -5.75 -5.04 -7.20 -8.98 -7.84 -6.37
+ 2.07 -6.72 -6.35 -6.27 -12.08 -5.46 -8.91 -3.50 -2.48 -5.54 -7.94 -5.99 -5.11
-0.47 -5.94 -7.28 -6.64 -11.46 -6.21 -8.46 -5.12 -3.96 -6.01 -8.08 -5.90 -5.98
+ 2.03 -3.69 -4.95 -3.38 -8.87 -3.07 -6.01 -2.87 -1.41 -4.36 -7.04 -4.05 -4.71
dipHF
dipBLYP
dip BP86
dipexP
0.00 3.20 4.15 3.35 4.42 1.88 1.90 1.65 1.69 2.01 2.15 3.21 1.79
0.00 2.61 3.79 2.75 3.91 1.67 1.57 1.53 1.45 1.79 2.02 2.79 1.73
0.00 2.62 3.83 2.76 3.93 1.68 1.57 1.55 1.44 1.80 2.03 2.81 1.75
0.00 2.69 3.92 2.88 3.76 1.72 1.74 1.31 1.30 1.70 1.85 2.98 1.47
Acknowledgement
We thank Dr. Maurizio Cossi for having provided us with the most recent version of the PCM codes developed by the Pisa group. This research was carried out with the contribution of the “Progetto Finalizzato Materiali Speciali per Tecnologie Avanzate” of the Italian CNR.
References [ 1 ] S. MiertuS, E. Scrocco and J. Tomasi, Chem. Phys. 55 ( 1981) 117.
A. Fortunelli, J. Tomasi/Chemical [2] J. Tomasi and M. Persico, Chem. Rev., in press. [3] B. Mennucci, M. Cossi and J. Tomasi, J. Chem. Phys., to be published; M. Cossi, B. Mennucci and J. Tomasi, Chem. Phys. Letters 228 (1994) 165. [4] E.L. Coitiho, J. Tomasi and 0. Ventura, J. Chem. Sot. Faraday Trans. 90 ( 1994) 1745, and references therein. [5] M. Persico and J. Tomasi, Croat. Chim. Acta 57 (1984) 1395. [6] EJ. Olivares de1 Valle and J. Tomasi, Chem. Phys. 150 (1991) 139. [7] T. Ziegler, Chem. Rev. 91 (1991) 651. [ 8 I J. Andzelm and E. Wimmer, J. Chem. Phys. 96 ( 1992) 1280. 191 B.G. Johnson, l?M.W. Gill and J.A. Pople, J. Chem. Phys. 98 (1993) 5612. [ 101 G.J. Laming, V. Termath and N.C. Handy, J. Chem. Phys. 99 (1993) 8765. [ 111 A.D. Becke, J. Chem. Phys. 98 ( 1993) 1372, 5648. [ 121 M.F. Ruiz-Mpez, E Bohr, M.T.C. Martins-Costa and D. Rinaldi, Chem. Phys. Letters 221 (1994) 109. [ 131 R.R. Contreras, J. Mol. Struct. THEOCHEM 310 (1994) 295. [ 141 J.L. Pascual-Ahuir, E. Silla, J. Tomasi and R. Bonaccorsi, J. Comput. Chem. 8 (1987) 778. [ 151 FM. Floris and J. Tomasi, J. Comput. Chem. 10 (1989) 616; EM. Floris, J. Tomasi and J.L. Pascual-Ahuir, J. Comput. Chem. 12 (1991) 784. [ 161 R.A. Pierotti, Chem. Rev. 76 (1976) 717; J. Langlet, I? Claverie, J. Caillet and B. Pullman, J. Phys. Chem. 92 (1988) 1617. [ 171 D.R. Lide, ed., Handbook of chemistry and physics, 72th Ed. (CRC Press, Boca Raton, 1991). [ 181 I! Hohenberg and W. Kohn, Phys. Rev. A 136 ( 1964) 864.
Physics Letters 231 (1994) 34-39
39
[ 191 W. Kohn and L. Sham, Phys. Rev. A 140 (1965) 1133. [20] M.J. Frisch, G.W. Trucks, M. Head-Gordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robb, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. DeFmes, J. Baker, J.J.F! Stewatt and J.A. Pople, GAUSSIAN 92/DFI (Gaussian Inc., Pittsburgh, 1993). [ 211 A.D. Becke, J. Chem. Phys. 88 (1988) 2547. [22] F?M.W. Gill, B.G. Johnson and J.A. Pople, Chem. Phys. Letters 209 (1993) 506. [23] A.D. Becke, Phys. Rev. A 33 (1988) 2786. [24] J.P Perdew, Phys. Rev. B 33 (1986) 8822. [25] C. Lee, W. YangandR.G. Parr, Phys. Rev. B 37 (1988) 785. ]26] P.A.M. Dirac, Proc. Cambridge Phil. Sot. 26 ( 1930) 376: S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58 (1980) 1200. [27] S. Huzinaga, J. Chem. Phys. 42 (1965) 293. [28] J. AImlGf and P.R. Taylor, Advan. Quantum Chem. 22 (1991) 301. [29] H.A. Carlson, T.B. Nguyen, M. Orozco and W.L. Jorgensen, J. Comput. Chem. 14 (1993) 1240. [30] M. Orozco, W.L. Jorgensen and EJ. Luque, J. Comput. Chem. 14 (1993) 1498. [31] I. Tut%n, E. Silla and J. Tomasi, J. Phys. Chem. 96 ( 1992) 9043; I. Alkorta. H.O. Villar and J.J. Perez, J. Comput. Chem. 14 ( 1993) 620, and references therein. [32] H. Hoshi, M. Sakurai, Y. Inowa and R. Chujo, J. Chem. Phys. 87 (1987) 1107; R. Cammi and J. Tomasi, J. Comput. Chem., in press. [33] R. Cammi and J. Tomasi, J. Chem. Phys. 100 ( 1994) 7495; 101 (1994) 3888.
CHEMICAL PHYSICS LETTERS Chemical Physics Letters 231 ( 1994) 34-39
The implementation of density functional theory within the polarizable continuum model for solvation Alessandro Fortunelli, Jacopo Tomasi Istituto di Chimica Quantistica
ed Energetica Molecolare
de1 CNR, Via Risorgimento
35, 56126 Piss, Italy
Received 1.5 September 1994; in final form 10 October 1994
Abstract The implementation of density functional theory within the polarizable continuum model for solvation is reported. Test calculations are performed on a set of representative compounds and the resulting free energies of hydration, dGhydrr are compared with experimental data and Hartree-Fock calculations. Two gradient-corrected functionals are considered and are found to assure an improved description of the solute-solvent electrostatic interactions with respect to the Hartree-Fock results. The mean square root deviation of the AGhydr values with respect to experimental values is found to be only 0.78 kcal/mol (respectively 0.77 kcal/mol) utilizing the Becke functional for exchange and the Perdew functional (respectively the Lee-Yang-Parr functional) for correlation, whereas it is 1.97 kcal/mol at the Hartree-Fock level.
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 [ 11, has proven to be a reliable tool for the description of electrostatic solute-solvent interactions (see Ref. [ 21 for a comprehensive review of the continuum models for solvation). In this approach, a quantum mechanical description of the solute and a quasi-continuum representation of the solvent are used. Since its proposal, the PCM has undergone a number of modifications and improvements. Among the most recent ones, we quote the possibility of treating non-homogeneous and/or non-isotropic solvents [ 31, which is of interest, e.g., in the biological field: 0009-2614/94/$07.00
the description of proteins, counter-ion effects, etc. The PCM was originally implemented at the Hartree-Fock (HF) level of quantum mechanical description of the solute [ 11, and the introduction of correlation effects has always followed the lines of standard ab initio quantum chemistry: MC-SCF [4], CI [ 51, MBPT [ 61, etc. methods have been explicitly considered within the approach. On the other hand, the merits of density functional (DF) theory, by far the most useful non-empirical alternative to conventional post-HF methods for studying physico-chemical properties of molecules, are nowadays well recognized [7-lo]. Gradientcorrected functionals have significantly increased the reliability of DF methods, and recent work has shown that small modifications of the functional form [lo] or partial inclusion of the HF exchange [ 111 can provide even better results. It is therefore of theoretical and practical interest to check the performances of DF-based approaches within the context of the polar-
@ 1994 Elsevier Science B.V. All rights reserved
.SSDIOOO9-2614(94)01253-9
A. Fortunelli, J. Tomb/Chemical
izable continuum model for solvation, to see whether they can represent also in this context a viable method for describing in an approximate way correlation effects within the solute. This is the aim of this Letter. The use of DF theory within continuum models for solvation has already been considered in the literature (see, e.g., Refs. [12,13]). In Ref. [12], in particular, DF theory is implemented within a self-consistent reaction field model which has several bearings with the PCM. The main differences lie in the fact that in Ref. [ 121 the electrostatic quantities of interest are obtained through a multipole expansion, whereas the PCM methodology is based on the boundary element method within the apparent surface charge approach (see Ref. [2] for details) and should therefore be more versatile in terms of the description of general cavities and an accurate representation of the interaction potential (we refer to Ref. [ 21 for a comparative discussion of these and other continuum approaches to solvation) .
2. The method 2.1. 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 solute charge distribution, generating a reaction field. The effect of this reaction field is introduced into the solute Hamiltonian by means of a perturbation operator (i& ) according to (“ria+Va)w.=EV,
(1)
where fia 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 solute/solvent interface. The reaction field generated by c+(s) is included into the Hamiltonian 7%~by means of the perturbation operator & as in (1). If the solute/solvent interface is divided into M elements, S,, small enough to consider the charge distribution (T(S) constant inside them, the perturbation operator can be further approximated as
Physics Letters 231 (1994) 34-39
M =gg+
(2)
where qi are the polarization charges on the cavity surface. Evaluation of the charge density at the different surface elements, I, is performed by solving the Laplace equation (3) where VT is the total electrostatic potential, which includes the contribution of both solvent 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. 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. [ 141 for more details). The present calculations were carried out at the largest surface cavity definition (level 5 in Ref. [ 141). Alternative definitions or refinements are discussed in the literature [ 21. 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. Following linear free response theory, the electrostatic contribution to the free energy ( AG,,) was computed according to the formula (see Ref. [2] for details)
where pnuci is the nuclear charge density, V, is the electrostatic potential generated by the solvation charges qi, fi and * (respectively fia and Wa) are the solute Hamiltonian and wavefunction in the solvent (respectively, in vacua) . According to this definition, the third term in (4) is included to account for the work involved in the polarization of the solvent. Apart from the electrostatic component calculated through the PCM, the total Gibbs free energy of hy-
36
A. Fortunelli. J. Tomasi/Chemical Physics Letters 231 (1994) 34-39
dration hChydr was evaluated by adding a steric contribution AGsk, accounting for the van der Waals interaction (distinguished into a repulsive AG,, and an attractive AGdisp term) and the work needed to disrupt the solvent structure and build the solute cavity (the cavitation free energy AGCaY). Other thermal contributions to Achydr were not considered explicitly. The van der Waals contributions were determined with the method proposed in Ref. [ 151, which gives the interaction as a sum of an atom-atom exponential term for the repulsive part and an atom-atom rP6 term for the dispersive part. The cavitation free energy was determined from the scaled particle theory according to Pierotti’s method, as modified by Langlet et al. [ 161. The only parameters utilized are the van der Waals radii, whose value is given above, the following numerical coefficients (which determine the c,~ coefficients, see Ref. [ 151 for more details) : H 1.O; C 1.5; N 1.lS; 0 1.36, and some data characteristic of the solvent (water) : number density, etc., which were taken from Ref. [ 171. 2.2. The density functional
approach
All DF computations were performed within the Kohn-Sham (KS) formalism [ 18,191 using the GAUSSIAN 92/DFI system of programs [20]. Among the characteristics of this code we mention the use of Gaussian basis functions, the avoidance of auxiliary functions, the implementation of pruned grids, and the availability of analytical first and second derivatives. The numerical integration procedure applied in the calculations was developed by Becke [ 211. On the basis of some test computations with large grids, the standard SG-1 grid [22] was selected for the present study. As is well known, the main problem with DF approaches is that the exact functional dependence of the exchange and correlation energies, E, and EC, on the density is not known, and one must resort to approximate expressions for them in practical calculations. The reliability of KS calculations thus strictly depends on the approximation for E, + EC. In the present work, an option of the program was selected which utilizes the functional of Becke [ 231 for exchange and either the Perdew [24] or the Lee-Yang-Parr [25] functional for correlation. These belong to the set of socalled non-local or (better) gradient-corrected func-
tionals, which have been shown [7] to provide improved results of formation energies, etc. with respect to the simplest approximation for E, + EC, namely the local-spin-density approximation [ 261. It must be noted that, since this is an exploratory work, only wellestablished functional forms have been considered (no hybrid forms as, for example, in Ref. [ 111) . All the computations were performed using the Huzinaga (5/3) (for hydrogen), (9, 5/5, 2) (for second row atoms) basis set [27], augmented by a single set of polarization functions, The contraction coefficients were taken from atomic calculations: only the innermost primitives were contracted and a scale factor of 1.2 was then introduced for the hydrogen s-functions [ 281. The GAUSSIAN 92/DFT set of programs was purposely modified and linked to the PCM routines developed by the Pisa group ’ (the program can also work with all the ab initio methods programmed in the GAUSSIAN 92/DFT code: restricted and unrestricted Hartree-Fock, Moller-Plesset perturbation, coupled cluster, quadratic configuration interaction, etc. calculations). A standard convergence procedure was applied: first the in vacua calculation was performed, then the electrostatic potential at the surface was evaluated and the corresponding polarization charges were calculated using the iterative solution of the self-consistent equations (see Ref. [ 21) : they were then introduced into the standard KS 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 lop6 Eb_ The sum of all the self-consistent cycles (including the in vacua part) was two to five times the number necessary for a simple in vacua calculation. The ratio between DF and HF CPU times was heavily dependent on the molecule under consideration; however, the DF CPU time never exceeded that of HF by a factor larger than 2-2.5. 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. 1 Courtesy
of Dr. Maurizio
Cossi.
A. Fortunelli, J. Tomasi/Chemical
We expect DF theory to assure an improved description of the electrostatic interactions in solution over the HF approach. As is well known [7], for molecules in vacua DF approaches are able to produce both dipole (in general, multipole) moments and static polarizabilities in better agreement with experiment than those of HE and of the same quality of the most refined post-HF methods. It is hardly necessary to be reminded that quantities such as multipole moments and polarizabilities are crucial in determining the extent of solute-solvent interactions.
3. The results In Table 1, the results of PCM-DF calculations are reported for a set of test molecules, together with those of HF calculations, the experimental data and the values of cavitation and van der Waals energies computed as explained above. The notation is as follows. AGr$!“d is the total Gibbs free energy of hydration according to the given methodology or experiment (method) : when from theoretical calculations, AGf;dtOd is the sum of an electrostatic component A Gsztiod and a steric component AGster, which, in turn, is the sum of the cavitation AGcav, dispersive AGdisr and repulsive AGrep van der Waals interaction free energies. As for the notation ‘method’, exp stands for experiment, HF stands for Hartree-Fock, BP86 (respectively, BLYP) stands for the DF approach using the functional of Becke [23] for exchange and the functional of Perdew [ 241 (respectively, Lee-YangParr [ 2.51) for correlation. In Table 1 the values of the in vacua dipole moments, dipmethod, are also reported as obtained from the various methodologies or experiment. The test molecules were chosen among simple compounds for which experimental data for AGt,,,dr were accessible, and each of them is representative of a class of substances. The geometries at which the computations were performed are the experimental ones (for the molecules in vacua) and were taken from Ref. [ 291 (or from Ref. [ 171 if not present in Ref. [ 291) . It may be questioned as to whether the choice of gasphase geometries is correct when considering solvation processes. However, previous experience (see, e.g., Ref. [ 3 1 ] ) suggests that AGhydr for the molecules considered should be insensitive to geometrical mod-
Physics Letters 231 (1994) 34-39
31
ifications in solutions. From the analysis of the results of Table 1 the following observations and conclusions can be drawn: (a) There is a substantial agreement between BP86 and BLYP methodologies, the discrepancies between the two sets of results always being smaller than the theoretical accuracy. (b) The DF values of AChydr are close to the experimental ones: the mean square root deviation of the theoretical versus the experimental values is 0.78 kcal/mol for the BP86 and 0.77 kcal/mol for the BLYP methodologies, while the mean absolute error is 0.68 and 0.66 kcal/mol, respectively. With the exception of CHsCN ( 1.16 kcal/mol for BPS6 and 1.05 kcal/mol for BLYP) and CHsNH2 ( 1.57 kcal/mol for BP86 and 1.73 kcal/mol for BLYP) molecules, all DF errors are smaller than 1 kcal/mol. The largest errors may be attributed to the presence of hydrogen atoms bound to electronegative elements, which give hydrogen bonding and thus probably require a smaller van der Waals radius for describing their electrostatic interactions. (c) The contribution of AG,,,, is crucial in bringing the theoretical results of AGhydr into close agreement with the experimental data. (d) As previously known (see, e.g., Ref. [7] ), the DF gas-phase dipole moments are in much better agreement with the experimental ones than those of HF: the mean square root deviation is 0.15 D for both BP86 and BLYP while it is 0.35 D for HF. (e) The DF results of AGt,,,dr are generally superior to the HF ones, for which the mean square root deviation with respect to the experimental values is 1.97 kcal/mol, while the mean absolute error is 1.73 kcal/mol. This is due to the improved representation of the electronic charge distribution and of its polarizability (as discussed above). The only exceptions are, again, CHjNH2 and (in a less significant way) CH30H. As a final conclusion, it can be asserted that, through the DF approach, AGhydr values are obtained which are close to the experimental values at low computational cost (minimum amount of CPU and/or disk-space facilities). It can finally be noted that here we have considered only calculations of a ‘standard’ type: we have computed just energies, at a fixed geometry, with an iterative convergence procedure and in a homogeneous,
38
A. Fortune&
J. Tomasi/Chemical
Physics Letters 231 (1994) 34-39
Table 1 Hydration free enthalpies (AC, in kcal/mol) and gas-phase dipole moments (dip, in D) for a set of test compounds. The free enthalpies are distinguished according to the various steric (ster) - further distinguished into van der Waals dispersive (disp), cavitation (cav) and van der Waals repulsive (rep) -, electrostatic (ele) and total (hydr) contributions and to the method used for their evaluation: experiment (exp), Hartree-Fock (HF), DF with the Becke functional [23] for exchange and the Perdew [24] (BP86) or the Lee-Yang-Parr [25] (BLYP) functional for correlation (see text for further details). The dipole moments are denoted accordingly. Experimental data of AGhydr were taken from Refs. [ 29,301, of the dipole moments from Ref. [ 171
CH3CH3 CH3CH0 CH3CN CH3COCH3 CHxCONH2 CH3COOCH.j CH3COOH CHjNH2 CH30CH3 CH3OH H20
HCN NH3
-8.63 -9.08 -8.53 -11.31 -10.95 -12.82 -10.48 -8.25 -10.15 -7.82 -4.91 -5.60 -5.60 ~@P86 ele
CHjCH3 CHjCHO CH3CN CH3COCHj CH3 CONH2 CHjCOOCHj CHjCOOH CHjNH2 CH30CH3 CHjOH H20
HCN NH3
-0.53 -5.96 -7.39 -6.25 -11.60 -6.24 -8.53 -5.28 -4.02 -6.09 -8.15 -6.00 -6.08
9.80
1.34
9.83 9.40 12.85 11.84 14.02 11.25 9.19 11.16 8.18 4.99 6.36 5.85
1.49 1.46 1.73 1.69 1.93 1.67 1.31 1.55 1.31 0.96 1.09 1.01
&BPS6 hydr
+ 1.98 -3.72 -5.06 -2.99 -9.02 -3.11 -6.08 -3.03 -1.47 -4.44 -7.11 -4.16 -4.8 1
2.51 2.24 2.33 3.27 2.58 3.13 2.44 2.25 2.56 1.67 1.04 1.85 1.26 AC?hydr +1.8 -3.5 -3.9 -3.8 -9.7 -3.3 -6.7 -4.6 -1.9 -5.1 -6.3 -3.2 -4.3
solvent. It is not difficult, however, to introduce in the code all the recent improvements of the PCM approach: simultaneous convergence of wavefunction and surface charges [ 321, analytical derivatives of the energy with respect to geometrical parameters or external fields [ 331, non-homogeneous and/or non-isotropic solvents [ 31, etc. (such extensions are in progress). isotropic
-0.44 -8.97 -8.68 -9.53 -14.66 -8.59 -11.36 -5.75 -5.04 -7.20 -8.98 -7.84 -6.37
+ 2.07 -6.72 -6.35 -6.27 -12.08 -5.46 -8.91 -3.50 -2.48 -5.54 -7.94 -5.99 -5.11
-0.47 -5.94 -7.28 -6.64 -11.46 -6.21 -8.46 -5.12 -3.96 -6.01 -8.08 -5.90 -5.98
+ 2.03 -3.69 -4.95 -3.38 -8.87 -3.07 -6.01 -2.87 -1.41 -4.36 -7.04 -4.05 -4.71
dipHF
dipBLYP
dip BP86
dipexP
0.00 3.20 4.15 3.35 4.42 1.88 1.90 1.65 1.69 2.01 2.15 3.21 1.79
0.00 2.61 3.79 2.75 3.91 1.67 1.57 1.53 1.45 1.79 2.02 2.79 1.73
0.00 2.62 3.83 2.76 3.93 1.68 1.57 1.55 1.44 1.80 2.03 2.81 1.75
0.00 2.69 3.92 2.88 3.76 1.72 1.74 1.31 1.30 1.70 1.85 2.98 1.47
Acknowledgement
We thank Dr. Maurizio Cossi for having provided us with the most recent version of the PCM codes developed by the Pisa group. This research was carried out with the contribution of the “Progetto Finalizzato Materiali Speciali per Tecnologie Avanzate” of the Italian CNR.
References [ 1 ] S. MiertuS, E. Scrocco and J. Tomasi, Chem. Phys. 55 ( 1981) 117.
A. Fortunelli, J. Tomasi/Chemical [2] J. Tomasi and M. Persico, Chem. Rev., in press. [3] B. Mennucci, M. Cossi and J. Tomasi, J. Chem. Phys., to be published; M. Cossi, B. Mennucci and J. Tomasi, Chem. Phys. Letters 228 (1994) 165. [4] E.L. Coitiho, J. Tomasi and 0. Ventura, J. Chem. Sot. Faraday Trans. 90 ( 1994) 1745, and references therein. [5] M. Persico and J. Tomasi, Croat. Chim. Acta 57 (1984) 1395. [6] EJ. Olivares de1 Valle and J. Tomasi, Chem. Phys. 150 (1991) 139. [7] T. Ziegler, Chem. Rev. 91 (1991) 651. [ 8 I J. Andzelm and E. Wimmer, J. Chem. Phys. 96 ( 1992) 1280. 191 B.G. Johnson, l?M.W. Gill and J.A. Pople, J. Chem. Phys. 98 (1993) 5612. [ 101 G.J. Laming, V. Termath and N.C. Handy, J. Chem. Phys. 99 (1993) 8765. [ 111 A.D. Becke, J. Chem. Phys. 98 ( 1993) 1372, 5648. [ 121 M.F. Ruiz-Mpez, E Bohr, M.T.C. Martins-Costa and D. Rinaldi, Chem. Phys. Letters 221 (1994) 109. [ 131 R.R. Contreras, J. Mol. Struct. THEOCHEM 310 (1994) 295. [ 141 J.L. Pascual-Ahuir, E. Silla, J. Tomasi and R. Bonaccorsi, J. Comput. Chem. 8 (1987) 778. [ 151 FM. Floris and J. Tomasi, J. Comput. Chem. 10 (1989) 616; EM. Floris, J. Tomasi and J.L. Pascual-Ahuir, J. Comput. Chem. 12 (1991) 784. [ 161 R.A. Pierotti, Chem. Rev. 76 (1976) 717; J. Langlet, I? Claverie, J. Caillet and B. Pullman, J. Phys. Chem. 92 (1988) 1617. [ 171 D.R. Lide, ed., Handbook of chemistry and physics, 72th Ed. (CRC Press, Boca Raton, 1991). [ 181 I! Hohenberg and W. Kohn, Phys. Rev. A 136 ( 1964) 864.
Physics Letters 231 (1994) 34-39
39
[ 191 W. Kohn and L. Sham, Phys. Rev. A 140 (1965) 1133. [20] M.J. Frisch, G.W. Trucks, M. Head-Gordon, P.M.W. Gill, M.W. Wong, J.B. Foresman, B.G. Johnson, H.B. Schlegel, M.A. Robb, E.S. Replogle, R. Gomperts, J.L. Andres, K. Raghavachari, J.S. Binkley, C. Gonzalez, R.L. Martin, D.J. Fox, D.J. DeFmes, J. Baker, J.J.F! Stewatt and J.A. Pople, GAUSSIAN 92/DFI (Gaussian Inc., Pittsburgh, 1993). [ 211 A.D. Becke, J. Chem. Phys. 88 (1988) 2547. [22] F?M.W. Gill, B.G. Johnson and J.A. Pople, Chem. Phys. Letters 209 (1993) 506. [23] A.D. Becke, Phys. Rev. A 33 (1988) 2786. [24] J.P Perdew, Phys. Rev. B 33 (1986) 8822. [25] C. Lee, W. YangandR.G. Parr, Phys. Rev. B 37 (1988) 785. ]26] P.A.M. Dirac, Proc. Cambridge Phil. Sot. 26 ( 1930) 376: S.H. Vosko, L. Wilk and M. Nusair, Can. J. Phys. 58 (1980) 1200. [27] S. Huzinaga, J. Chem. Phys. 42 (1965) 293. [28] J. AImlGf and P.R. Taylor, Advan. Quantum Chem. 22 (1991) 301. [29] H.A. Carlson, T.B. Nguyen, M. Orozco and W.L. Jorgensen, J. Comput. Chem. 14 (1993) 1240. [30] M. Orozco, W.L. Jorgensen and EJ. Luque, J. Comput. Chem. 14 (1993) 1498. [31] I. Tut%n, E. Silla and J. Tomasi, J. Phys. Chem. 96 ( 1992) 9043; I. Alkorta. H.O. Villar and J.J. Perez, J. Comput. Chem. 14 ( 1993) 620, and references therein. [32] H. Hoshi, M. Sakurai, Y. Inowa and R. Chujo, J. Chem. Phys. 87 (1987) 1107; R. Cammi and J. Tomasi, J. Comput. Chem., in press. [33] R. Cammi and J. Tomasi, J. Chem. Phys. 100 ( 1994) 7495; 101 (1994) 3888.