19570
J. Phys. Chem. B 2006, 110, 19570-19574
Using the Constrained DFT Approach in Generating Diabatic Surfaces and Off Diagonal Empirical Valence Bond Terms for Modeling Reactions in Condensed Phases Gongyi Hong, Edina Rosta, and Arieh Warshel* Department of Chemistry, UniVersity of Southern California, 3620 S. McClintock AVe., Los Angeles, California 90089-1062 ReceiVed: April 25, 2006; In Final Form: July 21, 2006
The empirical valence bond (EVB) model provides an extremely powerful way for modeling and analyzing chemical reactions in solutions and proteins. However, this model is based on the unverified assumption that the off diagonal elements of the EVB Hamiltonian do not change significantly upon transfer of the reacting system from one phase to another. This ad hoc assumption has been rationalized by its consistency with empirically observed linear free energy relationships, as well as by other qualitative considerations. Nevertheless, this assumption has not been rigorously established. The present work explores the validity of the above EVB key assumption by a rigorous numerical approach. This is done by exploiting the ability of the frozen density functional theory (FDFT) and the constrained density functional theory (CDFT) models to generate convenient diabatic states for QM/MM treatments, and thus to examine the relationship between the diabatic and adiabatic surfaces, as well as the corresponding effective off diagonal elements. It is found that, at least for the test case of SN2 reactions, the off diagonal element does not change significantly upon moving from the gas phase to solutions and thus the EVB assumption is valid and extremely useful.
Introduction Modeling reactions in solutions and in biological systems has remained a major challenge in computational chemistry (e.g., refs 1-3). Although significant progress has been accomplished by the use of macroscopic and/or simplified solvent models (e.g., refs 2-4), it is clear that microscopic approaches are crucial for further advances in the field. One of the most effective microscopic approaches has been the EVB method.5 However, this model involves the intuitive assumption that the off diagonal mixing elements (the Hijs) are similar in the gas phase and in solutions. This assumption remains untested and considered to be an ad hoc assumption. Before examining the validity of the assumption that the Hijs are transferable, it is important to clarify that this assumption is related to the idea that linear free energy relationships (LFERs) are valid and can be used with the same correlation coefficient in different environments. Now, the validity of LFERs for chemical processes in solutions has been the subject of many studies (e.g., refs 1, 6-10). Furthermore, experimental11,12 and theoretical13,14 studies have indicated that such relationships are also valid in enzymes, and more importantly, for the purpose of this work, can be used with similar correlation coefficients in the enzyme and in the reference solution reactions (e.g., ref 1). Thus, there are good qualitative reasons to assume that the off diagonal elements are transferable between different phases, but quantitative examinations of this issue are clearly long overdue. The present work describes an attempt to quantify the above hypothesis. In doing so, we face the problem that the transferability of the off diagonal terms cannot probably be proven in a rigorous way. Furthermore, we also realize that there is no current way of rigorously formulating the solvent dependence of the off diagonal term, despite some attempts in this * Corresponding author phone: 213-740-4114/7036; fax: 213-740-2701; e-mail:
[email protected].
direction.15 Fortunately however, it is clearly possible to try to validate the EVB assumptions by using a numerical approach, and such a strategy will be adopted here. In trying to establish the validity of the EVB assumptions numerically, we will have to choose a specific zero-order diabatic representation. This can be done by using different quantum mechanical approaches, but doing so by the frozen DFT (FDFT) or the constrained DFT (CDFT) methods16,17 is probably the most natural way (and currently perhaps the most rigorous18) of representing interacting fragments with full flexibility of allowing or restricting charge transfer between these fragments. Thus these approaches provide an effective way for examining the validity of the key EVB assumption (although alternative ways19 are also useful). That is, there are, of course different ways for obtaining H12 theoretically19,20 and estimating the effective H12 experimentally.21-23 However, the EVB intentionally avoids the rigorous VB formulation, since such a treatment suffers from major complications due to the overlap between the diabatic wave functions. Instead, the EVB chooses arbitrary diabatic states and then forces them (by using the proper H12) to reproduce “reality”, namely the adiabatic surface and charges. Of course, the selection of the diabatic states is not unique and rather arbitrary (except the requirement that the properties of the separate fragments at infinity will be reproduced by this diabatic state). However, taking the energies and charges of these states as those of the corresponding FDFT states provides a well-defined diabatic selection, without the need to specify a wave function and the resulting formal complications. Thus we have a formulation that can be examined rigorously within its framework. Using this formulation, we will try to show that the resulting H12 is transferable and the basic EVB assumption is, in fact, valid. Methods Our starting point is the generic reaction:
10.1021/jp0625199 CCC: $33.50 © 2006 American Chemical Society Published on Web 09/07/2006
DFT and Off Diagonal Empirical Bond Terms
A+B-CfA-B+C
J. Phys. Chem. B, Vol. 110, No. 39, 2006 19571
(1)
Where we can describe the reactant (r) and product (p) states by diabatic valence bond like wave functions:
The total energy of an N-electron system is defined within the Kohn-Sham formulation as follows:
E[F] ) T[F] +
Ψr ) Ψ [A] Ψ [B - C] Ψp ) Ψ [A - B] Ψ [C]
(2)
We can also consider the adiabatic ground state wave function, which describes the entire system. Comparing the diabatic and adiabatic energies of our system will allow us to examine the validity of the EVB approach. This point will be clarified below by briefly considering the EVB and CDFT/QM/MM methods. The EVB approach describes reactions by using diabatic surfaces and mixing them by effective off diagonal terms. In the simple 2-state case, e.g., i ) 1 or 2, we can write the following:
Hii ) i ) Rigas + Uiintra + UiSs + Uss
(3)
Hij ) A exp{- µ(R - R0)}
(4)
where Uiintra is the intramolecular potential of the solute system, which is represented by a molecular mechanics-type force field that describes the current bonding structure of the i state i. USs represents the interaction energies (van der Waals and electrostatic interactions) between the solute (S) and the environment (s), Uss is the force field of the solvent system. Finally Rigas is a constant that reflects the difference between the gas-phase energy of the diabatic states when all the corresponding fragments are at infinite separation.1 The off diagonal element H12 depends typically on the distance R between the A and C atoms. The ground-state energy, Eg of the system is obtained as the lowest eigenvalue of the 2 × 2 EVB Hamiltonian.
|
Eg - 1 H12 H12 Eg-2
|
(5)
Now the EVB approach assumes that H12 does not change significantly upon moving the solute from the gas phase to solutions and our task is to examine this approximation in a rigorous way. This can be done by using the CDFT approach to generate diabatic surfaces (see below) and then extracting the relevant Hij by relating the diabatic surfaces to the corresponding adiabatic surface (obtained by the full DFT method). Performing such calculations in the gas phase and in solution should allow us to extract the effect of the solvent on Hij and thus to validate our assumption. At any rate, once we obtain the quantum mechanical values of the diabatic 1 and 2, as well as the adiabatic Eg, we can obtain H12 through eq 5, as follows:1
H12 ) x(Eg - 1)(Eg - 2)
(6)
This constitutes our operational definition of H12. With the definition of eq 6, we can look for an ab initio treatment that can produce both the diabatic and adiabatic states, where the relationship between these states is defined the same way in the gas phase and in solution. One of the most unique ways of doing so is provided by the FDFT/CFDF method,17 and since the details on the frozen density functional method can be found elsewhere, we will provide below only a brief description.
∫Vext(r)F(r)dr + F(r)F(r′) 1 ∫ ∫ |r - r′| drdr′ + Exc[F(r)] 2
(7)
where Vext(r) is the external potential, Exc[F(r)] is a density functional of the exchange-correlation energy and N
F(r) )
∑ φi*(r)φi(r)
(8)
( )
(9)
i)1 N
T[F(r)] )
1
∑ ∫φi*(r) -2 ∇2 φi(r)dr i)1
In the FDFT/CDFT approach, we divide the system into subsystems, and apply Kohn-Sham formulation to each subsystem separately, taking into account the interaction with the other subsystems. For instance, if we divide the whole system into two subsystems with NI, NI ′ electrons, respectively, the total energy of the system is defined as follows:
ECDFT[FI, FI′] ) T[FI] + T[FI′] + Tsnadd[FI, FI′] + F(r)F(r′) 1 drdr′ + Vext(r)F(r)dr + Exc[F(r)] (10) 2 |r - r′|
∫∫
∫
∫FI(r)dr ) NI
(11)
∫FI′(r)dr ) NI′
(12)
F(r) ) FI(r) + FI′(r)
(13)
[FI,FI′] ) Ts[FI + FI′] - Ts[FI] - Ts[FI′] Tnadd s
(14)
where Ts[F] is defined as the orbital-independent kinetic energy functional of the electron density of the system.24 Fixing the electron density of region I′ and minimizing the total energy with respect to the orbitals of region I leads to the following modified Kohn-Sham equation: # (r)]φi(r) ) iφi(r) [-1/2∇2 + Veff,I
(15)
# where the effective potential Veff,I (r) is defined as follows:
Veff,I(r) ) VKS eff (r) +
δTnadd[FI(r), FI′(r)] δFI(r)
(16)
Here, the Kohn-Sham potential, VKS eff (r) given by the following:
VKS eff (r) ) Vext(r) +
F (r′)
F (r′)
∫ |r I- r′| dr′ + ∫ |r I′- r′| dr′ + Vxc[F(r)]
(17)
Once FI(r) is known, a counterpart equation for region I′ could be derived. These two coupled equations can be solved iteratively. In case a part of the system is described classically we use a QM(FDFT)-MM formulation where the effect of the MM part is incorporated in the FDFT Hamiltonian as described in ref 25.
19572 J. Phys. Chem. B, Vol. 110, No. 39, 2006
Hong et al.
The FDFT approach can be used to define the diabatic energies of the reactant state (H11) and product state (H22) of our system. This can be done by decomposing the density as either F1(r) ) FNA(r) + FNBC(r) (for one diabatic state) or F2(r) ) FNAB(r) + FNC(r) (for the other diabatic state). The specific treatment will be described in the next section. Since we are interested here in the validity of the EVB in studies of reactions in solutions, we will have to evaluate the free energy functional of the FDFT surface by the proper mapping procedure. The general approch was described in many of our works but will be reviewed here in order to give the reader a complete picture of our strategy. We start the evaluation of the activation free energy by running molecular dynamic simulation on a mapping potential of the form
Um ) (1 - λm)H11 + λmH22
(19)
where m designates an average of the quantity defined inside the bracket over the configurations generated by propogating trajectories over Um. The overall free energy profile is given by: m)n-1
∆G(λn) )
∑
m)0
δG(λm f λm+1)
(20)
The actual adiabatic ground-state free energy profile along the reaction coordinate, x ) 1 - 2 ) H22 - H11, is evaluated by the following1:
∆g(x′) ) ∆Gm - β-1 ln 〈δ(x - x′) exp[-β(Eg(x) - m(x))]〉m (21) Similar formulas are used to evaluate the diabatic free energies associates with H11 and H22. Results and Discussion In the present case, we consider as a general benchmark for charge-transfer reactions, the SN2 reaction:
ClA- + CH3ClB f ClACH3 + ClB
(22)
The VB resonance structures for the diabatic reactant (r) and product (p) states are, respectively
Ψr ) [ClA
CH3 - ClB]
Ψp ) [ClA - CH3
ClB]
[FCl-A, FCH3ClB] + ) T[FCl-A] + T[FCH3ClB] + Tnadd s 1 2
(FCl-A(r) + FCH3ClB(r))(FCl-A(r′) + FCH3ClB(r′))
∫∫ ∫Vext(r)(FCl
(23)
To represent the VB diabatic reactant state by the FDFT approach, we divide the whole system into the Cl-A ion and the CH3ClB subsystems, and express the gas-phase energy of this diabatic state as follows:
|r - r′|
-
drdr′+
(r) + FCH3ClB(r))dr + Exc[FCl-A(r) +
A
FCH3ClB(r))] (24) The modified Kohn-Sham equation (eq 15) was applied to obtain the densities of the subsystems. Similarly, we described the product state by the following:
H22 ) 2 ) EFDFT[FCl-B, FCH3ClA] [FCl-B, FCH3ClA]+ ) T[FCl-B] + T[FCH3ClA] + Tnadd s
(18)
Here λm gradually increases from 0.0 to 1.0 with a small constant increment and the system is driven gradually from reactant state to product state. The free energy contribution associated with the change of Um to Um′ (δG[λm f λm′]) can be obtained by using the free energy perturbation (FEP) method:
exp{-δG[λm f λm′]β} ) 〈exp{- (Um′ - Um)β}〉m
H11 ) 1 ) EFDFT[FCl-A, FCH3ClB]
1 2
∫∫
(FCl-B(r) + FCH3ClA(r))(FCl-B(r′) + FCH3ClA(r′)) +
∫Vext(r)(FCl
|r-r′| -
drdr′
(r) + FCH3ClA(r))dr +
B
Exc[FCl-B(r) + FCH3ClA(r))] (25) , Finally, we obtain the adiabatic ground-state energy Efull-dft g by treating the whole system at the same level with a regular DFT calculation. The same treatment is also implemented in solution, where the solvent effect is simulated by a QM/MM approach (e.g., ref 26), using the formulation of ref 25. The actual QM/MM treatment has been implemented in the program MOLARIS27,28 by using a very convenient “weak coupling” philosophy where any QM package can be called by MOLARIS. In this work, the QM region is treated by the DFT and CDFT/FDFT method that were implemented in the deMon program.29 The Becke88 exchange potential,30 the Perdew86 correlation potential,31 and the nonadditive kinetic energy functional of eq 14 were used in all the calculations. Standard 6-31+G* basis sets were used for expanding the molecular orbitals. The MM part (the water molecules in the classical region) was represented by the regular ENZYMIX force field.28 The simulation model includes an 18 Å explicit simulation sphere, surrounded by a surface region whose average polarization and radial distribution are determined by the surface constrained allatom solvent (SCAAS) model.32 The surface region is embedded in a bulk continuum region with a dielectric constant of 80. Long-range electrostatic effects are treated by the local reaction field (LRF) model.33 The diabatic QM(FDFT)/MM and adiabatic QM(DFT)/MM free energy profiles were obtained by using 21 windows of 1 ps each with 0.5 fs stepsize at 300°K. The runs were preequilibrated in each window. The calculated profiles are summarized in Figure 1 and show the expected trend in SN2 reactions,34 where the adiabatic barrier is much larger in solution due to the loss of solvation upon formation of the Cl-1/2‚‚‚CH3‚‚‚Cl-1/2 delocalized transition state. We also examined free energy calculations with longer MM relaxations, and obtained similar results, thus estabilished a reasonable convergence. The adiabatic barrier for our SN2 reaction in gas phase in known to be overestimated by about 9 kcal/mol due to the use of the DFTBP method.35 Assuming similar effect in solution, we expect that the barriers would be around 16 + 9 ) 25 kcal/mol with a better ab initio method. Now, since the calculated profile starts from configurations where the reactant molecules are in the same solvent cage, the calculated barrier should be increased by about
DFT and Off Diagonal Empirical Bond Terms
Figure 1. Diabatic and adiabatic FDFT energy profiles for the reaction, Cl- + CH3Cl f ClCH3 + Cl-, in the gas phase and in solution. The reaction coordinate is defined as the energy difference between the diabatic surfaces, ∆ ) 1 - 2. Note that the adiabatic barrier is underestimated by about 9 kcal/mol due to the use of the DFT-BP method.35 Also note that the gas-phase profile is evaluated using the same structures that are obtained in the simulations of the solution reaction.
J. Phys. Chem. B, Vol. 110, No. 39, 2006 19573
Figure 2. Plot of Hij of the reaction considered in Figure 1 in the gas phase and in solution. The data are obtained from the diabatic and the adiabatic curves of Figure 1, using H12(∆) ) x(1-Eg)(2-Eg).
despite the very large change in the corresponding diagonal energies. This suggests that the EVB off diagonal elements are indeed transferable between different phases. Concluding Remarks
2.3 kcal/mol (e.g., see ref 36) and be around 27.3 kcal/mol. This result is in a good agreement with the experimental barrier of about 26.6 kcal/mol.37 It is important to note that our calculations provide a complete free energy profile for an ab initio QM/MM simulation of a reaction in solution with a complete flexibility of the solute and the solvent. Such studies are extremely challenging and quite rare (see discussion in ref 38), even gas-phase free-energy calculations on this relatively simple system have significant difficulties.35 Also note that our free energy profile reflects the non equilibrium solvation effects (which are automatically obtained by eq 21) and the solute entropic effects, which are typically missing in PMF type calculations with fixed solute coordinates, for each mapping window (see discussion in ref 39). It is important to note that the gas-phase barrier is evaluated for the same structures as those used in the solution study. This means that the barrier is evaluated more or less from the minimum of the gas-phase profile (rather than the configuration of the separated fragments), where the corresponding gas-phase barrier is expected to be around 4 kcal/mol in the DFT-BP method,35 as compared to our calculated barrier of about 6 kcal/ mol. This 2 kcal/mol deviation reflects the fact that the gasphase maximum is calculated using the structures of the TS in solution. Inserting the free energy curves of Figure 1 into eq 6 allows us to obtain the effective H12 for both the gas phase and solution reaction for any point along the reaction coordinate between ∆ ) -60 and ∆ ) 60 kcal/mol. At larger absolute values of ∆ the calculations do not seem to be sufficiently reliable in view of the asymmetry of the results for positive and negative ∆, thus we did not evaluate H12 in this range. Obtaining fully symmetric results would require significantly longer runs and sampling with different initial conditions. The results obtained from the above analysis are summarized in Figure 2. The most striking conclusion from Figure 2 is that H12 is very similar in the gas phase and the solution cases,
This work explores the long-standing question of the validity of the EVB's key assumption that the off diagonal mixing terms are transferable between different environments. This examination exploited the ability of the FDFT approach to generate convenient and physically consistent diabatic states for QM/ MM treatments and thus to be used in numerical studies of the relationship between the diabatic and adiabatic surfaces in different environments. Obtaining these surfaces allows us to evaluate the corresponding off diagonal elements in different environments and thus to examine the effect of the environment on these crucial parameters. It has been found that at least for the test case of SN2 reactions the off diagonal element does not change significantly upon moving from gas phase to solution. This is quite significant since the SN2 reaction is a prototype of many charge transfer reactions and since the present study found that the off diagonal element does not depend significantly on the environment for the entire reaction coordinate and thus for a different degree of charge transfer. Of course, more studies will be required in the future; examining for example charge separation reactions, but the present study is the first one that actually examines the crucial problem addressed in this work. Now from the indications obtained so far, it seems that the EVB assumption is valid and extremely useful. The transferability of the off diagonal elements between different phases is important not only for quantitative EVB studies but also as a justification of the assumptions made for a long time by the physical organic chemistry community in LEFR studies, where it has been assumed implicitly that the correlation coefficients do not change significantly in different environments. Finally, it is important to note that even if in some cases, the off diagonal term changes upon moving from the gas phase to solution are larger than the current estimate, the changes upon moving from solutions to enzymes are expected to be significantly smaller since in both cases we deal with polar environments.
19574 J. Phys. Chem. B, Vol. 110, No. 39, 2006 Acknowledgment. This work was supported by the National Institutes of Health (NIH) grant GM 24492. We thank the High Performance Computing Center (HPCC) at the University of Southern California (USC) for computer time. References and Notes (1) Warshel, A. Computer Modeling of Chemical Reactions in Enzymes and Solutions; John Wiley & Sons: New York, 1991. (2) Cramer, C. J.; Truhlar, D. G. Structure and ReactiVity in Aqueous Solution; American Chemical Society: Washington, 1994; Vol. 568. (3) Tomasi, J.; Persico, M. Chem. ReV. 1994, 94, 2027. (4) Honig, B.; Sharp, K.; Yang, A.-S. J. Phys. Chem. 1993, 97, 1101. (5) Warshel, A.; Weiss, R. M. J. Am. Chem. Soc. 1980, 102, 6218. (6) Leffler, J. E.; Grunwald, E. Rates and Equilibria of Organic Reactions; Wiley & Sons: New York, 1963. (7) Marcus, R. A. Annu. ReV. Phys. Chem. 1964, 15, 155. (8) Albery, W. J. Annu. ReV. Phys. Chem. 1980, 31, 227. (9) Kreevoy, M. M.; Kotchevar, A. T. J. Am. Chem. Soc. 1990, 112, 3579. (10) Schweins, T.; Geyer, M.; Kalbitzer, H. R.; Wittinghofer, A.; Warshel, A. Biochemistry 1996, 35, 14225. (11) Fersht, A. R.; Leatherbarrow, R. J.; Wells, T. N. C. Nature (London) 1986, 322, 284. (12) Toney, M. D.; Kirsch, J. F. Science 1989, 243, 1485. (13) Warshel, A. Simulating the Energetics and Dynamics of Enzymatic Reactions; Pontificiae Academiae Scientiarum Scripta Varia, 1984; Vol. 55. (14) Warshel, A.; Schweins, T.; Fothergill, M. J. Am. Chem. Soc. 1994, 116, 8437. (15) Schmitt, U. W.; Voth, G. A. J. Phys. Chem. B 1998, 102, 5547. (16) Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1993, 97, 8050. (17) Hong, G.; Sˇ trajbl, M.; Wesolowski, T. A.; Warshel, A. J. Comput. Chem. 2000, 21, 1554.
Hong et al. (18) Wesolowski, T. A. CHIMIA 2002, 56, 707. (19) Mo, Y.; Gao, J. J. Chem. Phys. 2000, 104, 3012. (20) Shaik, S.; Shurki, A. Angew. Chem., Int. Ed. Eng. 1999, 38, 586. (21) Hush, N. S. Trans. Faraday Soc. 1961, 57, 557. (22) Hush, N. S. Electrochim. Acta 1968, 13, 1005. (23) Gwaltney, S. R.; Rosokha, S. V.; Head-Gordon, M.; Kochi, J. K. J. Am. Chem. Soc. 2003, 125, 3273. (24) Wesolowski, T. A. J. Chem. Phys. 1997, 106, 8516. (25) Olsson, M. H. M.; Hong, G.; Warshel, A. J. Am. Chem. Soc. 2003, 125, 5025. (26) Warshel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (27) Chu, Z. T.; Villa, J.; Strajbl, M.; Schutz, C. N.; Shurki, A.; Warshel, A. MOLARIS version beta9.05. University of Southern California, LosAngeles, 2004. (28) Lee, F. S.; Chu, Z. T.; Warshel, A. J. Comput. Chem. 1993, 14, 161. (29) St-Amant, A.; Salahub, D. R. Chem. Phys. Lett. 1990, 169, 387. (30) Becke, A. D. Phys. ReV. A 1988, 38, 3098. (31) Perdew, J. P. Phys. ReV. B 1986, 33, 8822. (32) Warshel, A.; King, G. Chem. Phys. Lett. 1985, 121, 124. (33) King, G.; Warshel, A. J. Chem. Phys. 1990, 93, 8682. (34) Shaik, S. S.; Schlegel, H. B.; Wolfe, S. Theoretical Aspects of Physical Organic Chemistry. Application to the SN2 Transition State; WileyInterscience: New York, 1992. (35) Yang, S. Y.; Fleurat-Lessard, P.; Hristov, I.; Ziegler, T. J. Phys. Chem. A 2004, 108, 9461. (36) Sˇ trajbl, M.; Floria´n, J.; Warshel, A. J. Am. Chem. Soc. 2000, 122, 5354. (37) Albery, J. K.; M. M. AdV. Phys. Org. Chem. 1978, 16, 87. (38) Rosta, E.; Klahn, M.; Warshel, A. J. Phys. Chem. B 2006, 110, 2934. (39) Shurki, A.; Warshel, A. AdV. Protein Chem. 2003, 66, 249.
