Benchmarking the Conductor-like Polarizable Continuum model

Benchmarking the Conductor-like Polarizable Continuum model

(Parte 1 de 3)

Benchmarking the Conductor-like Polarizable Continuum

Model (CPCM) for Aqueous Solvation Free Energies of Neutral and Ionic Organic Molecules

Yu Takano and K. N. Houk*

Department of Chemistry and Biochemistry, UniVersity of California, Los Angeles, 607 Charles E. Young DriVe East, Los Angeles, California 90095-1569

Received June 17, 2004

Abstract: The conductor-likepolarizablecontinuummodel (CPCM) using several cavity models is applied to compute aqueous solvation free energies for a number of organic molecules (30 neutral molecules, 21 anions, and 19 cations). The calculated solvation free energies are compared to the available experimental data from the viewpoint of cavity models, computational methods, calculation time, and aqueous pKa values. The HF/6-31+G(d)//HF/6-31+G(d) and the HF/6-31+G(d)//B3LYP/6-31+G(d) with the UAKS cavities, in which radii are optimized with

PBE0/6-31G(d),provideaqueoussolvationeffectsin best agreementwith availableexperimental data. The mean absolute deviations from experiment are 2.6 kcal/mol. The MP2/6-31++G- (d,p)//HF/6-31+G(d) with the CPCM-UAKS(HF/6-31+G(d)) calculation is also performed for the base-catalyzed hydrolysis of methyl acetate in water.


Many chemical and biological reactions occur in water, where the polar and ionic processesare much more favorable than in the gas phase. Many efforts have been devoted to the development of methods to compute reaction barriers and energetics occurring in condensed phases with experimental accuracy.1 Effective explicit water models become available for the description of chemical systems in liquid solution.1 However,with high-levelquantummechanics,only a limited number of solvent molecules can be included explicitly due to the high cost of the calculations.

The goal of this work is to determine which theoretical procedureprovidesthe most quantitativeestimateof aqueous solvation effects, so that the rates of chemical and biological reactions in water can be computed accurately. One of the most successful solvation models is the conductor-like polarizablecontinuum model (CPCM).2 Here we benchmark different variations of CPCM for the computation of solvation energies of neutral and ionic organic species and compare them to several other works. The CPCM method has also been applied to the computation of the alkaline hydrolysis of methyl acetate in aqueous solution.


Dielectric continuum theories1 are now widely used to describe hydration in conjunction with quantum mechanical calculations due to the relatively low cost of the calculation. CPCM2 and PCM3 are two of many successful solvation models. In their approaches, the solute interacts with the solventrepresentedby a dielectriccontinuummodel. The solute molecule is embedded into a cavity surrounded by a dielectric continuum of permittivity . The accuracy of continuum solvation models depends on several factors; the most important one is the use of proper boundary conditions on the surface of the cavity containing the solute. CPCM and PCM define the cavities as envelopes of spheres centered on atoms or atomic groups: a number of cavity models have been suggested. Inside the cavity the dielectric constant is the same as in vacuo, outside it takes the value of the desired solvent. Once the cavity has been defined, the surface is smoothly mapped by small regions, called tesserae. Each tessera is characterized by the position of its center, its area, and the electrostatic vector normal to the surface passing through its center. Recently, the CPCM method has been improvedand extendedin GAUSSIAN034a so that the cavity can be selected in a number of different ways.* Corresponding author e-mail:

10.1021/ct049977a C: $30.25 © 2005 American Chemical Society Published on Web 1/13/2004

In CPCM, the solvation free energy can be expressed1a

¢Gel is the electrostaticcomponentof ¢Gsolv. The Gel term is calculated using the CPCM self-consistent reaction field

(SCRF) method.2 The cavitation term, ¢Gcav, is calculated with the expression derived by Pierotti from the hard sphere theory5 and adapted to the case of nonspherical cavities.3b

The dispersion and repulsion terms, ¢Gdis and ¢Grep, are computed following Floris and Tomasi’s procedure,6 with the parameters proposed by Callet and Claverie.7 The qrot,g, qvib,g, qrot,s, and qvib,s are denoted the microscopic partition functions for rotational and vibrational states of the solute in the gas phase and in solution, respectively; nsolute,g and nsolute,s are the numeral densities of solute; and ⁄solute,g and

⁄solute,s are the momentum partition functions. The last term, P¢V, may be neglected since its value is normally less than

10-3 kcal/mol.8 The quantity, -RT ln(nsolute⁄solute), is a free energy correction to account for solute occupying the entire volume available in the reference state. For simple models such as isotropic solutions with no chemical association or dissociation processes, this contribution is equal to zero. The term involving the vibrational and rotational degrees of freedom, ln(qvib,g/qvib,s) and ln(qrot,g/qrot,s), is negligible.1a The last three terms in eq 1 are neglected in the PCM and

CPCM formulations.1a It, however,has been noted that things are actually more complex when one considers dimers or trimers held together by relatively weak interactions and for chemical association or dissociation processes.1a

