Combining implicit solvation models with hybrid ... - Semantic Scholar

JOURNAL OF CHEMICAL PHYSICS

VOLUME 117, NUMBER 10

8 SEPTEMBER 2002

Combining implicit solvation models with hybrid quantum mechanicalÕ molecular mechanical methods: A critical test with glycine Qiang Cuia) Department of Chemistry and Theoretical Chemistry Institute, University of Wisconsin, Madison, Madison, Wisconsin 53706

共Received 13 May 2002; accepted 17 June 2002兲 A combined approach to study reactions in solution in which the solute and a number of solvent molecules are described with a hybrid quantum mechanical/molecular mechanical 共QM/MM兲 method, and the bulk solvent is represented by a polarizable continuum model 共PCM兲 has been implemented. In this way, both short-range effects of the first-solvation shell and long-range electrostatics due to the bulk solvent can be taken into account. By carefully choosing the size of the solute–solvent cluster and the QM/MM partition, the current QM/MM/PCM approach can offer both computational efficiency and accuracy. The approach has been illustrated by two simple systems: water-dimer and glycine in water. The results demonstrated that the current approach offers a satisfactory description of solvation effects on the geometry and energetics of neutral and charged hydrogen-bonding systems. The method correctly produced the relative stability of the zwitterionic and neutral forms of glycine in solution, which was found to be a subtle issue in previous studies. The approach can be extended to study reactions in biomolecules in which part of the system is treated with QM/MM, and the bulk solvent plus part of the protein or nucleic acids are described with either a continuum or approximate microscopic representation. © 2002 American Institute of Physics. 关DOI: 10.1063/1.1499481兴

I. INTRODUCTION

solvent is represented at an approximate but microscopic level, such as the Langevin dipole approach of Warshel and co-workers.15 Integral equation theories 16 共such as the reference interaction site model, RISM17兲 have also been used with both classical and quantum mechanical descriptions of the solute, although the applications have been mostly restricted to relatively small molecular systems.18 Density functional theory has been found to be more reliable in the presence of high-valent ions 共such as magnesium兲19 than NLPB, but its implementation has been so far limited to simple solute/solvent systems.20 With careful parametrizations 共such as the choice of atomic radii兲, implicit solvent models can generate a qualitatively correct description for the influence of solvent molecules on the structure and energetics of the solute.1,4 It has been recognized for a long time, however, that the first several solvation shells around the solute tend to give rise to effects that are difficult to describe with a uniform continuum representation. These include, to name a few, the nonlinear response of the solvent when the solute is highly charged, the dispersion and hydrogen bonding interactions between solute and solvent, and the hydrophobic effect.1共b兲 One straightforward solution to this problem is to use the ‘‘supermolecule’’ or the ‘‘solute–solvent cluster’’ approach in which a number of solvent molecules are included explicitly to represent the first several solvation shells and the rest of the solvent molecules are represented by a dielectric continuum model. It should be emphasized that even with the explicit treatment of the first solvation shell, the contribution due to the bulk solvent can still make a qualitative difference on the result if the solute is highly charged or if solute prop-

Solvent molecules play an important role in determining the energetics and rate constants of chemical reactions that occur in solution1 and macromolecules.2 Solvation effects are also critical to the structural and dynamical properties of macromolecules.3 Although a full microscopic description of solvation is possible with molecular dynamics or Monte Carlo simulations in which the solvent molecules are represented explicitly, it is often time-consuming to do so and much attention has been focused on developing reliable implicit solvent models in which the solvent molecules are replaced by a structureless dielectric continuum.1 In the context of macromolecular simulations, most implicit solvent models divide the solvation effect into nonpolar and electrostatic contributions.4 The former is often described with terms related to the amount of exposed surface area,5 and the latter is often treated with the Poisson–Boltzmann equation 共PB兲6 or relevant simplifications 共such as the generalized Born model兲.7 When the solute is highly charged, the effect of counterions in solution can become crucial, for which nonlinear Poisson–Boltzmann 共NLPB兲8 equations based on the Debye–Hu¨kle theory9 has been found to be qualitatively adequate for low-valent ions.10 In the arena of quantum chemistry, popular implicit solvent models are based on numerical solution of the Poisson equation,11,12 and the influence due to counterions can be introduced with NLPB. Other interactions between the solute and the solvent molecules such as Pauli-repulsion13 and dispersions14 have also been considered. We note that there are also models in which the a兲

Electronic mail: [email protected]

0021-9606/2002/117(10)/4720/9/$19.00

4720

© 2002 American Institute of Physics

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Combining implicit solvation models with QM/MM methods

4721

erties that are sensitive to the environment 共e.g., nuclear magnetic resonance chemical shift21 or absorption spectra22兲 are of interest. An interesting example that will be described in more detail in the current work concerns the relative stability of the neutral and zwitterionic forms of glycine in water.23,24 Without the explicit treatment of the first solvation shell, several popular continuum models gave the wrong trend 共i.e., the neutral form being more stable兲.23 Even with an explicit description of the first solvation shell with either ab initio or the effective fragment potential 共EFP兲 method,25 the correct trend cannot be reproduced unless a reliable implicit model is used;23共b兲,23共c兲 e.g., combining the EFP treatment with the Onsager model26 still left the neutral form being more stable.23共a兲 When the ‘‘solute’’ is a macromolecule such as an enzyme, hybrid quantum mechanical/molecular mechanical 共QM/MM兲 methods27 can be used to describe the solute and solvent molecules. Very often only a small number of explicit solvent molecules are included in these calculations 共e.g., in the popular stochastic boundary approach28兲, because such a setup allows an efficient configurational sampling of the protein atoms. The limitation of such a setup has been recently emphasized by several authors.29 For example, it was found that a satisfactory estimate of reduction potential of FAD in cholesterol oxidase could only be obtained if effects of the bulk solvent are taken into account in the QM/MM free energy simulations.29共a兲 In that study, the key aspect is that the partial charges of the charged residues in the unsolvated system 共enzyme plus a small number of explicit water兲 were scaled down based on PB calculations to mimic the solvent screening effect;30 this avoids undesired structural distortions and overpolarization of the QM region. The effects of charge scaling and bulk solvation were then estimated by another set of Poisson–Boltzmann calculations for selected configurations from the molecular dynamics simulations. In these PB calculations, the QM atoms were replaced by a set of Mulliken charges. It is clearly of interest to examine the reliability of such an approach by explicitly combining QM/MM methods with PB calculations. Along this line, the QM/MM method was recently combined with a generalized Born model by Field and co-workers.31 In the current work, we have implemented a combined approach that will be referred to as QM/MM/PCM, in which the hybrid QM/MM method is used to describe the solute and first several solvation shells with the bulk solvent treated with an implicit solvation model. In particular, the polarizable continuum model 共PCM兲11 implemented in GAMESS32 was used for the latter. In such a way, effects due to both the first solvation-shell and long-range electrostatics in the bulk can be taken into account. By choosing the size of the solute–solvent cluster and the QM/MM partition judiciously, such a multilayer method can offer both accuracy and speed to suit the problem at hand. Such a QM/MM/PCM approach is conceptually similar to the EFP/Onsager and EFP/PCM models of Gordon and co-workers,23 although the latter is somewhat restricted in practical applications by the structural rigidity of the effective fragments;25 the advantage of EFP is that polarization of the discrete solvents is included explicitly. The current method can also be considered as a special

