Ab Initio Frozen Density Functional Calculations of Proton Transfer

In general one may define a free energy function Δga(X) that expresses the ... and the corresponding changes of X play a major role in chemical proce...
2 downloads 0 Views 401KB Size
Download PDF
15444

J. Phys. Chem. 1996, 100, 15444-15449

Ab Initio Frozen Density Functional Calculations of Proton Transfer Reactions in Solution Tomasz Wesolowski,† Richard P. Muller, and Arieh Warshel* Department of Chemistry, UniVersity of Southern California, Los Angeles, California 90089-1062 ReceiVed: April 9, 1996; In Final Form: July 10, 1996X

The recently developed frozen density functional theory (FDFT) is extended to ab initio free energy calculations of chemical reactions in solution. This method treats the solute-solvent system as a supermolecule but constrains the electron density of the solvent molecules. Unlike hybrid quantum mechanical/molecular mechanics techniques, FDFT represents the solvent quantum mechanically. The quality of the solute-solvent interaction potential is examined by generating clusters of a reacting system and several solvent molecules and comparing the supermolecule DFT energies to the corresponding FDFT energies. The FDFT potential surfaces for solute-solvent systems provide a good approximation of the supermolecule DFT surfaces and require, in some cases, several orders of magnitude less computation time (in particular if one treats many solvent molecules quantum mechanically). The ab initio free energy surface for the F- + HF f FH + Fproton transfer reaction in solution is calculated using the corresponding “classical” empirical valence bond (EVB) potential surface as a reference potential. The encouraging results indicate that FDFT can be used to study chemical reactions in solution, capturing the quantum mechanical aspects of the solvent, which is not possible using hybrid quantum mechanical/molecular mechanics approaches. Furthermore, the use of EVB as a reference potential is found to be an extremely effective way of obtaining ab initio free energies for chemical processes in solution or in clusters.

I. Introduction Quantum mechanical calculations of chemical processes in solution and proteins present a major challenge owing to the enormous size of the reacting system.1 Significant advances have been made in these calculations using hybrid quantum mechanical/molecular mechanics (QM/MM) models that treat a small part of the system (e.g., the solute) quantum mechanically and the rest of the system (e.g., the solvent) classically.1-20 Although these approaches are very useful, they leave open major questions about the approximations involved in treating the solvent classically. Representing the solvent on the same quantum mechanical level as the solute (a supermolecule treatment) could solve this problem, but such an approach is not practical even with the current generation of computers. Our attempts to address this issue have involved the development of a pseudopotential representation of the solvent21 and, more recently, the development of the frozen density functional approach (FDFT).22,23 This general idea can be further refined by removing the requirement of “frozen densities” as in approaches where the solvent density is allowed to change using criteria such as induced dipoles or by iteratively replacing the unfrozen sites.24 Such approaches are referred to here as constrained DFT (CDFT) methods. Other researchers have also developed DFT approaches where part of the electron density is frozen or constrained.25,26 The FDFT approach has been used in the preliminary studies that establish its validity in ab initio free energy perturbation (FEP) calculations of the solvation free energy. In the present work we extend the FDFT approach to studies of chemical reactions in solution, considering a proton transfer reaction as a test case. Quantum mechanical studies of proton transfer in solution have been reported by several workers who used hybrid molecular orbital/molecular mechanics (MO/MM) approaches,1-5,9 † Permanent address: University of Geneva, Physical Chemistry Department, 30, quai Ernest-Ansermet, CH-1211 Geneva 4, Switzerland. X Abstract published in AdVance ACS Abstracts, September 1, 1996.

S0022-3654(96)01068-4 CCC: $12.00