Computational Method

The CPCM-SCRF calculations2 at the HF/6-31+G(d) and B3LYP/6-31+G(d) levels were carried out on the stationary points to address solvation effects. A dielectric constant of 78.39 was utilized in order to simulate aqueous environment. The CPCM calculations were performed with tesserae of 0.2 Å2 average size. All structures were optimized at the HF/6- 31+G(d) and B3LYP/6-31+G(d) levels9-1 in the gas phase and at the B3LYP/6-31+G(d) level in the aqueous environment. All stationary points were characterized by frequency calculations at the same level. All calculations were carried out with GAUSSIAN034a and GAUSSIAN98.4b We used 1 mol L-1 as the standard state for both the gas phase and the solution for all thermodynamic properties.

In CPCM and PCM, the choice of cavities is important because the computed energies and properties depend on the cavity size. In this study, the UA0, UAHF, UAKS, UFF, PAULING, and BONDI cavities were used to evaluate the aqueous solvation effects using CPCM and PCM. The UA0 cavity is built up using the united atom topological model (UATM)8a applied on atomic radii of the universalforce field (UFF).9 By default, the UA0 model is chosen to build the cavity in GAUSSIAN03.14 The UAHF and UAKS cavities use UATM with radii optimized for the HF/6-31G(d) and PBE0/6-31G(d)10 levels of theory, respectively. The UAHF model is the default cavity of GAUSSIAN98. A set of the radii from UFF is used for making the UFF cavities. For the PAULING and BONDI cavities, each solute atom and group is assigned van der Waals values obtained from Pauling12b or Bondi12c atomic radii.

Test of the Reliability of CPCM

Aqueous solvation free energies for a number of organic molecules (30 neutrals, 21 anions, and 19 cations) computed using the CPCM method.2 We investigated the calculated solvation free energies compared to the available experimental data16 from the viewpoint of cavities, computational methods, calculation time, and aqueous pKa values. The computedaqueoussolvationeffectswere also comparedwith solvation energies computed using COSMO,2 clustercontinuummodel,17 SM5.42R,18 PCM,3 andIPCMmethods.19

Cavity Models. Table 1 summarizes the mean absolute deviations (MADs) of the aqueous solvation free energies calculatedwith seven cavitiesat the HF/6-31+G(d)//B3LYP/ 6-31+G(d) level9-1 from the experimental data16 for 70 neutral and charged molecules.

CPCM with the new cavities, UAKS and UAHF(G03), has improved accuracy in the aqueous solvation energies for a set of 70 organic species, though these methods still do not achieve the accuracy of experimental data (0.2 kcal/mol for neutral molecules20,21 and 2 kcal/mol for ions16). The MADs calculated by the CPCM-UAKS and CPCM-UAHF- (G03) are 2.61 kcal/mol (1.35, 3.21, and 3.93 kcal/mol for 30 neutrals, 21 anions, and 19 cations, respectively)and 2.84 kcal/mol (1.10, 3.92, and 4.93 kcal/mol), respectively. On the other hand, the CPCM-UA0 and CPCM-UFF methods fail for charged molecules with MADs of 9.64 (13.09) and 9.30 (15.12) kcal/mol for anion (cation) solutes, respectively. The PAULING cavitiesshow the best solvationfree energies (2.73 kcal/mol) for anion molecules but give the worst agreement with experiment for neutrals. For all cavities but the PAULING cavities,the calculatedsolvationfree energies for the neutral moleculesare much closer to the experimental results than those for the charged species.The large solvation energy errors for the ions is due to inadequate treatment of specific short-range interactions, probably associated with strong hydrogen bonds between the ions and first-shell water molecules. Dielectric continuum theory1 cannot account for short-range solute-solvent interactions such as hydrogen bond. In addition, since anions and cations have aqueous solvation free energies in the range of 60-110 kcal/mol in

