Structure and Reactivity in Aqueous Solution - American Chemical

Chapter 4

Solvation Free Energies from a Combined Quantum Mechanical and Continuum Dielectric Approach

Downloaded by UNIV OF PITTSBURGH on February 2, 2016 | http://pubs.acs.org Publication Date: September 29, 1994 | doi: 10.1021/bk-1994-0568.ch004

Carmay Lim, Shek Ling Chan, and Philip Tole Protein Engineering Network of Centres of Excellence, Departments of Molecular and Medical Genetics, Biochemistry, and Chemistry, MSB 4388, University of Toronto, 1 King's College Circle, Toronto, Ontario M5S 1A8, Canada

We propose a new method to compute solvation free energies by combining quantum mechanical and continuum dielectric meth ods. Instead of using atom-based partial charges and a dielec tric boundary derived from atomic radii to solve the Poisson or Poisson-Boltzmann equation, the solute charge distribution and the dielectric boundary are obtained from an electronic wavefunction determined using quantum mechanical methods. This results in several advantages over conventional applications of the contin uum dielectric method as well as molecular dynamics and Monte Carlo simulation techniques. First, the input charge distribution is treated more accurately than that in the conventional applica tion which uses atomic partial charges determined by non-unique procedures. Second, the electronic solute wavefunction provides a better description of the solute shape and thus, reduces the error in determining the dielectric boundary compared to a geometric dielectric boundary defined by arbitrary effective Born radii of the solute atoms. Third, the new method is generally a cost-effective means of obtaining relative free energies compared to free energy perturbation methods, which require parameterization and inten sive cpu-time.

During the past decade, a large variety of methods have been developed to treat solvent effects on the static and dynamic properties of molecules. These meth ods fall generally into two classes depending on their treatment of the solvent molecules; viz., "microscopic" methods that treat solvent molecules and their interactions with the solute explicitly, and "macroscopic" methods that treat sol vent as a continuous medium characterized by a dielectric constant. The "macro0097^156W/0568^50$08.00/0 © 1994 American Chemical Society

In Structure and Reactivity in Aqueous Solution; Cramer, Christopher J., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.

4. LIM ET AL.

Solvation Free EnergiesfromCombined Approach

51

scopic" methods have a distinct advantage over the "microscopic" methods, espe cially for large systems, in speeding up computation considerably by not taking solvent molecules into account explicitly.