case of the ONIOM-PCM approach of Morokuma and co-workers.33 The basic theory is briefly summarized in Sec. II, and several test calculations are presented in Sec. III. A few conclusions are summarized in Sec. IV.

II. THEORY AND IMPLEMENTATION

In this section, we first briefly review the theory of the polarizable continuum model 共PCM兲 used here, and then discuss the necessary modifications required in the extensions to the QM/MM framework. The current work used the matrix inversion implementation of PCM11共b兲 in GAMESS, although extension to the integral equation formulation of PCM11共c兲 will be straightforward. A. The polarizable continuum model

The solvation free energy of a solute molecule can be written as the sum of several contributions, ⌬G sol⫽⌬G el⫹⌬G rep⫹⌬G dis⫹⌬G cav ,

共1兲

where the first three terms on the right-hand side of Eq. 共1兲 represent the electrostatic, Pauli repulsion, and dispersion interactions, respectively, between the solute and the 共implicit兲 solvent molecules; the last term, ⌬G cav , represents the reversible work required to create a cavity of the molecular shape of the solute.34 In the following discussions, we focus on the treatment of the electrostatic term; procedures for the calculation of the other terms have been discussed in the literature.13,14 The electrostatic component tends to be the dominant contribution in solvation free energies especially when the solute is highly charged; this is particularly true for the supermolecular approach that will be discussed in Sec. II B, because the short-ranged Pauli repulsion and dispersion interactions between the solute and solvent are accounted for explicitly with QM/MM. To obtain the electrostatic solvation free energy at zero salt condition, the Poisson equation corresponding to the charge distribution of the solute immersed in a continuum medium has to be solved numerically. In the PCM model,11 the boundary element methodology was employed, in which the polarization effect in the continuum medium induced by the charge distribution of the solute ( ␥ (r៝ )) is represented by a set of ‘‘surface charges ( 兵 q i 其 )’’ on the solute–solvent boundary. The two sets of charge distributions are given in the following forms: N nuc

␥ (r៝ )⫽ 兺 Z j ␦ (R៝ j ⫺r៝ )⫹ 兺 P ␮ ␯ ␹ ␮ 共 r៝ 兲 ␹ ␯ 共 r៝ 兲 ,

共2兲

q i ⫽ ␴ 共 s៝ i 兲 a i ,

共3兲

j

␮␯

where Z i is the nuclear charge, P ␮ ␯ is the density matrix in the atomic orbital basis ( ␹ ␮ , ␹ ␯ ), ␴ (s៝ i ) is the polarized charge density on the ith solute–solvent boundary element 共‘‘tessera’’兲, characterized by vector of sជ i and area a i . With these notations, the total electrostatic free energy for the solute (G el) including the solute–solvent electrostatic interaction can be written as

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

4722

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Qiang Cui

G el⫽ 具 ⌿ sol兩 H 0 兩 ⌿ sol典 ⫹ 21 U ␥ ␴ ⫹V NN ⫽ 具 ⌿ sol兩 H 0 兩 ⌿ sol典 ⫹ 21 共 U ee⫹U eN⫹U Ne⫹U NN兲 ⫹V NN, 共4兲 where ⌿ sol is the solute wave function including the electrostatic polarization effect due to the continuum solvent, and the four U terms represent the electrostatic interaction between the solute electron, solute nuclei, and the surface charges induced by them. In the case of Hartree–Fock or density functional theory, the expressions can be made more explicit as follows: G el⫽ 具 ⌿ sol兩 H 0 兩 ⌿ sol典 ⫹ 21 共 tr PX⫹tr PJ⫹tr PY⫹U NN兲

⬘ , ⫹V NN⫽tr Ph⬘ ⫹ 21 tr P⌫ ⬘ ⫹V NN

共5兲

where the one- and two-electron matrices as well as the nuclear repulsion terms are modified to include the electrostatic interaction due to the implicit solvent, h⬘⫽h⫹ 21 共 J⫹Y兲 ,

共6a兲

⌫ ⬘⫽⌫⫹X,

共6b兲

⬘ ⫽V NN⫹ U NN , V NN

共6c兲

1 2

where J represents the interaction between nuclei-induced surface charge (q Ni )and solute electrons, Y represents the interaction between electron-induced surface charge 关 q ␮e ␯ (i) 兴 and solute nuclei, and X (U NN) is the interaction between electron 共nuclei兲-induced surface charge and the solute electron 共nuclei兲. The relevant expressions for J, X, Y, and U NN have been derived in details in Ref. 11. Here we only give the expressions for the one-electron terms 共J,Y兲 and U NN , which are of relevance to the later discussions on extensions to the QM/MM framework, N TS

J ␮␯⫽

兺i

N TS

V ␮i ␯ q Ni ⫽

兺i

N TS

Y ␮␯⫽

兺i

N TS

U NN⫽

冋

具 ␹ ␮兩

N TS

V Ni q ␮e ␯ 共 i 兲 ⫽ N TS

冋

1 兩 r៝ ⫺s៝ i 兩

N nuc

兺i 兺j

冋

N nuc

兺i V Ni q Ni ⫽ 兺i 兺j

册

兩 ␹ ␯ 典 q Ni ,

Zj