RT ln (qrot,gqvib,g qrot,sqvib,s

Table 1. Mean Absolute Deviations (MADs) of the Aqueous Solvation Free Energies of 70 Neutral and Charged Molecules at the HF/6-31+G(d)//B3LYP/ 6-31+G(d) Level Using CPCM with Several Cavitiesa total neutral anion cation a MADs are shown in kcal/mol.

Benchmarking the CPCM for Aqueous Solvation Free Energies J. Chem. Theory Comput., Vol. 1, No. 1, 2005 71

contrast to neutral molecules (0-10 kcal/mol), it is hard to achieve the 1 kcal/mol level of accuracy in prediction of the solvation free energy of ions.

Computational Methods. Table 2 shows the computational method dependence of MADs of the solvation free energy for the neutral and charged solutes. The geometries were optimizedin the gas phase,and the UAKS cavitieswere used. The MADs of the solvation free energies at the HF/ 6-31+G(d)//B3LYP/6-31+G(d) level become 0.7 kcal/mol smaller than those at the B3LYP/6-31+G(d)//B3LYP/6- 31+G(d) level. Especially, MADs of the solvation energies for the anion solutes are improved by 2.43 kcal/mol. There is a tendency that the cavity that reduces MADs of charged species increasesthe MADs of the neutral species.Optimization at the HF level provides MADs similar to those at the B3LYP level, indicating that optimized geometries using HF and B3LYP are close to each other.

Calculated MADs of the neutral and charged solutes with the geometries optimized in water are shown in Table 3. The geometries optimized in water make the MADs of charged species small while those for neutrals large. As a whole, optimizationin water shows MADs similarto those in vacuo, implying that the optimized structures in water are similar to those in vacuo. In Table 4, the optimized geometrical parameters for some examples in vacuo and water are reported. It is apparent that the effect of reoptimization in water on the geometrical parameters are very small. However, some charged species failed to optimize geometries in water due to the dissociation of a proton, especially using the UAKS, UAHF, and PAULING cavities (Tables S10- S15).

The 6-31G, 6-31G(d), 6-31+G(d), 6-31+G(d,p), and 6-311+G(2d,p) basis sets were utilized to investigate the dependenceon basis sets for aqueous solvation free energies. Even when basis sets are enlarged up to 6-311+G(2d,p), the MAD values from the experiment are very similar as shown in Table 5, indicating that diffuse and polarization functions of basis sets hardly change the aqueous solvation effects.

Calculation Time. Figure 1A shows the CPU time required to calculate the hydration energy using the previous (GAUSSIAN98)4b and the present (GAUSSIAN03)4a versions of the code. The present version provides a remarkable decrease of computational time because of the introduction of the fast multipolemethod to computethe solvationcharge.

In Figure 1B, we compare CPCM to PCM with respect to the CPU time for the estimate of the aqueous solvation free energies in water. Although PCM is faster than CPCM for smallmolecules,CPCM is fasterwhen the moleculesbecome larger. In the CPCM approach, the electrostatic problem related to solute-solvent interaction can be solved with a

Table 2. Mean Absolute Deviations (MADs) of the Aqueous Solvation Free Energies of 70 Neutral and Charged Molecules at the HF/6-31+G(d)//B3LYP/ 6-31+G(d), B3LYP/6-31+G(d)//B3LYP/6-31+G(d), and HF/6-31+G(d)//HF/6-31+G(d) Levels Using CPCM with the UAKS Cavitiesa

HF//HF HF//B3LYP B3LYP//B3LYP a MADs are shown in kcal/mol.

Table 3. Mean Absolute Deviations (MADs) of the Aqueous Solvation Free Energies of 70 Neutral and Charged Molecules at the HF/6-31+G(d)//B3LYP/ 6-31+G(d) Level Using CPCM with the UAKS Cavitiesa,b water vacuo a MADs are shown in kcal/mol. b The geometries were optimized in the water environment and in vacuo.

Table 4. Optimized Geometrical Parameters and Aqueous Solvation Free Energies for Some Sample Molecules in Vacuo and in Water

CH3O- CH3OH CH3NH3+ vacuo water vacuo water vacuo water

a ¢Gsolvs are shown in kcal/mol. b Bond lengths are shown in Å. c Angles are shown in degrees.

Table 5. Mean Absolute Deviations (MADs) for the Aqueous Solvation Free Energies of 70 Neutral and Charged Molecules at the HF and B3LYP Levels Using CPCM-UAKS with Five Different Basis Setsa total neutral anion cation total neutral anion cation

a The geometries were optimized in vacuo at the B3LYP/6-31+G(d) level.

72 J. Chem. Theory Comput., Vol. 1, No. 1, 2005 Takano and Houk

much simpler formalism than in PCM. This simpler formalism can be a faster estimate of the aqueous solvation free energies when larger systems are studied.

Aqueous pKa. Aqueous pKa values were evaluated using the B3LYP/6-31+G(d)//B3LYP/6-31+G(d) level of theory and inclusion of solvent effects at the HF/6-31+G(d) and

B3LYP/6-31+G(d) levels. At present, many theoretical pKa predictionshave been reported.17b,2 The vast majority of pKa calculations use the direct definition shown in eq 2 with combination of the experimental data of the proton.1c,22a-f

Since this reaction, however, involves the formation of charged species starting from neutral molecules, the procedure using eq 2 yields large errors due to imbalance of the charged and neutral species and requires the value of the experimental solvation free energy of the H+ ion. By comparison, reactions that conserve the number of charged species are more suitablefor accuratecalculationsof changes in solvation free energies. Pliego and Riveros used the proton-transferreactionbetween the HA acids and hydroxide anion as shown in eq 3 along with the continuum-cluster method17a and compared their calculated values with the pKa values calculated with the other solvation models. Since our objective is the benchmarking of the CPCM method, we made use of the reported results by Pliego and Riveros for comparison of the CPCM-UAKS and performed the pKa calculations according to their procedure.17b,22g,h

(Parte 1 de 3)