Under the class of "microscopic" methods, one of the earliest approaches was to employ quantum-mechanical methods to study the solute and a limited number of solvent molecules (1). This approach does not yield the energetics associated with the long-range bulk solvent forces and simply increasing the number of solvent molecules becomes computationally impractical. A breakthrough for modeling solution chemistry came with the development of condensed-phase molecular dynamics and Monte Carlo simulation procedures us ing empirical energy potentials. The procedure (2) for modeling solution reactions involves (i) determination of the minimum energy reaction path in vacuum via ab initio calculations (#), (ii) determination of the solute-solvent non-bonded pa rameters as a function of the reaction coordinate by fitting to ab initio results of the reacting complex with a water molecule in various orientations, (iii) deter mination of the corresponding potential of mean force in solution from a series of fluid simulations using free energy statistical perturbation or thermodynamic integration techniques (4). Recently, condensed-phase molecular dynamics and Monte Carlo methodology based on combined quantum and molecular mechanical potentials have been de veloped to simulate reactions in solution directly. Various combinations have been made: ab initio/Car-Parrinello methods with molecular dynamics (5) and semiempirical techniques with molecular dynamics (6) or Monte Carlo (7). These various approaches have in common that the system is partitioned into quan tum mechanical and molecular mechanical regions and the combined quan tum/molecular mechanical potential and forces are calculated at each step for use in classical dynamics simulations. A n approach bridging "microscopic" and "macroscopic" methods has been pro posed whereby a quantum mechanical solute is surrounded with a grid of point dipoles representing the average solvent polarization, and solvent outside this re gion is treated as a dielectric continuum (8). This approach has been further simplified by incorporating continuum solvent models (e.g., a generalized Bora term (9) or a sum of surface element terms (10)) in the solute Hamiltonian and carrying out self-consistent field calculations. The continuum solvent model can also be employed on its own to compute solvation free energies by solving a finitedifference Poisson equation (11,12); it has been used to estimate solvent effects on gas-phase activation free energy profiles (IS, 14). Unfortunately, the "microscopic" and "macroscopic" methods described above suffer from the disadvantage that extensive parameterization is required; e.g., pa rameters for van der Waals and electrostatic solute-solvent interactions in simula tions using empirical energy potentials (2), parameters for van der Waals solutesolvent interactions in simulations based on combined quantum and molecular mechanics potentials (5), parameters for the effective Born radius and accessible In Structure and Reactivity in Aqueous Solution; Cramer, Christopher J., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.


52

STRUCTURE AND REACTIVITY IN AQUEOUS SOLUTION

surface tension in quantum mechanical calculations employing a semiempirical Hamiltonian augmented by a generalized Bora term (9). However, experimental data are not always available to calibrate the parameters; this is often the case for short-lived transition states or reaction intermediates. Furthermore, parameteri zation against ab initio calculations on the reacting complex with a single water molecule may not yield sufficiently accurate parameters. "Microscopic" meth ods have the further disadvantage that free energy perturbation calculations are computationally demanding as simulations must run long enough to permit ade quate sampling of the relevant configuration space. If there are multiple minima with significant intervening barriers, the problem may be difficult to detect and overcome. Moreover, a simulation with combined quantum mechanical/molecular mechanical potentials takes two to three times longer than a corresponding sim ulation using empirical energy functions. Here, a new method for computing solvation free energies is proposed. The new procedure combines quantum mechanical and continuum dielectric methods by employing the quantum mechanical charge density directly in solving the finitedifference Poisson or Poisson-Boltzmann equation for the electrostatic contribu tion to the solvation free energy. To test the reliability of the new method, it was employed to calculate solvation free energies of monatomic ions and polyatomic molecules for which experimental data are available for comparison. These test calculations serve to calibrate (i) the basis set and level of quantum mechanical calculation required to yield the most accurate charge density, and (ii) the cut-off value of the charge density for determining the dielectric boundary. The method is outlined in the next section. This is followed by a discussion of the results of monatomic ions and polyatomic molecules. The performance of the method and its potential application to model reactions in solution are discussed in the last section. Method The continuum dielectric method sets up a grid in space and solves for the electric potential at each grid point given a molecule residing in a dielectric medium. In the standard continuum dielectric model (11, 12) the charge distribution on the grid is determined from partial atomic charges, and the dielectric boundary is derived from atomic radii of the solute atoms. The essence of the new procedure is to employ directly a charge density grid generated from quantum mechanical calculations and to derive the dielectric boundary from this electron distribution. The grid of electron charge densities was generated from the solute wavefunction using the ab initio program GAUSSIAN 92 (15). For monatomic ions, the basis sets employed were 6-31G*/6-31+G* for first and second row elements and ST0-3G* for third row elements. Unless stated otherwise, the geometries for polyatomic ions were fully optimized at the Hartree Fock (HF) level with the 6-31G* basis, and the charge density grids were computed at the second-order M0ller-Plesset perturbation (MP2) level. The grid for the continuum dielectric calculation was checked to ensure that the dielectric boundary was at least 1.5 Â from the grid edge. Unless stated In Structure and Reactivity in Aqueous Solution; Cramer, Christopher J., et al.; ACS Symposium Series; American Chemical Society: Washington, DC, 1994.


4. LIM ET AL.


53

otherwise, the grid-spacing (or grid-eye) was 0 . 2 Â x 0 . 2 Â x 0 . 2 Â ; it was the same as that used in the ab initio calculations so that the electron density output of GAUSSIAN 92 could be fed directly into the input of the continuum dielectric program. As GAUSSIAN 92 gave the electron densities at the grid points, when the density at each grid point was multiplied by the grid-eye volume and subsequently summed, the result might deviate from the total electron charge by as much as 10%. To correct this, the nuclear charges were scaled down so as to preserve the net charge of the molecule/ion. For example, if - Q is the total electron charge of the molecule/ion, - q is the sum of the electron charge on the grid, and Ν is the total nuclear charge, then all nuclear charges are scaled by a factor of (N-Q-f q)/N; this preserves the ratio of the charges on the nuclei. The scaled nuclear charges were then distributed between the eight nearest grid points by the procedure of Wachspress (16). The physical boundary of the molecule was defined using the electron density distribution: all grid points with an electron density above a certain cutoff value were considered to be inside the molecule. The low dielectric region was then defined as the space inaccessible to contact by a 1.4 Â solvent sphere rolling over the physical boundary. Thus, the definition of the dielectric boundary is similar to that in the conventional application of the continuum dielectric method. The difference lies in the definition of the physical boundary of the molecule, which is determined by a set of atomic radii in the conventional application. Here, the need to calibrate the atomic radii for different atom types is bypassed and replaced by a single cutoff value for the determination of the physical molecular boundary. The low dielectric region was assigned a dielectric constant of 1, whereas the high dielectric region was characterized by a dielectric constant of 1 for vacuum and 80 for water. The numerical procedure for solving the electric potential on the grid is the same as that in conventional applications. The electric potential φ at every grid point was initialized using Coulomb's Law with a relative permittivity of the solvent medium. Then Poisson equation was solved for the values of φ on the grid by an over-relaxation algorithm (17). The work W of assembling the electric charges is given by

th

where ç is the charge associated with the i grid point and φι is the electric po tential at that point; the summation goes over all the grid points. The difference in the work calculated in vacuum and solution yields the electrostatic part of the solvation free energy. t

Since each solvation free energy calculation required the difference between the work term in vacuum and solution, it took two iterative cycles for one free energy calculation. For each iterative cycle, 900 steps were performed; the maximum


54


Table L Effect of Diffuse Functions on the Electrostatic Solvation IVee E n e r g y Solvn free energy (kcal/mol) Charge Density Basis 6-3lG* e/Â 6-31+G** 106 Na+ 51x51x51 0.0002 101 Na+ 106 100 0.0002 51x51x51 CH NHJ 68.7 62.6 71x71x71 0.0005 ImH+ 78.0 47x97x97 0.0001 75.8 88.8 74.6 51x51x51 0.0025 (OH)81.2 94.0 (OH)51x51x51 0.0045 86.2 CH3COO- 61x61x61 0.0045 81.7 CH NH 8.2 11.0 61x61x61 0.0015 CH NH 6.1 71x71x71 0.0005 7.9 3.2 CH3COOH 61x61x61 0.0015 2.5 T h e charge densities were obtained from a single point MP2 calculation corresponding to HF/6-31G* geometries. The charge densities were computed at the H F level. Species

Grid Size

Cut-off

a

3

6

6


3

3

2

3

2

a

6

change in the magnitude of the potential at a grid point in the last step dropped to the order of several micro-volts. The cpu time required for each iterative cycle is strongly dependent on the linear dimension of the grid; if the linear dimension of the grid is doubled, the number of points in the grid is eight times the original. On an I B M 340, each iterative cycle took about 2 to 3 hours for a 5 0 Â x 5 0 Â x 5 0 Â grid and 5 to 6 hours for a 6 0 Â x 6 0 Â x 6 0 Â grid. R e s u l t s and Discussion Dependence on Charge Densities. The effect of diffuse functions in evaluat ing the charge densities on the electrostatic solvation free energy is summarized in Table I. Table I shows that for charged species, the electrostatic solvation free en ergies computed with MP2/6-31+G* charge densities are generally smaller than those calculated with MP2/6-31G* charge densities. This is because the 6-31+G* electronic charge densities of the ions are more spread out than the 6-31G* values, hence the same charge density cutoff defines a physical boundary that is further away from the nucleus; i.e., a larger effective Born radius, and thus a smaller solvation free energy. Dependence on Charge D e n s i t y Cut-off. The dependence of the electro static solvation free energy on the cut-off value employed for the charge densi ties is summarized in Table II for monatomic ions and Table III for polyatomic molecules. The computed and experimental solvation free energies correspond


4. LIM ET AL.


55

Table II. Effect of Charge Density Cut-off on the Electrostatic Solvation free Energy: Monatomics Species

Basis'*

Solvation free energy (kcal/mol) Charge density cut-off ( e / À ) O.52~ 5" IP" expt _c -(115-124) Li+ 6-31+G* -121 -136 -152 _c Be -563 6-31+G* -607 -721 -780 _c Na+ -(90-98) 6-31+G* -92 -101 -108 _c -443 Mg 6-31+G* ^121 -476 -511 _ C Na+ -(90-98) ST0-3G* -82 -92 -101 _c -443 Mg ST0-3G* -353 -400 -445 _c K+ -(74-75) -72 ST0-3G* -79 -65 _c Ca + ST0-3G* -363 -284 -327 -362 T h e charge densities were obtained from a single point M P 2 calculation corresponding to HF/6-31G* geometries; a 50 Â x 50 Â x 50 Â grid was used for the dielectric calculations. ^Experimental values are taken from Latimer W . M . ; Pitzer, K . S.; Slansky, C. M . ; J. Chem. Phys. 1939, 7, 108; R. Gomer, G . Tryson J. Chem. Phys. 1977, 66, 4413; H . L . Friedman and C. V . Krishnan Water: A Comprehensive Treatise, Eds: F . Franks, Plenum Press, N Y 1973, 3, 1. Value has not been computed for this cut-off. 3

4

4

4

4

6

2 +


2 +

2 +

2

e

c

to a common standard state of 1 M for both gas and solution. As the charge density cut-off increases, the computed solvation free energy increases. This can be rationalized using the Born formula for the heat of hydration W (18),

2r V

t) 0