៝ ⫺s៝ 兩 兩R j i Zj

៝ ⫺s៝ 兩 兩R j i

册

册

q ␮e ␯ 共 i 兲 ,

q Ni ,

共7a兲

共7b兲

共7c兲

where N TS is the number of tesserae. The surface charges 关 q Ni and q ␮e ␯ (i)] can be determined with the matrix expression for the boundary element solution of the Poisson equation, DA⫺1q⫽⫺„Ene⫹EnN…⫽

冉

冊

⳵ Ve ⳵ VN ⫹ , ⳵n ⳵n

共8兲

where n is the norm vector of the tesserae and EneÕN, VeÕN are the electrostatic field and potentials at tesserae due to the solute electron and nuclei 关see Eq. 共7兲兴. The matrix D contains geometrical properties of the solute cavity and dielectric properties of the implicit solvent,11 and A is a diagonal matrix containing the areas of the tesserae. With Eqs. 共5兲–共8兲, the Hartree–Fock–Roothan 共Kohn– Sham兲 equations can be derived and solved self-consistently to obtain the solute wave function 共density兲, free energy, and molecular properties35 in solution.

To perform geometry optimization and molecular dynamics simulations, analytical gradients of the solvation free energy is required. This has been discussed in detail in Ref. 36, and only relevant formulas will be given here. The derivative of the electrostatic free energy 关Eq. 共5兲兴 with respect to a nuclear displacement, a, is given in the following form taking into account the fact that the PCM-SCF equations are variational,36,37

⬘a , G ael⫽tr Ph⬘ a ⫹ 21 tr P⌫ ⬘ a ⫺tr SaW⫹V NN

共9兲

where W is the eigenvalue-weighted density matrix W ⫽PF⬘P, and h⬘a, ⌫⬘a, and Sa are the integral derivatives of h⬘, ⌫⬘ and S with respect to a. The expressions for the integral derivative of the PCM related one-electron quantities 共J,Y兲, and U NN can be derived with chain rules in a rather straightforward manner,36

冋冉 冊 兺 再冉 冊 兺 冋冉 冊

N TS X

J ␮A␯ ⫽

兺i

N TS X

Y ␮A␯ ⫽

i

N TS

X

A ⫽ U NN

i

冉 冊册 冋 册冎 冉 冊册

⳵ V ␮i ␯ N ⳵ q Ni q i ⫹V ␮i ␯ ⳵XA ⳵XA

⳵ V Ni e ⳵ q ␮e ␯ 共 i 兲 q ␮ ␯ 共 i 兲 ⫹V Ni ⳵XA ⳵XA ⳵ V Ni N ⳵ q Ni q i ⫹V Ni ⳵XA ⳵XA

共10a兲

,

共10b兲

,

共10c兲

.

While the derivatives for the electrostatic potentials are straightforward, the derivatives for the q’s have to be obtained by taking derivatives of Eq. 共8兲,

冉 冊冉 冊 冉 冊冉 冊 冉 冊

冉 冊

⳵ ENn ⳵ qN ⳵ A ⫺1 N ⳵ D⫺1 N ⫽⫺AD⫺1 ⫺ D En ⫺A En , ⳵XA ⳵XA ⳵XA ⳵XA 共11a兲 ⳵ q␮e ␯ ⳵ E␮e ␯ ⳵ A ⫺1 e ⫽⫺AD⫺1 ⫺ D E␮ ␯ ⳵XA ⳵XA ⳵XA ⫺A

⳵ D⫺1 e E␮ ␯ . ⳵XA

共11b兲

The expressions for the derivatives of cavity-related quantities such as D⫺1 , A and the vectors defining the tessera and its norm 共n兲 are given in Ref. 36. B. Extension to the hybrid QMÕMM framework

In the QM/MM/PCM method, the solute and a number of solvent molecules 共‘‘solute–solvent cluster’’兲 are described with QM/MM and the bulk solvent is treated with PCM. A typical application will describe the solute with QM, the first several solvation shells with MM solvents, and the rest with a dielectric continuum. In such a way, both the firstsolvation shell effects 共such as short-range Pauli repulsions, dispersion interactions, and solvent entropic terms related to hydrophobic interactions兲 and long-range electrostatic effects in the bulk can be taken into account with modest computational cost. The small number of explicit solvent molecules makes the calculations of quantities that require configurational samplings much faster compared to a full-scale microscopic simulation 共see Sec. III B兲. For negatively charged solute species with diffuse charge distributions, the QM/

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Combining implicit solvation models with QM/MM methods

MM/PCM model also serves as a systematic way to extend the cavity size 共determined by the space occupied by the MM atoms兲 to avoid the charge-escape problem that causes the violation of Gauss’s theorem.38 For enzymatic applications, the active site is described with QM, the rest of the system plus a limited number of explicit solvent molecules are described with MM, and the bulk solvent is described with an implicit solvent model. The MM atoms in the current implementation are represented by fixed partial charges, and they interact with the QM atoms through electrostatic and van der Waals terms. With such a scheme, the total electrostatic free energy of the system in solution 关Eq. 共5兲兴 is modified to be

where hQM/MM is the direct electrostatic interaction between the QM electrons and MM atoms in the familiar form

G el⫽tr P关 h⬘ ⫹hQM/MM⫹ 21 共 JQM/MM⫹YQM/MM兲兴 QM/MM 1 QM/MM ⬘ ⫹E MM ⫹ 12 tr P⌫ ⬘ ⫹ 共 V NN ⫹ 2 U NN el ⫹V NN MM/QM 1 MM ⫹ 12 U NN ⫹ 2 U NN 兲 ,

共12兲

N el N MM

⫽

QM/MM

h

兺d 兺c

⫺Q c

៝ ⫺r៝ 兩 兩R c d

.

4723

共13兲

QM/MM MM/QM MM , U NN , and U NN The terms JQM/MM, YQM/MM, U NN represent the additional electrostatic interactions between the solute–solvent cluster and polarized surface charges due to the presence of MM atoms. They are interactions between: QM electrons with MM induced surface charges, MM with QM electron induced surface charges, QM nuclei with MM induced surface charges, MM with QM nuclei induced surface charges, and MM with MM induced surface charges. is the direct electrostatic interaction beThe quantity E MM el QM/MM is that between QM nutween the MM atoms, and V NN clei and MM atoms. A subtle issue is that for macromolecuis often calculated with certain cutoff lar systems, E MM el