the hybrid valence bond/classical empirical valence bond (EVB) method,1,7-9,27 hybrid molecular orbital/reaction field methods,16 and a recent DFT/classical study.28 However, none of these studies treats the solvent quantum mechanically, and this challenging problem is the subject of the present study. II. Methods Studies of chemical processes in solution require a convenient method for describing the potential surface of the solute-solvent system and an efficient method for evaluating the free energy of this system along the relevant reaction coordinate. Both of these requirements present a major challenge when the calculations are done on an ab initio level and when the quantum mechanical features of the solvent are taken into account. The present work incorporates the quantum mechanical features of the solvent in the potential energy surface of the system using the FDFT approach. The specific formulation will be outlined in section II.A below. Having a reliable ab initio surface for reactions in solution does not necessarily mean that the corresponding free energy barriers can be evaluated accurately. The solute-solvent system has enormous dimensionality, and the relevant reaction coordinate must involve both the solute and solvent coordinates. Free energy perturbation approaches that use the solute coordinate as the mapping parameter do not allow one to sufficiently sample transition state configurations to find the proper transition state (see discussion in ref 5). What is needed for effective calculation of the free energy is a mapping potential that forces both the solute and the solvent coordinates to their corresponding transition state values. Section II.B will describe the use of the EVB potential surface for this purpose. A. Potential Energy Surface of the Solute-solvent System. 1. FDFT Potential Energy Surface for Reactions in Solutions. The potential energy surface of the solute-solvent system is represented by the recently developed frozen density functional method22,23 (see also ref 26 for a closely related method). This method treats the solute-solvent system as a © 1996 American Chemical Society

Proton Transfer Reactions

J. Phys. Chem., Vol. 100, No. 38, 1996 15445

Figure 2. Potential energy curve for the proton transfer process in an isolated (FHF)- where the F‚‚‚F distance is 3.4 Å. Because the potential surface is symmetric, we report only the results for proton transfer from the ground to the transition state. Figure 1. Schematic description of the solvated (FHF)- system studied in this work. The electronic density of the solvent molecules is kept frozen, while the electron density of the (FHF)- solute system is allowed to change to minimize the total energy of the system.

supermolecule in the DFT formulation29 but freezes the electron density of the solvent molecules and computes the effect of these solvent molecules on the solute Hamiltonian at each nuclear configuration of the system. The effective potential used for the evaluation of the solute electron density reflects both the solute and the solvent electron densities and includes a special nonadditive kinetic energy functional that couples these densities. The FDFT Hamiltonian includes both the long range solute-solvent interaction and the solute-solvent electronic overlap contribution. For more details of the model and its implementation see refs 22 and 23. The validation of the FDFT surface in our specific case will be given in the results section. The present work considers proton transfer in the solvated (FHF)- complex (Figure 1) in the simple case where the proton moves along the straight line connecting both fluorine atoms. Accurate description of the (FHF)- solute system requires a relatively large basis set (i.e., including d and f functions for the F and d function for the H30), which is significantly larger that the used in our study (we used a Whitten31 basis set with s and p functions only). However, the purpose of this work is to explore the performance of the FDFT method in solvation problems rather than to address well-known issues of the quality of solute potential energy surfaces. Thus, we replace our FDFT potential surface by the potential gas E′q(R) ) Eq(R) + E h gas q (R) - Eq (R)

(1)

is the where q designates “quantum mechanical” and E h gas q energy of the isolated (FHF) (gas phase) solute obtained with a large basis set. In the present case we obtain E h qgas by the deMon program of Salahub32 using double-ζ plus polarization and triple-ζ plus polarization on the hydrogen and fluorine atoms, respectively. These calculations reproduce a gas phase proton transfer barrier of 30.9 kcal/mol for the F‚‚‚F distance of 3.4 Å (see Figure 2), in a good agreement with the potential energy surface obtained with high level ab initio calculations of Yamashita et al.30 2. EVB Potential Energy Surfaces for Simulating Reactions in Solution and for Efficient Mapping of ab Initio Surfaces. This work uses the EVB potential surface to map the proper solutesolvent reaction coordinate in our FEP calculations. This type