which shows that the solvation free energy of a monatomic ion with charge q depends inversely upon the Born radius r. Thus, as the cutoff value of the electron density is increased, the volume defined as inside the ion decreases; i.e., the effective Bora radius decreases, and so the solvation free energy increases. In general, the same variation in cutoff has a smaller effect on the solvation free energy of a polyatomic ion than a monatomic ion. This is because the volume of a polyatomic ion is generally bigger than that of a monatomic ion and hence, for the same absolute difference in the positioning of the physical boundary, the effect on the overall volume mass is smaller for polyatomic ions. Comparison with experiment. Although experimental solvation free energies of neutral species may be directly measured, "experimental'' values of charged




Table III. Effect of Charge Density Cut-off on the Electrostatic Solvation Free Energy: Polyatomics Species

Basis*

Grid Size

Solvation free energy (kcal/mol) Cut-off R ( e / A ) R expt Â RJ -104 H 0+ 6-31G* -109 -104 61 -73 PH+ -72 -77 3-21G* 61 -70 CH NHj -72 -69 6-31G* 61 -(96-106) -94 6-31G* -101 (OH)51 -(77-82) -86 CH COO6-31G* -98 61 -76 (SH)6-31G* -78 -87 61 -6.7 C H C O O H 6-31G* -4.1 -2.3 61 -6.3 H 0 6-31G* 51 -10.3 -6.4 -4.6 CH NH -10.6 -6.1 6-31G* 61 -0.7 (SH) -3.5 6-31G* 51 -6.5 "The charge densities were obtained from a single point MP2 calculation corresponding to HF/6-31G* geometries. R i equals 0.0025, 0.0100, and 0.0030 for cations, anions and neutrals, respectively. R equals 0.0015, 0.0045, and 0.0005 for cations, anions and neutrals, respectively. Experimental values are taken from R. G . Pearson / . Am. Chem. Soc. 1986, 108, 6109. A . Ben-Nairn, Y . J. Marcus J. Chem. Phys. 1984, 81, 2016. 3