FIG. 1. Optimized structures of glycine in the gas phase and representative optimized structures of glycine– water clusters. ‘‘N’’ stands for the neutral form, while ‘‘Z’’ stands for the zwitterionic form. The gas phase results are from HF/6-31⫹⫹G共d,p兲 calculations, and C s symmetry was assumed for the zwitterionic form. The cluster results are from HF/6-31 ⫹⫹G共d,p兲/MM calculations, in which the water molecules were treated with a modified TIP3P model. The configurations of the water molecules prior to optimization were collected from a short 共10 ps兲 molecular dynamics calculations for glycine in a 16 Å water droplet. The starting geometries for the optimization of neutral glycine– water clusters are the optimized zwitterionic glycine–water clusters. See the text for more details.

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

4724

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Qiang Cui

schemes to reduce the computational cost,39 and thus it seems appropriate to introduce similar cutoff schemes when MM calculating U NN . Whether this should be introduced for QM/MM MM/QM QM/MM QM/MM , Y , U NN , U NN depends on whether a J QM/MM QM/MM and V NN ; e.g., no cutoff was cutoff is used in h used in most QM/MM methods implemented in CHARMM, except for the original work of Field et al.27共a兲 Since MM atoms are represented by fixed partial charges in the current model, the expressions for all the PCM-related quantities take the same form as those involving nuclear charges, i.e., Eq. 共7兲. This also applies to the derivative calculations, e.g., Eqs. 共10兲 and 共11兲, except that MM atoms do not contribute to terms related to the derivatives of the atomic basis functions 共␮,␯兲 due to the lack of basis functions on them. For example, only the second term in Eq. 共10a兲 has to be evaluated for the MM derivative of JQM/MM. The QM/MM/PCM method was implemented into the simulation package CHARMM,40 which was previously interfaced with GAMESS for ab initio QM/MM calculations.41 The PCM model in GAMESS was modified to include electrostatic QM/MM MM/QM , U NN , and solvation terms (JQM/MM, YQM/MM, U NN MM U NN ) and their analytical first derivatives. Thus the QM/ MM/PCM approach can be used for single point solvation free energy calculations, Monte Carlo simulations and geometry optimization as well as molecular dynamics calculations for reactions in solution. With the recent direct implementation of PCM in GAMESS,4 the QM/MM/PCM approach can be straightforwardly extended to study reactions in large systems such as enzymes. For these large systems, although MC or MD simulations can in principle be carried out, practical applications are likely to be limited to reaction path type of calculations;43 the method can also be used to obtain solvation corrections in a similar fashion as in previous studies,29 but with a more explicit description on effects of long-range electrostatics on the QM region.

III. TEST CALCULATIONS

In this section, we test the current implementation of QM/MM/PCM with two simple but representative molecular systems, i.e., water-dimer and the two solvated forms 共neutral and zwitterionic兲 of glycine in water. Although the waterdimer model is not an example in which the entire first solvation shell is included, it is chosen here to examine the reliability of QM/MM/PCM for describing solvation effects on hydrogen bonding interactions. The QM/MM/PCM results are compared with various calculations with pure QM or pure MM both in the gas phase and in solution with implicit solvent calculations. In the QM/MM calculations, one of the water molecules in the water-dimer and the first solvation shell of glycine 共eight water molecules兲 were treated with MM, and the remaining solute was treated with QM. The QM level chosen here is Hartree–Fock with double-zeta plus polarization quality basis set 共see tables for specific basis sets兲, and explicit MM solvent molecules are treated with a modified44 TIP3P model.45 The van der Waals radii for the QM atoms in QM/MM calculations were chosen from the 46 CHARMM22 force field based on similar atom types. To define the cavity that contains the solute or the solute–solvent cluster in the continuum model calculations, the Pauling set of radii47 was used. In the case of water-dimer, different sets of radii 共Pauling, CHARMM van der Waals, and atomic Born48兲 were tested for the MM water. All QM and QM/MM implicit solvent calculations were carried out with the PCM module in the CHARMM/GAMESS package,41 and pure MM solvation calculations were made with the finite-difference Poisson–Boltzmann module in CHARMM.49 As to the geometries, the water-dimer was fully optimized in all the calculations, including the QM/MM/PCM descriptions. To sample the configurations of the water molecules surrounding the glycine, a short QM/MM molecular

TABLE I. Optimized structure and solvation energy of water-dimer at different level of calculations.a Levelb

R H 共Å兲

D R O-H 共Å兲

A 共Å兲 R O-H

␪ O-H-O 共°兲

␶共°兲

⌬G ele 共kcal/mol兲

QM

2.039 1.978 1.787 关1.748兴 2.044 2.027 共2.027兲 关2.052兴 2.088 2.057 共2.065兲 关1.980兴

0.947 0.952 0.976 关0.982兴 0.971 0.975 共0.974兲 关0.973兴 0.950 0.952 共0.952兲 关0.954兴

0.944 0.946 0.961 关0.976兴 0.943 0.944 共0.944兲 关0.944兴 0.960 0.967 0.968兲 关0.963兴

172.2 176.4 172.2 关177.3兴 174.8 174.9 共175.9兲 关169.9兴 175.7 175.1 共176.1兲 关177.2兴

117.6 122.0 162.8 关175.6兴 133.9 130.5 共129.9兲 关132.1兴 164.4 178.1 共177.0 关174.9兴s

••• ⫺11.4 ••• 关⫺15.6兴 ••• ⫺9.8 共⫺8.7兲 关⫺6.5兴 ••• ⫺9.7 共⫺9.5兲 关⫺7.0兴

MM QM/MM

MM/QM

a

The QM level is HF/6-31G共d,p兲, and MM is a modified version of TIP3P in CHARMM. In the QM/MM 共MM/QM兲 calculations, the hydrogen bond acceptor 共donor兲 is treated with QM. For the definition of geometrical parameters, see Scheme I. b For each level, the first row contains the values in the gas phase, and the rest are obtained from implicit solvent calculations; the latter were carried out with PCM for both pure QM and QM/MM calculations, and with finite difference Poisson–Boltzmann for pure MM. In the implicit solvent calculations, the QM part is always described with the Pauling set of atomic radii; for the MM atoms, three sets of radii were used: Pauling 共without parentheses: r O⫽1.5 Å, r H⫽1.2 Å兲, CHARMM 共with parentheses: r O⫽1.768 Å, r H⫽0.22 Å兲, Born radii from Ref. 48 共with brackets: r O⫽2.2 Å, r H⫽0.0 Å兲. With pure MM, only the results with Born radii are in reasonable agreement with QM results and thus only these are included. See text for other details.

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Combining implicit solvation models with QM/MM methods