of potential energy surface is extensively discussed elsewhere,1,5,8 and we will only mention here points that are specific to the F- + HF f FH + F- reaction. This reaction may be described using three resonance configurations:1

F-

ψ1′ ) FsH ψ2′ ) F-

HsF

(2)

ψ3′ ) F- H+ FHere, however, we consider only two effective structures: ψ1, obtained by mixing ψ1′ with ψ3′; ψ2, obtained by mixing ψ2′ with ψ3′. The corresponding effective VB Hamiltonian is given by

[

 H Hevb ) H1  12 12 2

]

(3)

where the relevant matrix elements are given by (1) (1) + USs + Uss 1 ) M(R1) + Vnb(R2) + VQQ (2) (2) 2 ) Vnb(R1) + M(R2) + VQQ + USs + Uss

(4)

H12 ) A e-µR3 and R1, R2, and R3 designate the F(1)‚‚‚H, H‚‚‚F(2), and F(1)‚‚‚F(2) distances, respectively, M is a Morse potential, Vnb is a (i) nonbonded interaction function, VQQ is the Coulombic interac(i) tion between the fragments of resonance structure i, USs is the solute-solvent classical force field in the ith resonance configuration, and Uss is the solvent-solvent force field. The parameters taken for the solute potential functions are given in Table 1.1 of ref 1. The solute charges in each resonance structure are obtained from the gas phase charges of the isolated fragments: -1.0e, 0.447e, -0.447e and -0.447e, 0.447e, -1.0e for F(1)‚‚‚HF(2) and F(1)H‚‚‚F(2), respectively. The solute-solvent potential surface was represented using the standard ENZYMIX force field33 with the given solute charges and with r* and * of 2.501 Å and 110.0 kcal/mol for the nonbonded interactions involving fluorine atoms, and the solvent-solvent force field was represented by the SCAAS model as implemented in the ENZYMIX program.33

15446 J. Phys. Chem., Vol. 100, No. 38, 1996

Wesolowski et al.

B. Calculations of ab Initio Free Energy Surfaces. 1. General Considerations. The determination of rate constants for chemical reactions in solution requires the calculation of the corresponding activation free energies. This can be done by evaluating the free energy along the proper reaction coordinate, which is referred to as the ground state free energy function. The evaluation of the ground state free energy for chemical reactions in solution is far from trivial, since such processes involve enormous numbers of solvent coordinates and strong coupling between the solute-solvent coordinates. In general one may define a free energy function ∆ga(X) that expresses the probability of finding the system at a specified value, Xn, of a given reaction coordinate, X, for a potential surface Ea, as

exp{-β∆ga(Xn)} )

za(Xn) Za

(5)

where β ) 1/kBT (kB is the Boltzmann constant) and the partition functions za and Za are defined by

za(Xn) ) ∫exp{-βEa(Xn)} dγn Za ) ∫za(X) dX

(6)

and γn represents the coordinates perpendicular to X at Xn. It is useful to write za(Xn) as

za(Xn) ) ∫δ(X - Xn) exp{-βEa(X)} dΓ

(7)

where dΓ ) dγ dX. The probability za(Xn)/Za can be evaluated, in principle, by running very long trajectories and evaluating the ratio between the time the system is at Xn to the total time. However, such an approach might be very ineffective when Xn corresponds to high-energy regions. Thus, it is frequently advantageous to rewrite eq 5 as5,34,35

exp{-β∆ga(Xn)} )

)

( )( )

(

Zm za(Xn) Za Zm

)

∫exp{-β(Em - Ea) - βEa} dΓ ∫ exp{-βEa} dΓ

×

( )

za(Xn) (8) Zm

) 〈exp{-β(Em - Ea)}〉Ea

( )

where we have multiplied and divided by the partition function Zm corresponding the potential surface Em, and where 〈〉Ea denotes a MD average over the potential surface Ea, and we used the standard approach of replacing a configuration integral by an average over a trajectory. We may further manipulate this expression using5,34,35