3

c

d

e

3

3

3

3

3

3

3

3

3

3

3

3

2

3

3

2

3

2

6

c

e

d


4. LIM ET AL.

Solvation Free Energies from Combined Approach

57

species are commonly determined indirectly from the experimental values of the solvation free energy AG of the corresponding neutral species and the proton, the gas-phase proton affinity or basicity ( A G ) and solution p K or free energy (AGgin) using Scheme I: S

g M

AG ->

a

g a 8

BH(g)

B(g)

-AG (BH) j Downloaded by UNIV OF PITTSBURGH on February 2, 2016 | http://pubs.acs.org Publication Date: September 29, 1994 | doi: 10.1021/bk-1994-0568.ch004

-» AG l s

B(s)

H+(g) +

I -AG.(B)

g

BH(s)

+

j -AG (H ) 8

+

H+(s)

(Scheme I)

n

Since experimental gas basicities are accurate to only 2-5 kcal/mol, and there are at least five different experimental values for the proton solvation free energy -(254 to 261) kcal/mole with a spread of 7 kcal/mol, the uncertainty in the "experimental" solvation free energies of charged species is at least ± 5 kcal/mol. On the other hand, the computed solvation free energies are subject to errors due to uncertainties in the charge densities and cut-off employed and the neglect of the non-polar contribution to the solvation free energy. 3