dynamics simulation at 300 K was carried out in which a zwitterionic form of glycine was solvated by a 16 Å water sphere that is subject to the deformable boundary condition.50 Eleven snapshots were collected from a 10 ps run, and only the eight closest water molecules 共‘‘first solvation shell’’兲 were kept in the subsequent calculations. These structures were then optimized with the QM/MM method in the gas phase. The optimized structures were then used as the starting geometry for optimizing the neutral glycine–water clusters. Selected optimized structures are shown in Fig. 1. The QM/MM/PCM single point energy calculations were then carried out at these QM/MM optimized glycine–water clusters. Such a protocol was chosen to allow a fair comparison with the recent work of Gordon and co-workers on the same system with a QM/EFP/Onsager or QM/EFP/PCM model.23 For the same reason, only the electrostatic solvent contributions have been considered in the current work; this also makes it easier to compare with Poisson–Boltzmann results for pure MM models 共see below兲.

TABLE II. The energy difference between the zwitterionic and neutral forms of glycine with different calculations.

A. Water-dimer

As shown in Table I, the hydrogen bonding distance becomes shorter by 0.06 Å upon solvation, and the solvation energy is 11.4 kcal/mol at the HF/6-31G共d,p兲 level with the PCM model. With the pure molecular mechanics models, the hydrogen bond distance is much shorter 共1.787 versus 2.039 Å兲 compared to Hartree–Fock results, because the TIP3P model45 was parametrized for condensed phase properties of water. The MM/PB results are very sensitive to the choice of atomic radii. The solvation free energy ranged from ⫺15 kcal/mol to more than ⫺50 kcal/mol with different choices, and only the geometry obtained with the atomic Born radii48 was found to be reasonable. With this set of radii, the hydrogen bond distance is 1.748 Å, and the solvation free energy is 15.6 kcal/mol. In both the gas phase and solution, the optimized water-dimer have more planar structures at the pure MM level 共scheme 1兲, with ␶ close to 170° compared to the value of 120° in the full QM and QM/PCM calculations.

4725

Model

⌬E(E/Z) vacuum

⌬G(E/Z) PCMa

No explicit waterb 8 water S1c 8 water S2c 8 water S3c 8 water S4c 8 water S5c 8 water S6c 8 water S7c 8 water S8c 8 water S9c 8 water S10c 8 water S11c Avg 共SD兲d

⫹31.1 ⫺1.3 ⫹2.1 ⫺4.5 ⫺1.3 ⫺12.4 ⫺9.6 ⫺8.5 ⫺7.4 ⫺8.4 ⫺7.1 ⫺3.6 ⫺5.6 共4.3兲

⫹3.8 ⫺6.2 ⫺6.4 ⫺5.6 ⫺5.8 ⫺7.2 ⫺7.7 ⫺8.1 ⫺6.5 ⫺7.4 ⫺8.6 ⫺7.0 ⫺7.0 共0.9兲

a

In the PCM results, only the electrostatic contributions have been taken into account. b The calculations were performed at the HF/6-31⫹⫹G共d,p兲 level; C s symmetry was imposed for the geometry optimization of the zwitterionic form in vacuum, which otherwise converges to the charge-neutral form. c The glycine was treated by HF/6-31⫹⫹G共d,p兲, and the water molecules were treated with a modified flexible TIP3P model 共Refs. 44 and 45兲. The eleven configurations for the glycine–water clusters were collected from a short QM/MM molecular dynamics simulation. For more calculation details, see the text. For representative structures, see Fig. 1. d The numbers without parentheses are the average values, and those with parentheses are the standard deviations.

hydrogen bond acceptor is QM, and is closer to the pure MM results 共⬃170°兲 when the hydrogen bond donor is treated with QM. The solvation free energy is close to 10 kcal/mol when the Pauling set of radii was used for both QM and MM atoms, which is rather close to the value 共11.4兲 from full QM calculations. With other radii for the MM atoms, the solvation free energy becomes smaller, in a trend similar to that found for pure MM calculations. The results demonstrate that the QM/MM/PCM method gave a rather satisfactory description of solvation effects on the geometry and energetics of hydrogen bonding interactions.

B. Solvated glycine in water

SCHEME 1.

With the QM/MM calculations, the results vary somewhat depending on whether the hydrogen bond donor or acceptor is treated with QM. The solvation results depend also on the choice of atomic radii for the MM atoms 共the Pauling radii were always used for the QM atoms兲, although the variations are much smaller than pure MM results. The hydrogen bond distance is shorter by 0.02–0.03 Å upon solvation, compared to the values of 0.06 and 0.04 Å from full QM and full MM calculations, respectively. The ␶ angle 共⬃130°兲 is close to the full QM results 共⬃120°兲 when the

As mentioned in Sec. I, the relative stability of the neutral and zwitterionic form of glycine in water was found to be very sensitive to the treatment of solvation.23 Although it is known that the zwitterionic form is more stable than the neutral form in solution,51 all implicit solvent models tested in Ref. 23共a兲 were found to favor the neutral form when no explicit solvent molecule was included in part because all three solvation models 共Onsager,26 isodensity PCM,52 SM5.2R53兲 included only linear polarization effects 共the SMx model uses the surface-tension terms to empirically account for first-solvation shell effects兲; a similar result was found here with the PCM model in GAMESS 共Table II兲. When eight water molecules were included to explicitly solvate the glycine with either Hartree–Fock or EFP but neglecting the bulk solvation effects, the neutral form is also more stable by 3–5 kcal/mol.23 Even if the bulk solvation effects were introduced for the glycine–water cluster with a simple solvation

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

4726

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