(

)

∫δ(X - Xn) exp{-β(Ea - Em + Em)} dΓ ∫exp{-βEm} dΓ

exp{-β∆ga(Xn)} ) 〈exp{-β(Em - Ea)}〉Ea × 〈δ(X - Xn) exp{-β(Ea - Em)}〉Em ) exp{-β∆Gafm} × 〈δ(X - Xn) exp{-β(Ea - Em)}〉Em (10) This expression is particularly effective when Em is chosen as the potential that forces X to spend most of its time at Xn. The 〈〉Ea average is evaluated by standard free energy perturbation approaches using a series of mapping potentials (see below). The proper coordinate X should reflect the solute and solvent reaction coordinate and can be obtained in a general way by considering the energies 1 and 2 of the two valence bond states given by (ψ1) and (ψ2). Note that the energy difference 1 2 is negative when the system is in a state corresponding to (ψ1), goes through zero at the transition state, and is positive when the system is in the state corresponding to (ψ2). We may therefore monitor the progress of the reaction FH + F- f F+ HF via the energy difference 1 - 2. What is significant about the definition X ) 1 - 2 is that it provides a unique way of representing the many dimensional solvent coordinates because the electrostatic part of X is essentially the product of the solvent reaction field and the dipole moment associated with the change of the solute charges upon transfer from the reactant to the product state.8 Since the solvent fluctuations and the corresponding changes of X play a major role in chemical processes in solution, it is crucial to be able to sample these fluctuations. Unfortunately, the solvent component of the reaction coordinate involves a concerted motion of many coordinates and one cannot simply force all water molecules to rotate in a given direction in the same way that a single bond in a gas phase calculation is forced to stretch. Thus, it is important to use an effective mapping procedure that will use some parametric representation of the change in the solute charges during the reaction to map the solvent coordinates. Such a mapping approach is conveniently provided by the EVB approach as described below. 2. Using the EVB as a Reference Potential. To obtain reasonable convergence in free energy calculations of activation barriers, it is essential to force the system to spend significant time at the transition state. This should be done by using a mapping potential that can polarize the solvent to the desired configuration. Perhaps the most effective way of achieving this is via a valence bond mapping potential4,8 m ) (1 - λm)1 + λm2 ref

za(Xn) Zm

exp{-β∆ga(Xn)} ) 〈exp{-β(Em - Ea)}〉Ea ×

and finally obtain

(9)

(11)

where λm is a variable that changes from 0 to 1 to control the progress of the reaction from 1 to 2, the previously mentioned potential surfaces of the two valence bond resonance structures given in eq 4. Using this potential as Em in eq 10, we obtain8

exp{-β∆g(Xn)} ) exp{-β∆Gref(λ0 f λm)} × m )}〉mref (12) 〈δ(X - Xn) exp{-β(Eq - ref m where ref is the mapping potential from eq 11, λm is the value of λ that keeps the system closest to Xn, and Eq is the given quantum mechanical potential surface, either the ab initio surface or the ground state EVB surface obtained by solving eq 3. The free energy ∆Gref(λ0 f λm) is obtained using standard free

Proton Transfer Reactions

J. Phys. Chem., Vol. 100, No. 38, 1996 15447

energy perturbation (FEP) approaches with a gradual change of λm in eq 11 in m + 1 discrete steps, using8,36

exp{-β∆Gref(λ0 f λm)} ) m

n+1 n 〈exp{-β(ref - ref )}〉 ∑ n)0

n ref

(13)

The reaction coordinate X is usually taken as the energy difference between 2 and 1. The EVB mapping technique described in eqs 11-13 essentially maps around a thermodynamic cycle that allows one to express the quantum free energy barriers as † 0 - ∆greffq ∆g†q ) ∆Gref(λ0 f λ†) + ∆greffq

(14)