For the monatomic ions in Table I I , a 0.0010 e / Â cutoff of the M P 2 / S T O 3G* charge densities gives solvation free energies in excellent agreement with experiment, whereas a 0.5xl0~ e / Â cutoff of the MP2/6-31+G* charge densities yields solvation free energies that are within 7% of the experimental values. The STO-3G* basis requires a larger cut-off than the 6-31+G* basis as the latter results in a larger solvation free energy than the former for the same cutoff (Table II). For the polyatomic molecules in Table I I I , a charge density cut-off of 0.0015 e / Â for cations, 0.0045 e / Â for anions, and 0.0005 e / Â for neutrals, yields solvation free energies in accord with experiment. As expected, the agreement with experiment for neutral species is not as good as the charged ions since for neutral molecules, the non-polar contribution (around 1-2 kcal/mol) is no longer negligible relative to the electrostatic component of the solvation free energy. 4

3

3

3

3

Concluding Remarks A p p l i c a t i o n s . The present method can be applied to study reactions in aqueous solutions by computing the solvation free energies of reaction complexes. Thus, if X and Y are stationary points on the gas-phase ab initio potential energy surface and Δ Α | is the free energy difference from ab initio calculations, the present method can be used to obtain the solvation free energies of X and Y and thus, μ

Δ Α ^ from Scheme I I .


58


ΔΑ, X(gas)

Y(gas) |-ΔΑ,(Υ)

-ΔΑ,(Χ) I X(sln)

(Scheme II)

Y{sln)


t ΔΑ,'sin By piecing together consecutive Δ Α ^ , a qualitative free energy profile is ob tained in solution. Note that for methods which require parameterization, there are generally no experimental data available on short-lived reaction intermedi ates/transition states to calibrate parameters. The above procedure automatically yields the separate gas-phase and solvation contributions to the observed solution activation free energy. If the X and Y species involved have gas-phase activation and solvation free energy errors that are similar and cancel, the proposed method may yield reasonably accurate solvation and solution activation free energy differences. In such cases, the outcome of two competing reaction pathways in solution and substituent effects may be predicted. The present method can readily be applied to study reactions in organic solvents as well as the ionic-strength and temperature dependence of reaction rates. Advantages. The present methodology is not only versatile (as described above), it also possesses several advantages over existing methods, in principle. First, since the charge density is directly employed in the calculations, the need to fit solute partial charges empirically and errors due to ambiguous ways of extract ing partial charges (e.g., Mulliken population analysis) are eliminated. Second, as the dielectric boundary is determined by a single charge density cut-off value, the need to fit the Born radii empirically (0, 11) for various atom types (e.g., C, Η, 0, Ν, P, S) and errors due to arbitrary ways of defining the solute Born radii (e.g., in terms of van der Waal's, ionic or covalent radii) are eliminated. Third, if the number of grid-points considered is not excessive, the new method is a cost-effective means of obtaining relative free energies compared to expensive free energy perturbation methods. Note that the limitations in the quantum me chanical calculations; e.g., the basis set size and C P U requirements, are inherent in all existing methods. Potential Improvements Accuracy may be improved if the charge density out side the physical boundary of the molecule is set to zero and the nuclear charge rescaled accordingly, as described in the Method section. Another interesting avenue is to explore means of using the ab initio electrostatic potential grid di rectly instead of the quantum mechanical charge density distribution in solving for the electrostatic potentials in vacuum and in solution; however, the charge density grid will still be needed to compute the electrostatic self-energy (see equa tion in Method section). This may speed up convergence as well as potentially