model such as Onsager26 or SM5.2R,53 the neutral form was still found to be more stable by several kcal per mole. Only when the bulk solvation effect was described with a more sophisticated treatment like the Isodensity-PCM52 the zwitterionic form became more stable than the neutral form by about 8 kcal/mol,23共a兲 which is close to the experimental estimate of 7.7 kcal/mol.51 The problem was also studied by Troung and co-workers24共a兲 with the generalized COSMO model,54 and it was found that it is crucial to optimize the geometry in solution. Without any explicit solvent molecules, however, the zwitterionic form was found to be more stable by only about 3 kcal/mol.24共a兲 A similar value was found by Tomasi and co-workers23共c兲 using B3LYP/6-31G共d兲 with the Integral Equation Formalism 共IEF兲-PCM.11共c兲 QM/MM calculations for the glycine–water cluster without any implicit solvent gave structures qualitatively similar to those in Refs. 23 and 24 共see Fig. 1兲. Although the geometry of the neutral form 关Figs. 1共aN兲, 1共bN兲, 1共cN兲兴 does not change much upon solvation, the geometry for the zwitterionic form in the water cluster 关Figs. 1共bZ兲 and 1共cZ兲兴 is rather different from that of the isolated glycine 关Fig. 1共aZ兲兴. For example, both C–C and C–N distances become shorter by 0.02–0.03 Å upon solvation by the water cluster, and the O–C–O angle becomes smaller by 4°–5°. All these are consistent with the trends observed in the geometry optimizations with the generalized COSMO calculations,24共a兲 which were found to be important for reproducing the correct relative stability of the two forms of glycine. The hydrogen bonds between water molecules vary from 1.7 to 2 Å, which are likely to be too short compared to typical ab initio results because these water molecules are treated with the TIP3P model which was parametrized to reproduce the property of liquid water.45 The hydrogen bond distances between the QM glycine and MM water molecules, however, are on the order of 2 Å, which is quite consistent with those found in previous studies23 in which the water is represented by either Hartree–Fock or effective fragment potentials.25 Unlike the full QM or QM/EFP results for the glycine–water cluster in Ref. 23, the QM/MM calculations predicted that the zwitterionic form is more stable in most structures sampled here, most likely because of the more extensive conformational sampling. On average, the zwitterionic form was found to be more stable by 5.6 kcal/mol, which is already fairly close to the experimental estimate of 7.7 kcal/mol.51 There are, however, substantial variations in the results depending on the structure of the cluster, and the rms fluctuation is 4.3 kcal/ mol for the eleven selected configurations. In one case, the neutral form was found to be more stable by 2.1 kcal/mol, because the ⫺NH⫹ 3 group is not well solvated by the water molecules included 关see Fig. 1共bZ兲兴. The electrostatic interaction energy between the MM water molecules also makes substantial contributions to the relative stability of the two forms, due to the fact that certain water molecules rearrange their orientations when one proton is moved from the ⫺NH⫹ 3 group to the carboxylate 关e.g., W 1 and W 2 in Fig. 1共cZ兲/ 共cN兲兴. When effects of the bulk solvent were included with PCM, the QM/MM/PCM calculations predicted that the zwitterionic form is more stable than the neutral form by 7.0

Qiang Cui

kcal/mol on average, which is in very good agreement with the experimental value of 7.7 kcal/mol.51 It is interesting to note that the fluctuations in the QM/MM/PCM calculations are much reduced compared to the gas phase QM/MM results, which is 0.9 and 4.3 kcal/mol, respectively. This reflects that fact that the glycine was better solvated when the continuum model was used. For example, the zwitterionic form of S2 is less stable than the neutral form by 2.1 kcal/ mol because the ⫺NH⫹ 3 group is not well solvated in the gas phase cluster; the correct trend is reproduced when PCM is used to describe the solvents surrounding the glycine–water cluster 共Table II兲. In S5, the zwitterionic form is 12.4 kcal/ mol more stable than the neutral form; upon solvation with PCM, the energy difference between the two became 7.2 kcal/mol 共Table II兲, which is closer to the experimental value. Due to the small number of configurations sampled here and the fact that only the electrostatic component of the solvation energy was included, the close agreement with the QM/MM/PCM result has to be considered as fortuitous. For example, it was argued that electron correlation might be important for the relative stability as well, although the magnitude given in Ref. 23共b兲 could be an overestimation considering the limited configuration sampling. A more solid analysis is left for the future, which requires a QM/MM/ PCM implementation with correlated QM approaches and more extensive configurational sampling. Finally, we briefly comment on the computational cost of QM/MM/PCM calculations. When only electrostatic solvation free energy is calculated, introducing MM atoms does not cause a large increase in computational cost. For example, one solvation free energy calculation for isolated glycine at the HF/6-31G共d兲 level 共a basis set smaller than the one used in the above-discussed calculations such that pure QM/PCM calculations can be done quickly兲 took about 0.5 min on an Athlon machine at 1.2 GHz, and one HF/MM/ PCM calculation with eight MM water molecules took about 1.5 min. By contrast, the calculation took more than 13 min when all atoms were treated with Hartree–Fock. The efficiency of the current implementation is essentially limited by the speed and scaling of the continuum model, especially for gradient calculations. With the recent direct implementation of the PCM model,42 the QM/MM/PCM approach can be straightforwardly extended to macromolecules, which is our long-term goal that inspired the current work. IV. CONCLUSIONS

A combined approach in which the solute and a number of solvent molecules are described with the hybrid QM/MM method and the bulk solvent is represented with an implicit solvation model was implemented into the simulation package CHARMM. A typical application of such an approach would treat the solute molecule with a QM approach, the first 共several兲 solvation shell共s兲 with MM models, and the bulk solvent with a continuum model. When the system is a macromolecule such as an enzyme, hybrid QM/MM methods can be used to describe the solute plus a number of explicit solvent molecules, leaving the bulk as an implicit solvent

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Combining implicit solvation models with QM/MM methods