Here, 0 0 0 } ) 〈δ(X - X0) exp{-β(Eq - ref )}〉ref (15) exp{-β∆greffq

and

Figure 3. Comparison of DFT and FDFT potential energy curves for the proton transfer in a cluster of (FHF)- and five water molecules. The DFT calculations treat the entire system as a supermolecule, whereas the FDFT calculations freeze the densities of the water molecules.

† † † (16) } ) 〈δ(X - X†) exp{-β(Eq - ref )}〉ref exp{-β∆greffq

One advantage to this approach is that the computationally expensive mapping from the reference potential to the ab initio potential need only be done at the reactant and transition states. A similar idea has already been implemented in solvation studies.6,23 III. Results and Discussion A. Reliability of the FDFT Potential Energy Surface. Since the FDFT approach is aimed at providing the ab initio free energy of very large clusters, it cannot be validated by comparing its results to those obtained by the full DFT treatment, since these calculations are not feasible at this time. It is obviously not useful to compare the FDFT results to those of hybrid QM/MM approaches, since we are addressing the much more challenging problem of representing the solvent quantum mechanically rather than classically. We focus here on establishing the assumption that the supermolecule DFT potential surface can be represented by the FDFT approximation. As a first step, we refer the reader to our earlier studies of these issues where we showed that the FDFT model can reproduce accurately the DFT solute-solvent potential surface. In the second step of our analysis we compare the FDFT and DFT proton transfer potential surfaces in supermolecular systems composed of the (FHF)- ion and clusters of five and six water molecules. Obviously, we could not afford to perform such a comparison on all possible geometries of our arbitrary test cases, and consequently, the potential surface obtained has nothing to do with the free energy surface of proton transfer in the complete system with many solvent molecules that can only be evaluated via the FDFT approach. The arbitrary coordinates of the water molecules in the supermolecular cluster test case are generated by running MD trajectories on the reactant state F(1)H‚‚‚F(2)- and selecting at random times the solvent molecules that are closest to the solute. As a consequence of using the reactant trajectories, we generate conformations where the solvent molecules interact favorably with the F(1)H‚‚‚F(2)- charge distribution. The DFT energies are evaluated using the deMon program of Salahub et al.32 with a local exchange-correlation functional. The FDFT energies are calculated using the same exchange-correlation functional and the FDFT program developed in our lab previously and discussed elsewhere,22 representing the water molecules with their frozen gas phase electron density.

Figure 4. Comparison of supermolecule DFT and FDFT potential energy curves for a proton transfer in a cluster of (FHF)- and six water molecules.

The results of these calculations are summarized in Figures 3 and 4, which depict the contribution of the solute-solvent interaction to the potential surface for proton transfer in each of the clusters studied. As it is clear from the figures, the asymmetric potential is due to the asymmetric field from the solvent molecules. Both FDFT and deMon/DFT methods place the minimum at the point where the proton is attached to the F(1) atom owing to the reactant geometry used to generate the solvent configurations. The difference between both energies is most pronounced when the proton approaches the F(2) atom. This difference reflects the fact that the frozen water molecules cannot change their polarization when the net charge of the F(2) atom changes. One can argue, however, that configurations with the largest deMon/DFT and FDFT differences are not accessible at room temperature MD simulations. The agreement between the FDFT and DFT energies is much better when there is more uniform charge distribution around the solute, as demonstrated in Figure 4. Overall, FDFT provides a reasonable approximation for the supermolecule DFT results, considering the fact that it does not involve any adjustment of parameters (the van der Waals repulsion is generated by the ab initio surface and not by an empirical force field). The agreement clearly can be improved by adding induced dipole terms to the FDFT densities. As a further check of the reliability of the FDFT approximation, we consider MD trajectories of the entire solvent system and compare the ab initio FDFT solute-solvent interaction

15448 J. Phys. Chem., Vol. 100, No. 38, 1996

Figure 5. Fluctuations of the EVB and FDFT energies of a solvated (FHF)- system with a F‚‚‚F distance of 3.4 Å along a highly constrained “trajectory” that involves a change of the proton position from a F‚‚‚H distance of 1.0-1.7 Å and a corresponding change in the EVB charges. The trajectory reflects an arbitrary procedure that involves a change of geometry and charges every 100 time steps (where each step represents 1 fs). Note that the trajectories used for the mapping procedure are much longer than those in the arbitrary procedure of this figure. Regardless of the short time and the arbitrary procedure, the figure indicates that the two potential surfaces are quite similar.

potential to the corresponding EVB classical solute-solvent interaction term. The comparison was done by a rather arbitrary procedure running a constrained trajectory where the position of the proton and the corresponding EVB charges are changed (in a rather short time interval) from their values at the reactant state to their values at the transition state. This procedure is clearly improper for any free energy calculation (since no equilibration is obtained), but it allows one to compare the EVB and FDFT potential surfaces at different solvent configurations relevant for the proton transfer reaction. The results, summarized in Figure 5, demonstrate that the FDFT energies follow the same trends as the classical EVB energies. This is very encouraging, considering the fact that the FDFT energy is an entirely ab initio quantity without any adjustable van der Waals parameters. B. FDFT Free Energy Surface. With reasonable confidence in the reliability of the FDFT potential energy surface we now evaluate the free energy surface for a proton transfer in a solvated (FHF)- system. These calculations fix the F‚‚‚F distance at 3.4 Å, since our main interest in this work is in the computational aspect of the model and not in an analysis of a specific reaction. The results of our free energy calculation are summarized in Figure 6. As seen from the figure, the solvent effect leads to an increase in the activation barrier for the proton transfer process, following the trend from other proton transfer reactions.1,5 However, in the (FHF)- case, it is possible in principle that the stabilization of the F-H+F- VB configuration by the solvent will result in some stabilization of the transition state. In fact in our calculations we find that the solvent contribution to the reaction barrier reaches its maximum before the formation of the symmetric transition state. This effect might still be due to the neglect of the solvent polarization and is outside the scope of the present work, which focuses simply on the ability to evaluate activation free energies with a given ab initio potential energy surface and not on the reliability of the relevant free energy. This issue will be studied using systems where the relevant activation barriers are better characterized experimentally. Regardless of the reliability of our free energy calculations (which still involve significant convergence problems), we have demonstrated here that one can obtain converging free energy

Wesolowski et al.

Figure 6. FDFT free energy surface for a proton transfer in a solvated (FHF)- system. The calculations correspond to a F‚‚‚F distance of 3.4 Å. Because the potential surface is symmetric, we report only the results for proton transfer from the ground to the transition state. The solid line represents a smooth curve fit through these points.

Figure 7. Projection of the dipole moment of the solvated (FHF)- on the F-F axis along a trajectory of the solvent molecules. The trajectory here uses the same arbitrary procedure described in Figure 5; different segments of the trajectory correspond to different positions of the proton in the (FHF)- complex.

surfaces for ab initio calculations that treat the solvent quantum mechanically. It is also important to note that the relevant activation barrier can be evaluated by using the thermodynamic cycle of eq 14. This, and the use of the EVB reference potential, make our approach much more effective than the direct use of the ab initio FDFT potential energy surface in the MD simulations. Another useful check of our calculations involves assessing the effect of the solvent on the solute electron distribution as expressed in the change in the solute dipole moment with solvation, shown in Figure 7, which compares the solute dipole moment in solution and in vacuum for different fixed values of F(1)‚‚‚H distance. The solute dipole moment changes with changing environment for each conformation of the solute. The solute dipole moment systematically increases relative to its gas phase value upon solvation, and, as expected, the largest effect occurs at the asymmetric geometry.

Proton Transfer Reactions

J. Phys. Chem., Vol. 100, No. 38, 1996 15449

IV. Concluding Remarks

F(1)-‚‚‚HF(2). We can further define the ground state energy Eg as the energy that corresponds to the frozen (FHF)- system. With the FDFT 1, 2 , and Eg we can define the off-diagonal element H12 using the identity

The present work has demonstrated that the FDFT method can offer a practical way for the evaluation of free energy surfaces of chemical reactions in solution representing both the solute and the solvent on an ab initio level. This advance is due to two features of our treatment: (i) the representation of the system by the FDFT formulation, which allows us to treat very large systems on an ab initio level and (ii) the use of the EVB potential surface as a reference potential in the quantum mechanical umbrella sampling/FEP calculations, which allows us to properly sample the solvent configuration space and to evaluate correctly and efficiently the solvent contribution to the activation free energy. The present FDFT formulation can easily be upgraded to allow the frozen densities to change. This CDFT formulation can be achieved by performing the “freeze and thaw” cycle24 or by adding inductive effects and allowing charge transfer between the solvent molecules, using EVB or charge equilibration considerations. It is also possible to consider charge transfer between the solute and its neighboring solvent molecules by iteratively changing the identity of the active space.24 That is, it is possible to perform DFT calculations on selected solvent molecules by holding the densities of the solute and other solvent molecules fixed. Such a procedure (e.g., ref 24) will resemble the “divide and conquer” approach37 but has a formally more rigorous coupling between the active and inactive spaces, since our nonadditive kinetic energy functional is formally rigorous. It is also important to mention our philosophy and aims are quite different than those expressed by most “divide and conquer” practitioners. We are concerned with chemical reactions in solution and in enzymes, and thus have no interest in the evaluation of the absolute energy of a molecule 15 Å from the reactive region (since it is inconceivable that a change here is related to the problem we address). Focusing on the evaluation of the exact energy of the entire system rather than its change during the reaction is probably an unnecessary obsession. The key ideal of FDFT, and its predecessor the pseudopotential model,21 is the replacement of hybrid QM/MM solvation treatments by more quantum mechanical approaches, rather than the introduction of a faster way of performing the unnecessary and presently impossible task of evaluating quantum mechanically the absolute free energy of the entire solutesolvent system. It is useful to note that the FDFT approach offers what is probably the most effective way of resolving the “connectingatom” problem in hybrid QM/MM approaches. In such methods one has to construct a connection region that reflects correctly the quantum chemistry of the first neighbor classical atoms. Although the use of hybrid orbitals2 or natural orbitals38 can reduce the problem, frozen densities can solve the problem, particularly when the frozen density modifications are included in the model. Another promising aspect of the FDFT method is the ability to explore the effect of the solvent on diabatic matrix elements such as the EVB coupling term H12 and the diabatic energies 1 and 2 from a first principles approach. The selection of the diabatic states is always a mathematical convenience rather than a reflection of any absolute reality, and we may chose the 1 and 2 that provide the most practical utility. For example, it is very convenient to define 1 as the energy of the density F1(r1°) obtained at the reactant equilibrium configuration F(1)H‚‚‚F(2)-. Similarly we may define 2 as the energy of F2(r2°) obtained at the product equilibrium configuration

[

]

Eg - 1 H12 H12 Eg - 2 ) 0