4. LIM ET AL.


59

decrease the dependence of the electrostatic self-energy on the accuracy of the charge density. Acknowledgments This work was supported by the Protein Engineering Network of Centres of Ex cellence of Canada.


Literature Cited (1) Goldblum, Α.; Perahia, C.; Pullman, A. FEBS Lett. 1979, 91, 213. (2) Madura, J.; Jorgensen, W. L. J. Am. Chem. Soc. 1986, 108, 2517-2527. (3) Hehre, W. J.; Radom, L.; Schleyer, P. v.R.; Pople, J. A. Ab Initio Molecular Orbital Theory. John Wiley and Sons, (1986). (4) Beveridge, D. L.; DiCapua, F. M. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 431-492. (5) Field, M. J. J. Phys. Chem. 1991, 95, 5104-5108. (6) Field, M.J.;Bash, P. Α.; Karplus, M. J. Comp. Chem. 1990, 11, 700-733. (7) Gao, J. J. Phys. Chem. 1992, 96, 537-540. (8) Luzhkov, V.; Warshel, A. J. Comp. Chem. 1992, 13, 199-213. (9) Cramer, C. J.; Truhlar, D. G. J. Computer-Aided Molecular Design 1992, 6, 629-666. (10) Ford, G. P.; Wang, B. J. Am. Chem. Soc. 1992, 114, 10563-10569. (11) Lim,C.;Bashford, D.; Karplus, M. J. Phys. Chem. 1991, 95, 5610-5620. (12) Gilson, M . K.; Honig, Β. H. Proteins: Structure, Function and Genetics 1988, 4, 7-18. (13) Lim,C.;Tole, P. J. Phys. Chem. 1992,114,7245. (14) Tole, P.; Lim, C. J. Phys. Chem. 1993, 97, 6212. (15) Frish, M. J.; Head-Gordon, M.; Trucks, G. W.; Foresman, J. B.; Schlegel, Η. B.; Raghavachari, K.; Robb, M.; Binkley, J. S.; Gonzalez, C.; Defrees, D. J.; Fox, D.J.;Whiteside, R. Α.; Seeger, R.; Melius, C. F.; Baker, J.; Martin, R. L.; Kahn, L. R.; Stewart, J. J. P.; Topiol, S.; Pople, J. A. Gaussian Inc.; Pittsburgh, PA 15213, USA. (16) Wachspress, Ε. I. Iterative Solution of Elliptic Systems, Prentice Hall, Englewood Cliffs, (1966). (17) Press, W. H.; Flannery, B. P.; Teukolsky, S. Α.; Vetterling, W. T. Numerical Recipes. The Art of Scientific Computing; Cambridge University Press, Cam bridge, (1986). (18) Born, Μ. Ζ. Z. Phys. 1920, 1, 45. (19) Chan, S. L.; Lim, C. J. Phys. Chem. 1994, 98, 692. RECEIVED April 25, 1994


Structure and Reactivity in Aqueous Solution - American Chemical

Recommend Documents