model. The advantage of such a multilayer treatment is that both effects due to the first solvation shell 共e.g., Pauli repulsion, dispersion, solvent entropic contributions related to hydrophobic interactions as well as nonlinear response of the solvent to the solute charge distribution兲 and long-range electrostatics due to the bulk solvent can be taken into account. By judiciously choosing a solute–solvent supermolecule and an appropriate QM/MM partition, such a QM/MM/PCM approach can offer both computational accuracy and efficiency. With only a small number of explicit solvent molecules, calculations for quantities that require configurational samplings can achieve faster convergence compared to full microscopic simulations. The method was tested with two simple but representative systems: water-dimer and two forms of glycine in water. It was shown that the QM/MM/PCM model gave satisfactory descriptions on the effect of solvation on the structure and energetics of these systems. It is particularly encouraging that the QM/MM/PCM calculation correctly reproduced the relative stability of the neutral and zwitterionic forms of glycine in solution, while many previous attempts gave the opposite trend unless sophisticated solvation models were used.23,24 Although the current implementation is most suited for applications toward medium size systems in solution due in part to the steep scaling of the matrix-inversion formulation of PCM11共b兲 used here, it can be extended to macromolecular systems with the direct implementation of PCM.42 A similar philosophy can also be used to combine QM/MM methods with other approaches that deal with solvations, such as the generalized Born model7,31 or density functional theory,19,20 depending on the problem at hand. For macromolecular applications, it is also of great interest to combine QM/MM with simplified representations for some of the protein or nucleic acid atoms in addition to the bulk solvent, such as in the methods pioneered by Warshel et al.15 and the related work of Roux and co-workers.55 These developments are currently in progress.

