J. Phys. Chem. B 2008, 112, 12491–12497
12491
VB/MMsThe Validity of the Underlying Approximations Avital Sharir-Ivry and Avital Shurki†,* Department of Medicinal Chemistry and Natural Products, The Lise Meitner-MinerVa Center for Computational Quantum Chemistry, School of Pharmacy, The Hebrew UniVersity of Jerusalem, Jerusalem 91120, Israel ReceiVed: March 27, 2008; ReVised Manuscript ReceiVed: May 12, 2008
Two recently developed methods VB/MM and density embedded VB/MM (DE-VB/MM) are compared, and their respective approximations are examined. The two methods combine valence-bond (VB) calculations with molecular mechanics (MM) and aim to allow VB analysis of reactions in large biological environments. Furthermore, the two methods utilize two major approximationssregarding both the overlap and the reduced resonance integral between the various VB configurations. The difference between the two methods, however, is that VB/MM employs these approximations for the overall interaction of the reacting fragments with their surrounding, whereas DE-VB/MM employs the approximations only with regards to the van der Waals (VdW) interactions whereas the electrostatic interactions are calculated rigorously at the quantum level. The approximations that lay the grounds for the two methods involve the assumption that the overlap between the VB configurations and the respective reduced resonance integral are both invariant to the environment. Similar approximations are utilized in several other VB-based QM/MM methods. However, although extensively used, these approximations were never rigorously proved. Here, we exploit the development of the DE-VB/ MM method to numerically examine the approximations by calculating the accurate as well as the approximated values of overlap and reduced resonance integral for systems where the environment involves only electrostatic interactions. The quality of the approximations is examined together with their effect on the absolute energies, the wave function, and the overall energetics. Three test cases are chosen, the dissociation of CH3F and LiF and the identity SN2 reaction. It is shown that the approximations are usually good with the exception of cases where extreme changes are expected in the wave function. Furthermore, the impact of the approximations on the overall wave function and the overall energetics is found to be quite small. It is concluded that VB/ MM, where the approximations are used more extensively, can serve as the first method of choice. Introduction Chemical reactions in biological systems such as solvents and proteins comprise the majority of reactions in nature and as such are of great interest. Thus, the development of new methods that would aid in understanding their reactivity patterns is important. Current methods that address these questions combine quantum calculations with molecular mechanics (e.g., refs 1–19). Many of these methods are based on valence bond (VB) ideas because of their advantage in shedding light on reactivity.16–41 Yet, usually, these methods are semiempirical, and parametrization, therefore, governs the quality of the results.16–35 Several attempts to use ab initio VB and avoid the need for parametrization were also suggested.38–43 These methods, however, have not yet been tested on protein environments or simply were not developed to do so. Recently, we have presented two new hybrid ab initio VB/ molecular mechanics (MM) methods: the VB/MM41 and the density-embedded VB/MM (DE-VB/MM).42 Both methods consider the effect of the environment on each one of the VB configurations separately. The overall adiabatic state is then obtained as a result of mixing of these VB configurations. Furthermore, both methods utilize two approximations in order to evaluate the off-diagonal elements of the VB matrix in the * Corresponding author. E-mail:
[email protected]. Phone: +9722-675-8696. Fax: +972-2-675-7076. † Affiliated with the David R. Bloom Center for Pharmacy at the Hebrew University.
new environment (solvent or protein). These include the overlap integral, Sij, between VB configurations φi and φj and their respective reduced resonance integral, βij. Thus, both methods assume Sij and βij to be invariant with respect to the environment. These assumptions are based on similar assumptions that have laid the foundations of various established methods, such as the empirical VB (EVB).17,18,21 In these methods, Sij, the overlap between the diabatic states, is considered to be zero, and Hij, the resonance integral, is assumed to be invariant with respect to the environment. The reduced resonance integral βij is identical to Hij in the absence of overlap, thus, leading us to deduce our above-described assumptions with respect to Sij and βij. The main difference between VB/MM and DE-VB/MM is rooted in the manner these methods consider the electrostatic interactions of the reacting fragments with the environment. VB/ MM employs the point charge (mechanical) embedding scheme,44–46 thus calculating classically the electrostatic interactions of the environment with each VB diabatic state separately.41 DE-VB/MM, on the other hand, employs the more sophisticated electrostatic (density) embedding scheme,44,46,47 by incorporating the charges of the environment in the ab initio quantum Hamiltonian.42 Therefore, DE-VB/MM utilizes the above-described approximations only with regards to the van der Waals (VdW) interactions and their effect on the overall energy. The reliability of both methods depends largely on the validity of the assumptions used in the approximations. The invariance
10.1021/jp802667y CCC: $40.75 2008 American Chemical Society Published on Web 09/04/2008
12492 J. Phys. Chem. B, Vol. 112, No. 39, 2008
Sharir-Ivry and Shurki
of neither the overlap nor the reduced resonance integral has been proven. Recently, a numerical evaluation of the value of Hij by utilizing the FDFT method for the identity SN2 reaction demonstrated that the gas-phase value of Hij serves as a good approximation for its corresponding value in solution.48 The relationship between Hij and both Sij and βij suggests that the approximation may also be good for the latter two (Sij and βij). Yet, the deduction regarding the invariance of Sij and βij to the environment based on the invariance of Hij is not straightforward. Furthermore, the approximation regarding Hij was shown to be good in only one example. Therefore, examination of the approximation regarding both Sij and βij is still required. The present work presents a detailed examination of the validity of the approximations regarding both the overlap and the reduced resonance integral (Sij and βij). Utilizing the DEVB/MM approach while disregarding the VdW interactions provided us with accurate values of both Sij and βij in an environment of point charges. This, in turn, allowed us to compare these values with the original gas-phase values and examine the approximation when only electrostatic interactions are involved. Three test cases are considered here, to cover a variety of situations: dissociation of LiF and CH3F and the identity SN2 reaction of chloride anion with methyl chloride. It is shown that the approximations are usually very good with an exception of cases where extreme changes take place in the wave function. The first part of the paper describes the basic principles behind the two methods (VB/MM and DE-VB/MM). The approximations are then examined, and their effect on the energies is shown. The last part of the paper compares the results obtained by the two methods, VB/MM and DE-VB/MM, with emphasis on the effect of the environment on both the wave function and the overall energetics. It is concluded that the simpler embedding scheme used in VB/MM captures most of the effects and, therefore, can serve as a prime method to study chemical reactivity. Theory and Methodology Both the VB/MM and the DE-VB/MM methods were discussed in detail elsewhere.41,42 Thus, here, the two methods are described briefly with emphasis on the various approximations and the way they are utilized. Any VB calculation involves several diabatic VB configurations, φi, that span the system’s wave function. For a system with N diabatic VB configurations, the resulting energy is approximated by the lowest eigenvalue of the following N secular equations N
∑ (Hij - εSij)cj ) 0
(1)
j
where Hij and Sij are the Hamiltonian and overlap matrix elements between diabatic configurations φi and φj and cj are the respective coefficients of the diabatic configurations φj in the overall adiabatic wave function. VB/MM and DE-VB/MM are designed to study reactivity in large systems. Therefore, the calculations account for both the reacting fragments and the effect of the environment (solution/protein). The diagonal elements, Hii, signify the energy of the respective diabatic configurations, (φi in this case). Their interaction (electrostatic and VdW) with the environment is simply added to their gas-phase value in a form which depends on the method. At the VB/MM level, both electrostatic and VdW interactions are calculated classically. In the DE-VB/MM,
however, only the VdW interaction is treated classically, whereas the electrostatic interaction is obtained by incorporating a one-electron integral term which represents the electric field caused by the environment’s charges, into the core Hamiltonian. Thus, we note that in DE-VB/MM, the energy obtained from the quantum Hamiltonian already accounts also for the electrostatic interaction with the environment. The basic approximation that lays the grounds for the two methods is that both the overlap integral Sij and the reduced resonance integral βij are invariant to the environment. Namely,
Sijenv ) Sijg
(2)
βijenv ) βijg
(3)
g env where Sgij, Senv ij , βij, and βij are the overlap and the reduced resonance integrals in gas phase and in the new environment, and βij is defined by
1 βij ) Hij - (Hii + Hjj)Sij 2
(4)
where Hkk is the energy of φk and Hkl is the resonance integral between the two diabatic VB configurations φk and φl. We note in this respect that βij is identical to Hij in the absence of overlap (when Sij ) 0). Applying the approximation of the reduced resonance integral (eq 3) and expanding it by using eq 4, we obtain
Hijenv -
1 env (H + Hjjenv)Sijenv ) Hijg - 21 (Hiig + Hjjg)Sijg 2 ii
(5)
Implementation of the approximation of the overlap (eq 2) and reordering the terms of eq 5 result in
Hijenv ) Hijg +
1 env (H + Hjjenv - Hiig - Hjjg)Sijg 2 ii
(6)
The difference in energy of the ith diabatic VB configuration in the new environment compared to the gas phase defines the interaction of that configuration with its surrounding:
∆Hiiint ≡ Hiienv - Hiig
(7)
By using this definition, eq 6 can be written in a more compact way:
Hijenv ) Hijg +
1 (∆Hiiint + ∆Hjjint)Sijg 2
(8)
Accordingly, VB/MM approximates the resonance integral in the various environments by using eqs 7 and 8. In the DEVB/MM method, on the other hand, the effect of the electrostatic interaction with the environment on the energy of the diabatic states and on the resonance integral between them is already considered at the quantum level. Consequently, DE-VB/MM requires employing the approximation described above only with regards to the VdW interactions between the reacting fragments and the environment. Therefore, in this case, ∆Hint ii involves only the VdW interaction rather than the overall interactions as defined in eq 7 and will be marked accordingly, ∆HVdWint . ii Furthermore, Hgij and Sgij in eq 8 are replaced by the corresponding resonance and overlap integrals as obtained from the quantum g+ele g+ele Hamiltonian, HQM,ij and SQM,ij , which account also for the electrostatic interaction with the environment and do not
VB/MMsThe Validity of the Underlying Approximations
J. Phys. Chem. B, Vol. 112, No. 39, 2008 12493
SCHEME 1
represent only the gas phase values. Thus, eq 9 presents the resulting approximated Henv ij in the DE-VB/MM method: g+ele Hijenv ) HQM,ij +
1 g+ele (∆HiiVdWint + ∆HjjVdWint)SQM,ij 2
(9)
Finally, in both methods, the secular equations for the system are solved, and the classical energy of the surrounding is added to give the overall energy of the system. Computational Details The quantum calculations were performed with the GAUSSIAN 03 program package49 to obtain one-electron integrals and the Xiamen VB (XMVB) program package50,51 for the various VB calculations. The classical calculations including the dynamics and all the force-field-related calculations utilized the MOLARIS package.52,53 Additionally, an interface program that communicates between the two programs was used. Three different test cases were examined: the bond dissociation of CH3-F and Li-F and the identity SN2 reaction of methyl chloride with a chloride ion. The gas phase results involved the breathing orbitals VB (BOVB) level of calculation54,55 and were then used to get the solution results at the BOVB/MM and the DE-BOVB/MM levels. In the BOVB, each VB structure has its own set of orbitals.54,55 Thus, a minimal set of VB structures was chosen for each system, and the orbitals of each structure were optimized in the presence of the orbitals of the other structures. More specifically, two VB structures were considered for the dissociation reactions, the covalent configuration, φcov (drawings 1 and 3 in Scheme 1 for CH3F and LiF, respectively) and the chemically significant ionic configuration, φion (drawings 2 and 4 in Scheme 1). For the SN2 reaction, only three VB configurations were considered, the reactants covalent configuration, φRcov (drawing 5 in Scheme 1), the products covalent configuration, φPcov (drawing 6 in Scheme 1), and the tri-ionic (ionic hereafter) configuration, φion (drawing 7 in Scheme 1). The quantum calculations employed the 6-31G* basis set, and the geometries were optimized at the Hartree-Fock (HF) level of calculation. Throughout the BOVB calculations, innercore electrons of the various atoms (1s for lithium carbon and fluorine and 1s, 2s, and 2p for chlorine) were frozen at the HF level, whereas the remaining electrons were explicitly calculated in the VB scheme. The classical simulations involved all atom calculations. The system was spherical and involved three major regions. The first region involved the reacting fragments and was considered the quantum region. This region was treated at the EVB level during the dynamics and at the BOVB level for the actual energetics of chosen geometrical configurations. The second region involved the solvent molecules up to an 18 Å radius and was presented with the standard ENZYMIX force field,53 where the solvent molecules were subjected to the surfaceconstrained all-atom solvent boundary56 and the local reaction field model.57 The last region involved a 2Å shell of Langevin dipoles and a dielectric continuum to simulate the bulk solvation.
The dissociation/reaction coordinates were achieved by following the potential of mean force (PMF) protocol,58–61 fixing the coordinates of the reacting fragments while relaxing the rest of the system. Each point along the dissociation/reaction paths was first relaxed for 50 ps with fixed coordinates for the reacting fragments and with some initial guess for the wave function (normally taken from a gas-phase VB calculation). A singlepoint calculation by the method of interest (VB/MM or DEVB/MM) was then carried out, and the resulting new wave function was introduced as a guess for a new relaxation of the solvent. The process was repeated until changes in the wave function and in particular the weights of the various VB configurations were below 0.05. Next, each point along the dissociation curves of LiF and CH3F was simulated for 80 ps and for 50 ps along the SN2 reaction path. Geometry configurations were collected every 1 ps, and their respective energy was calculated at both the VB/MM and DE-VB/MM levels. Finally, the results were analyzed by using the PMF procedure combined with the free energy perturbation (FEP) theory,62,63 to obtain the free energy of the dissociation/reaction profiles. Examination of the approximation regarding both the overlap and the reduced resonance integral was carried on in the following manner: the gas-phase values Sgij and βgij calculated for the reacting fragments were compared with the correspondg+ele g+ele ing quantities, SQM,ij and βQM,ij calculated for the same reacting g+ele fragments in an electric field. The latter quantities, SQM,ij and g+ele βQM,ij , were taken from the quantum step of the DE-VB/MM calculation. The electric field, in that case, is produced by introducing the environment’s partial charges, as external point charges, in a one-electron integral to the quantum Hamiltonian. The partial charges were determined by the ENZYMIX force field.53 Finally, in order to study the effect of the approximations on the energetics, a comparison of the eigenvalues obtained with and without the approximations was carried out. Eigenvalues taken from the quantum step of the DE-VB/MM calculation were considered as the accurate ones, εacc, for a system that involves only electrostatic interactions between the reacting fragments and their environment. In order to obtain eigenvalues that emerge from the approximations and that will be comparable with those discussed above, Henv ii values were taken from the quantum step of the DE-VB/MM calculation, and the interaction energy was defined by using eq 7. We note again that the interaction energy in this case consists only of electrostatic interaction. Equation 8 was then used to calculate Henv ij . Eventually, the secular equations with the newly defined matrix elements were solved, and eigenvalues resulting from the use of the approximations, εapx, were obtained. The three different test cases will be discussed. However, in order to avoid a plethora of graphs, the corresponding graphs will be given only for the SN2 reaction in which there are three significant VB configurations to show that the approximations are good also in the more general case and not only in the specific case of 2 by 2. The results for both dissociation curves are given in the Supporting Information. Results and Discussion Examination of the Approximations. Examination of the approximations involved primarily the comparison of the overlap and the reduced resonance integrals of the reacting fragments calculated in gas phase, Sgij and βgij, with the corresponding quantities calculated for the reacting fragments in an electric g+ele g+ele g+ele g+ele field, SQM,ij and βQM,ij . The calculation of SQM,ij and βQM,ij does not involve these approximations, and thus, they were regarded
12494 J. Phys. Chem. B, Vol. 112, No. 39, 2008
Sharir-Ivry and Shurki
R Figure 1. (a) Overlap and (b) reduced resonance integral values between φcov , 5, and φion, 7, VB configurations of various geometries along the SN2 reaction. Black bars represent the overlap values in the gas phase, and green bars present the values in solution.
TABLE 1: Overlap and Reduced Resonance Integral Values between Various VB Configurations of the SN2 Reaction in Reactants’ Complex and TS Geometries (Comparison of Gas Phase with Solution Values) reactants gas 1 2 3
Scov R, cov P Scov R, ion Scov P, ion βcov R, cov P βcov R, ion βcov P, ion
4 5 6 a
solution
Overlap 0.131 0.123 -0.490 -0.493 -0.251 -0.236
TS gas
solution
0.126 -0.365 -0.364
0.118 -0.360 -0.359
Reduced Resonance Integrala -50.5 -51.5 -45.9 137.9 137.2 75.4 49.4 54.0 75.4
-45.1 75.6 75.4
Values are in kcal/mol.
as accurate values of overlap and reduced resonance integrals in an environment which includes only electrostatic interactions. Comparison of the overlap and the reduced resonance integral values in the gas phase with those obtained in solution is presented in Table 1 and Figure 1. Entries 1-3 in the table as well as Figure 1a correspond to the overlap values, and entries 4-6 and Figure 1b correspond to the reduced resonance integral. In the table, the three entries for each quantity account for the three possible combinations of the VB configurations. Moreover, the values are given for two different geometries along the reaction coordinate, the reactant’s complex, and the transition state (TS). Figure 1, on the other hand, presents overlap and reduced resonance integral values only between φRcov, 5, and φion, 7, for various geometries along the reaction coordinate.64 The overlap and reduced resonance integral values between the other configurations along the different geometries as well as the overlap and reduced resonance integral in the two dissociation cases are given in Figure 1S in the Supporting Information. By looking at the results displayed here and in the Supporting Information and comparing the overlap and the reduced resonance integral, similar overall trends are found, as expected. Namely, a large overlap leads to a large reduced resonance integral, and vice versa. Furthermore, large differences between the gas-phase and solution overlap lead in turn to large differences in the respective values of the reduced resonance integral. In other words, a good approximation of the overlap usually leads to good approximation of the reduced resonance integral, and vice versa. By looking separately at the results for each quantity, it is seen that the approximations are usually very good. More specifically, the overlaps and reduced resonance integrals in solution are almost identical to their corresponding values in
Figure 2. Absolute energy difference between eigenvalues calculated with and without the approximations for various geometries along the SN2 reaction.
gas phase, both in the SN2 reaction and in the CH3F dissociation case, reaching a maximum difference of 10% and typically being 3% and less. The two approximations are good, regardless of the geometry or the energy difference between the respective diabatic VB configurations, which can be very small (e.g., at the crossing point between two configurations) or very large. Furthermore, in the SN2 reaction, the approximations are shown to be good for any combination of VB configurations, suggesting that they are not limited to systems that are described with only two configurations but are valid even when an interplay between several VB configurations is involved. The LiF example displays, however, a different behavior, which is less ideal (see Figures 1S-g and 1S-h in the Supporting Information). Here, for various geometries along the dissociation curve, the solution values of both the overlap and the reduced resonance integrals are somewhat different from their corresponding values in gas phase, and the difference does not seem negligible. We believe, however, that this case is an exception. The nature of LiF changes considerably when moving from the gas phase to solution. While going through a homolytic dissociation into the respective radicals in gas phase, LiF dissociates into ions in solution. Hence, the LiF wave function, which is governed by the covalent contribution at large distances in the gas phase, should be dominantly ionic in solution. These major changes in the description of the wave function cause the approximations to be less valid in this case. Thus, using gas-phase values of overlap and reduced resonance integral instead of the actual values in the new environment (e.g., solution) seems like an excellent approximation for most cases where the environment affects the system but still keeps its general electronic nature. Yet, for systems the wave
VB/MMsThe Validity of the Underlying Approximations
J. Phys. Chem. B, Vol. 112, No. 39, 2008 12495
Figure 3. Weights of the various VB configurations that compose the overall wave function of the identity SN2 reaction. (a) gas-phase weights and (b) weights in solution calculated with the VB/MM (blue) and the DE-VB/MM (green) methods.
function of which undergoes extreme changes due to the surrounding, the approximations are expected to be less good. Examination of the Eigenvalues. The approximations regarding the overlap and the reduced resonance integral are designated to allow the calculation of the energetics, which is usually the quantity of interest. Therefore, comparison of the eigenvalues obtained with and without the approximations was carried out. This was done, as mentioned in the Theory and Methodology section, when considering only the electrostatic interactions with the environment. Figure 2 displays the energy difference ∆ε between eigenvalues that were calculated with and without the approximations (∆ε ) εacc - εapx). The values are given in kcal/mol for various geometries along the reaction coordinate of the SN2 reaction. ∆ε values for the two dissociation reactions are illustrated in Figure 2S in the Supporting Information. By looking at the results, one can see that for both the SN2 reaction and the dissociation of CH3F, the energy differences are very small (