(17)

H12 ) x(EFDFT - 1(F1))(EFDFT - 2(F2))

(18)

which leads to

This approach is currently being studied in our lab. We expect to find that the average value of H12 will have a similar value in the gas phase and in solution. If this is found to be true, then one may obtain a unique method of obtaining the EVB diagonal energies by fitting them to 1(F1) and 2(F2). Acknowledgment. This work was supported by NIH Grant GM24492 and DOE Grant DE-FG03-94ER61945. The authors thank Professor D. R. Salahub for providing a copy of the deMon program. References and Notes (1) Warshel, A. Computer modeling of chemical reactions in enzymes and solutions; John Wiley: New York, 1991. (2) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (3) Warshel, A. J. Phys. Chem. 1979, 83, 1640. (4) Warshel, A. J. Phys. Chem. 1982, 86, 2218. (5) Muller, R. P.; Warshel, A. J. Phys. Chem. 1995, 99, 17516. (6) Luzhkov, V.; Warshel, A. J. Comput. Chem. 1992, 13, 199. (7) Warshel, A.; Weiss, R. M. J. Am. Chem. Soc. 1980, 102, 6218. (8) Hwang, J.-K.; King, G.; Creighton, S.; Warshel, A. J. Am. Chem. Soc. 1988, 110, 5297. (9) Aqvist, J.; Warshel, A. Chem. ReV. 1993, 93, 2523. (10) Kollman, P. A.; Hayes, D. M. J. Am. Chem. Soc. 1981, 103, 2955. (11) Rinaldi, D.; Rivail, J.-L. Theor. Chim. Acta. 1973, 32, 57. (12) Tunon, I.; Martins-Costa, M. T. C.; Millot, C.; Ruiz-Lopez, M. F. Chem. Phys. Lett. 1995, 241, 450. (13) Tapia, O.; Goscinski, O. Mol. Phys. 1975, 29, 1653. (14) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (15) Fortunelli, A.; Tomasi, J. Chem. Phys. Lett. 1994, 231, 34. (16) Longo, E.; Stamato, F. M.; Ferreira, R.; Tapia, O. J. Theor. Biol. 1988, 112, 783. (17) van Duijnen, P.; Tole, Th.; Broer, B.; Nieuwport, R. Int. J. Quantum Chem. 1980, 17, 651. (18) Bash, P. A.; Field, M. J.; Karplus, M. J. Am. Chem. Soc. 1987, 109, 8092. (19) Gao, J.; Luque, F. J.; Orozco, M. J. Chem. Phys. 1993, 98, 2975. (20) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991, 113, 8305. (21) Vaidehi, N.; Wesolowski, T. A.; Warshel, A. J. Chem. Phys. 1992, 97, 4264. (22) Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050. (23) Wesolowski, T; Warshel, A. J. Phys. Chem. 1994, 98, 5183. (24) Wesolowski, T. A.; Weber, J. Chem. Phys. Lett. 1996, 248, 71. (25) Stefanovich, E. V.; Truong, T. N. J. Chem. Phys. 1996, 104, 2946. (26) Cortona, P. Phys. ReV. B. 1991, 44, 8454. (27) Kim, H. J.; Hynes, J. T. Int. J. Quantum Chem. Symp. 1990, 24, 821. (28) Wei, D.; Salahub, D. R. J. Chem. Phys. 1994, 101, 7633. (29) Kohn, W.; Sham, L. J. Phys. ReV. 1965, 140, A1133. (30) Yamashita, K.; Morokuma, K.; Leforestier, C. J. Chem. Phys. 1993, 99, 8848. (31) Whitten, L. J. Chem. Phys. 1966, 44, 359. (32) St.-Amant, A.; Salahub, D. R. Chem. Phys. Lett. 1990, 169, 387. (33) Lee, F. S.; Chu, Z. T.; Warshel, A. J. Comput. Chem. 1993, 14, 161. (34) Hwang, J.-K.; Warshel, A. J. Am. Chem. Soc. 1987, 109, 715. (35) King, G.; Warshel, A. J. Chem. Phys. 1990, 93, 8682. (36) Valleau, J. P.; Torrie, G. M. In Modern Theoretical Chemistry; B. J. Berne, Ed.; Plenum: New York, 1977; Vol. 5, p 169. (37) Lee, C.; Yang, W. J. Chem. Phys. 1992, 96, 2408. (38) Thery, V.; Rinaldi, D.; Rivail, J.-L.; Maigret, B.; Ferenczy, G. G. J. Comput. Chem. 1994, 15, 269.

JP961068X