共a兲 J. D. Jackson, Classical Electrodynamics 共Wiley, New York, 1962兲; 共b兲 B. Honig and A. Nicholls, Science 268, 1144 共1995兲. 7 共a兲 W. C. Still, A. Tempczyk, R. C. Hawley, and T. Hendrickson, J. Am. Chem. Soc. 112, 6127 共1990兲; 共b兲 M. Schaefer and M. Karplus, J. Phys. Chem. 100, 1578 共1996兲. 8 B. Jayaram, K. A. Sharp, and B. Honig, Biopolymers 28, 975 共1989兲. 9 P. Debye and E. Hu¨ckel, Phys. Z. 24, 305 共1923兲. 10 V. K. Misra, K. A. Sharp, R. A. Friedman, and B. Honig, J. Mol. Biol. 238, 245 共1994兲. 11 共a兲 S. Miertus, E. Scrocco, and J. Tomasi, Chem. Phys. 55, 117 共1981兲; 共b兲 R. Cammi and J. Tomasi, J. Comput. Chem. 16, 1449 共1995兲; 共c兲 B. Mennucci, E. Cance`s, and J. Tomasi, J. Phys. Chem. B 101, 10506 共1997兲; 共d兲 E. Cance`s, B. Mennucci, and J. Tomasi, J. Chem. Phys. 107, 3032 共1997兲. 12 A. Klamt and G. Schu¨u¨rmann, J. Chem. Soc., Perkin Trans. 2 1993, 799 共1993兲. 13 F. M. Floris, A. Tani, and J. Tomasi, Chem. Phys. 169, 11 共1993兲. 14 共a兲 B. Linda, Adv. Chem. Phys. 12, 225 共1967兲; 共b兲 D. Rinaldi, B. J. Costa Cabral, and J. L. Rivail, Chem. Phys. Lett. 125, 495 共1986兲; 共c兲 F. Floris and J. Tomasi, J. Comput. Chem. 10, 616 共1989兲. 15 共a兲 A. Warshel and S. T. Russell, Q. Rev. Biophys. 17, 283 共1984兲; 共b兲 J. Floria`n and A. Warshel, J. Phys. Chem. B 101, 5583 共1997兲. 16 J.-P. Hansen and I. R. McDonald, Theory of Simple Liquids, 2nd ed. 共Academic, New York, 1986兲. 17 共a兲 H. C. Andersen and D. Chandler, J. Chem. Phys. 57, 1918 共1972兲; 共b兲 D. Chandler and H. C. Andersen, ibid. 57, 1930 共1972兲; 共c兲 L. R. Pratt and D. Chandler, ibid. 67, 3683 共1977兲; 共d兲 F. Hirata and P. Rossky, ibid. 78, 4133 共1983兲. 18 See, for example, 共a兲 B. C. Perng, M. D. Newton, F. O. Raineri, and H. L. Friedman, J. Chem. Phys. 104, 7153 共1996兲; 共b兲 H. Sato, F. Hirata, and S. Kato, ibid. 105, 1546 共1996兲; 共c兲 D. Beglov and B Roux, ibid. 104, 8678 共1996兲. 19 See, for example, 共a兲 C. N. Patra and A. Yethiraj, Biophys. J. 78, 699 共2000兲; 共b兲 J. Phys. Chem. B 103, 6080 共1999兲. 20 共a兲 J. K. Percus, in The Equilibrium Theory of Classical Fluids, edited by H. L. Frisch and J. L. Lebowitz 共Benjamin, New York, 1964兲; 共b兲 D. Chandler, J. D. McCoy, and S. J. Singer, J. Chem. Phys. 85, 5971 共1986兲; 共c兲 R. Evans, in Fundamentals of Inhomogeneous Fluids, edited by D. Henderson 共Dekker, New York, 1992兲; 共d兲 A. Yethiraj and C. E. Woodward. J. Chem. Phys. 102, 5499 共1995兲. 21 Q. Cui and M. Karplus, J. Phys. Chem. B 104, 3721 共2000兲. 22 B. Mennucci, J. Am. Chem. Soc. 124, 1506 共2002兲. 23 共a兲 P. Bandyopadhyay and M. S. Gordon, J. Chem. Phys. 113, 1104 共2000兲; 共b兲 P. Bandyopadhyay, M. S. Gordon, B. Mennucci, and J. Tomasi, ibid. 116, 5023 共2002兲; 共c兲 L. Gontrani, B. Mennucci, and J. Tomasi, J. Mol. Struct.: THEOCHEM 500, 113 共2000兲. 24 共a兲 T. N. Troung and E. V. Stefanovich, J. Chem. Phys. 103, 3709 共1995兲; 共b兲 F. R. Tortonda, J. L. Pascual-Ahuir, E. Silla, and I. Tunon, Chem. Phys. Lett. 221, 260 共1996兲. 25 P. N. Day, J. H. Jensen, M. S. Gordon, S. P. Webb, W. J. Stevens, M. Krauss, D. Garmer, H. Basch, and D. Cohen, J. Chem. Phys. 105, 1968 共1996兲. 26 J. G. Kirkwood, J. Chem. Phys. 2, 351 共1934兲. 27 共a兲 M. J. Field, P. A. Bash, and M. Karplus, J. Comput. Chem. 11, 700 共1990兲; 共b兲 V. Gogonea, D. Suarez, A. van der Vaart, and K. W. Merz, Curr. Opin. Struct. Biol. 11, 217 共2001兲; 共c兲 J. Gao, in Reviews in Computational Chemistry, edited by K. B. Lipkowitza and D. B. Boyd 共VCH, ˚ qvist and A. Warshel, Chem. Rev. Weinheim, 1996兲, Vol. 7, p. 119; 共d兲 J. A 93, 2523 共1993兲. 28 C. L. Brooks III and M. J. Karplus, J. Mol. Biol. 208, 159 共1989兲. 29 共a兲 M. Formaneck, G. Li, X. Zhang, and Q. Cui, J. Compt. Theor. Chem. 共in press兲; 共b兲 A. R. Dinner, X. Lopez, and M. Karplus, Theor. Chem. Acc. 共in press兲. 30 T. Simonson, G. Archontis, and M. Karplus, J. Phys. Chem. B 101, 8349 共1997兲. 31 E. Pellegrini and M. J. Field, J. Phys. Chem. B 106, 1316 共2001兲. 32 M. W. Schmidt, K. K. Baldridge, J. A. Boatz et al., J. Comput. Chem. 14, 1347 共1993兲. 33 T. Vreven, B. Mennucci, C. O. de Silva, K. Morokuma, and J. Tomasi, J. Chem. Phys. 115, 62 共2001兲. 34 R. A. Pierotti, Chem. Rev. 76, 717 共1976兲. 35 R. Cammi, B. Mennucci, and J. Tomasi, J. Chem. Phys. 110, 7627 共1999兲. 36 共a兲 R. Cammi and J. Tomasi, J. Chem. Phys. 101, 3888 共1994兲; 共b兲 E.

ACKNOWLEDGMENTS

The current research is supported by a start-up fund from the Department of Chemistry and College of Letters and Science at University of Wisconsin, Madison. Q.C. wishes to thank Professor J. H. Jensen for interesting discussions related to recent developments of PCM in GAMESS and A. Van Wynsberghe and M. Formaneck for critically reading the manuscript. He also thanks Professor Martin Karplus for support in the initial stage of this work. For reviews, see 共a兲 J. Tomasi and M. Persico, Chem. Rev. 94, 2027 共1994兲; 共b兲 C. J. Cramer and D. G. Truhlar, ibid. 99, 2161 共1999兲. 2 See, for example, A. Warshel, Computer Modeling of Reactions in Enzymes and Solution 共Wiley, New York, 1997兲. 3 See, for example, C. L. Brooks III, M. Karplus, and M. B. Pettitt, Adv. Chem. Phys. 71, 1 共1987兲. 4 See, for example, 共a兲 T. Lazaridis and M. Karplus, Proteins: Struct., Funct., Genet. 35, 133 共1999兲; 共b兲 B. Roux and T. Simonson, Biophys. Chem. 78, 1 共1999兲. 5 共a兲 C. Tanford, Proc. Natl. Acad. Sci. U.S.A. 76, 4175 共1979兲; 共b兲 K. Sharp, A. Nicholls, R. Friedman, and B. Honig, Biochemistry 30, 9696 共1991兲. 1

4727

6

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp

4728

J. Chem. Phys., Vol. 117, No. 10, 8 September 2002

Cance`s and B. Mennucci, ibid. 109, 249 共1998兲; 共c兲 E. Cance`s, B. Mennucci, and J. Tomasi, J. Chem. Phys. 109, 260 共1998兲. 37 Y. Yamaguchi, H. F. Schaefer III, Y. Osamura, and J. D. Goddard, New Dimension in Quantum Chemistry: Analytical Derivative Methods in Ab Initio Molecular Electronic Structure Theory 共Oxford University Press, New York, 1994兲. 38 See, for example, D. Mennucci and J. Tomasi, J. Chem. Phys. 106, 5151 共1997兲. 39 P. J. Steinbach and B. R. Brooks, J. Comput. Chem. 15, 667 共1994兲. 40 B. R. Brooks, R. E. Burccoleri, B. D. Olafson, D. J. States, S. Swaminathan, and M. Karplus, J. Comput. Chem. 4, 187 共1983兲. 41 P. D. Lyne, M. Hodoscek, and M. Karplus, J. Phys. Chem. A 103, 3462 共1999兲. 42 J. H. Jensen 共private communication兲. 43 See, for example, 共a兲 Q. Cui and M. Karplus, J. Am. Chem. Soc. 122, 2284 共2001兲; 共b兲 J. Phys. Chem. B 106, 1768 共2002兲.

Qiang Cui E. Neria, S. Fischer, and M. Karplus, J. Chem. Phys. 105, 1902 共1996兲. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey, and M. L. Klein, J. Chem. Phys. 79, 926 共1983兲. 46 A. D. MacKerell, Jr., D. Bashford, M. Bellott et al., J. Phys. Chem. B 102, 3586 共1998兲. 47 Handbook of Chemistry and Physics, 78th ed. 共CRC, New York, 1997兲. 48 M. Nina, D. Beglov, and B. Roux, J. Phys. Chem. 101, 5239 共1997兲. 49 W. Im, D. Beglov, and B. Roux, Comput. Phys. Commun. 109, 1 共1998兲. 50 C. L. Brooks III and M. Karplus, J. Chem. Phys. 79, 6312 共1983兲. 51 共a兲 P. Haberfield, J. Chem. Educ. 57, 346 共1980兲; 共b兲 R. C. Gaffney and J. S. Pierce, J. Am. Chem. Soc. 99, 4293 共1977兲. 52 J. B. Foresman, T. A. Keith, K. B. Wiberg, J. Snoonian, and M. J. Frisch, J. Phys. Chem. 100, 16098 共1996兲. 53 J. Li, J. Cramer, and D. G. Truhlar, Chem. Phys. Lett. 288, 293 共1998兲. 54 T. N. Truong and E. V. Stefanovich, Chem. Phys. Lett. 240, 253 共1995兲. 55 W. Im, S. Berne`che, and B. Roux, J. Chem. Phys. 114, 2924 共2001兲. 44 45

Downloaded 17 Sep 2002 to 128.104.71.49. Redistribution subject to AIP license or copyright, see http://ojps.aip.org/jcpo/jcpcr.jsp