6304
J. Phys. Chem. 1996, 100, 6304-6309
Effect of Hydration on the Barrier to Internal Rotation in Formamide. Quantum Mechanical Calculations Including Explicit Solvent and Continuum Models
J. Simon Craw, Jonathan M. Guest, Matthew D. Cooper, Neil A. Burton, and Ian H. Hillier* Department of Chemistry, UniVersity of Manchester, Manchester, M13 9PL, U.K. ReceiVed: September 28, 1995; In Final Form: December 7, 1995X
The effect of hydration upon the barrier to rotation in formamide has been predicted using both continuum models and those that consider the solvent molecules explicitly. The continuum calculations are based upon the polarizable continuum model, while the simulation calculations employ a Monte Carlo (MC) treatment. Variants of the MC method are studied employing classical and quantum mechanical (both semiempirical and ab initio) descriptions of the solute. All methods suggest a different transition state structure in water compared to the gas phase and show the importance of solute polarization. The calculated free energy barrier to rotation is in good agreement with the experimental data in all cases.
∆Gsoln ) G(B) - G(A)
I. Introduction The rate of a chemical process is often strongly dependent on the nature and polarity of the solvent, particularly if there are intermediates along the reaction path having bonding and charge distribution significantly different from the reactants. However, most ab initio quantum chemical methods refer to isolated molecules, and although quantities such as barriers to rotation in the gas phase can now be routinely calculated to a high level of accuracy, the effects of solvation on such barriers have proved a more difficult problem. The gross effects of solvation, i.e. the polarization of the solute by the field of the solvent, can now also be included in quantum mechanical (QM) calculations, using the self-consistent reaction field (SCRF) model of Rivail et al.1-3 and the polarizable continuum model (PCM) of Tomasi et al.4 Continuum models place the solute inside a cavity and represent the solvent as a dielectric continuum which interacts with the charge distribution of the solute, producing a reaction field. Solvation energies thus calculated are dependent on the size and shape of the cavity and do not include potentially important explicit solute-solvent interactions (hydrogen bonding). Classical molecular dynamics (MD) and Monte Carlo (MC) liquid simulations seek to address these problems. By introducing molecular mechanical (MM) potentials to model solute-solvent and solvent-solvent interactions, explicit water molecules can be used. Most commonly the pairwise additive Coulomb plus Lennard-Jones (LJ) potential is adopted, requiring a charge, LJ radius, and LJ well depth parameters for each center, and results in a solute unable to polarize in the field of the solvent, although Kollman et al.5 have recently extended this model to include polarizable molecules. In this paper we examine various methods of estimating the effect of the solvent on the change in free energy for the process ∆Gsoln
Aaq 98 Baq
(1)
where A and B correspond to different structures of the same molecule, i.e., X
Abstract published in AdVance ACS Abstracts, March 15, 1996.
0022-3654/96/20100-6304$12.00/0
(2)
) Ggas(B) + ∆Gsolv(B) - Ggas(A) - ∆Gsolv(A)
(3)
) ∆Ggas + ∆∆Gsolv
(4)
Classical MC free energy perturbation (FEP) calculations are traditionally used to obtain ∆∆Gsolv. The results of such simulations are dependent on the atomic charges used in the empirical potential function. The use of charges derived from gas phase electronic structure calculations in liquid simulations may underestimate the solvation energy because such charges do not model the polarization of the solute by the solvent. A means of including the average effects of solvent polarization on the solute is to use ab initio QM calculations in which the solvent field is represented with a dielectric continuum to obtain solvation energies directly or, alternatively, to use the effective charges from QM continuum calculations in classical simulations. Scha¨fer et al.6 have used SCRF “polarized” charges in MC simulations to obtain the solvent contribution to the free energy of tautomerization of 3-hydroxypyrazole and found the use of polarized charges yields improved agreement with experiment. An alternative, and perhaps more satisfactory approach, is to include solute polarization by the solvent at a molecular level using a hybrid QM-MM method. Here the solute is described by either a semiempirical or ab initio QM method, and the solvent is modeled using MM. We have incorporated the hybrid QM-MM potential, where the semiempirical (AM1) method is used to model the solute in combination with the TIP3P water model for the solvent, as parametrized by Field, Bash, and Karplus,7 into a MC scheme which allows us to perform FEP simulations and obtain ∆∆Gsolv directly. We have also implemented an ab initio QM-MM approach, where a Hartree-Fock (HF) wave function is used to model the solute and MM used to model the solvent. To render the calculation computationally feasible, a full ab initio calculation is performed on a representative subset of the MC ensemble, and a thermodynamic cycle is used to obtain ∆∆Gsolv. The cycle we use is © 1996 American Chemical Society
Hydration and the Barrier to Rotation in Formamide
J. Phys. Chem., Vol. 100, No. 15, 1996 6305
Figure 1. Ground and transition states of formamide.
TABLE 1: Experimental Activation Energy (Ea, kcal/mol), Frequency Factor (log10 A), Free Energy of Solvation (∆Gsoln at 298 K, kcal/mol, Wia Eq 9), and Dielectric Constant (E) solventa 11
dioxane acetone11 MPK12 H2O11 neat11 diglyme12
Ea
log10 A
∆Gsoln
16.8 ( 1.0 16.9 ( 1.9 19.2 ( 0.2 21.3 ( 1.3 18.9 ( 1.0 19.7 ( 0.2
12.9 ( 0.6 12.1 ( 1.3 13.9 14.9 ( 1.3 12.9 ( 0.6 14.2
18.0 ( 1.4 17.8 ( 2.6 17.4 ( 0.2 18.4 ( 2.2 18.8 ( 1.3 17.8 ( 0.2
2.2 20.7 22.0 78.4 109.0
a MPK ≡ methyl propyl ketone, diglyme ≡ diethyleneglycol dimethyl ether.
AQM
∆GQM
98
∆GQMfMMV AMM
∆GMM
98
A number of high-level ab initio calculations of the gas phase barrier have been reported. These have identified the charge redistribution occurring on bond rotation and predict barriers for rotations Via TS1 of 14-15 kcal/mol. Thus, even solvents of very low dielectric constant (e.g. dioxane) produce a significant solvent contribution to ∆Gsoln. The effect of solvent on the barrier to rotation in both formamide and other amides such as N-methylacetamide (NMA) has been studied using both classical20,21 and hybrid semiempirical QM-MM MC simulations,22 as well as ab initio supermolecule13,16 and continuum methods.17-19,23,24 The greater solvation energy of the planar form compared to the nonplanar transition states and its origin are well-known. However, the different approaches often lead to different quantitative conclusions and frequently underestimate the magnitude of the barrier height due to hydration. For these reasons, we present here a comparison of a series of classical and quantum mechanical techniques for calculating solvent effects on one molecular property, the barrier to internal rotation in formamide.
BQM
II. Computational Considerations
v∆GMMfQM
While it is generally straightforward to locate stationary points on the potential energy (PE) surface for small gas phase structures, for condensed phase systems the problem is naturally more complex. Within continuum models, which involve fairly simple extensions to traditional ab initio codes, such stationary points can easily be obtained. However, the situation is more complex when solvent molecules are included explicitly, and in these models, the effect of solvation on a gas phase PE surface is generally considered. We here describe both continuum and calculations including explicit solvent molecules. A. Gas Phase and Continuum ab Initio Calculations. The gas phase barriers to rotation Via TS1 and TS2 were calculated at the MP2/6-31G** stationary points using a highly correlated wave function and large basis set (CCSD(T)/cc-pVTZ) with zero-point and thermal corrections evaluated using MP2/631G** harmonic frequencies at 298 K and 1 atm. Free energy barriers to rotation of 15.5 kcal/mol (TS1) and 16.9 kcal/mol (TS2) were obtained. In view of the small energy separation between these structures, it is not clear which is favored in water. In an attempt to resolve this question we have used the PCM method as implemented in Gaussian 9425 to obtain the stationary structures in aqueous solution. The calculations were carried out at the HF/6-31G** level using a cavity determined selfconsistently from the 0.001e isodensity surface and a dielectric constant of 78.4. B. Classical Monte Carlo Calculations. The classical MC simulations were performed with the BOSS program,26 i.e. isothermal-isobaric (NPT) ensembles with one solute molecule (with the internal degrees of freedom constrained at the MP2/ 6-31G** geometry) solvated by 212 TIP3P27 water molecules. Simulations were carried out at 25 °C and 1 atm pressure. Equilibration runs were performed for typically 2 × 106 configurations, followed by data collection for 2 × 106 configurations. The FEP simulations were performed both in the forward and backward directions using increments of 0.1 for the coupling parameter. The results quoted in the tables represent the average of the forward and backward runs. The potential function in our classical MC calculations comprises Coulombic and LJ terms and uses the OPLS28-31 combination rules,
BMM
(5)
i.e
∆GQM(AfB) ) ∆GQMfMM(A) + ∆GMM(AfB) + ∆GMMfQM(B) (6) Such a cycle has previously been used by Wesolowski and Warshel8 and by Gao9 to calculate free energies of solvation. Our procedure involves several hundreds of thousands of ab initio energy evaluations and is ideally suited to computers with massively parallel architecture. In this paper we choose to study the barrier to internal rotation of formamide, a prototype for the peptide (HNCO) bond and for which there have been several experimental10-12 and theoretical13-19 studies. The rotational barrier about the central C-N bond in formamide is solvent sensitive due to partial double-bond nature. The ground state of formamide is planar, and there are two possible transition states (TS1 and TS2) corresponding to a pyramidal nitrogen center (Figure 1). In the gas phase it is clear that the energetically favored transition state is the one involving the stabilizing intramolecular O‚‚‚H interaction (TS1). However it is not clear whether TS1 or TS2 is preferred in polar solvents due to the competition between intra- and intermolecular interactions. In Table 1 we have collected together experimentally determined values of the activation energy (Ea) and frequency factor (log10 A) in various solvents derived from 1H NMR experimental measurements.10,12 From these data we have obtained ∆Gsoln at 298 K Via
∆H ) Ea - RT
(7)
∆S ) R ln[Ah/kT] - R
(8)
and the usual expression relating ∆H, ∆S, and ∆G, i.e.
∆Gsoln ) Ea - RT ln[Ah/kT]
(9)
values of which are also listed in Table 1. There is clearly very little variation of ∆Gsoln with respect to the dielectric constant () of the solvent, all values being between 17 and 19 kcal/ mol.
Aij ) xAiAj, Bij ) xBiBj
(10)
6306 J. Phys. Chem., Vol. 100, No. 15, 1996
Craw et al.
TABLE 2: Potential-Derived Charges for Planar Formamide (au) center
qRHF g
H C N Ha Hb O
-0.01 0.75 -0.96 0.40 0.43 -0.61
TABLE 6: CHARMM Parameters Used To Model Formamide (E in kcal/mol, Rmin in Å)
RHF qaq
qMP2 g
MP2 qaq
center
type
Rmin
0.01 0.83 -0.96 0.42 0.44 -0.74
-0.02 0.62 -0.88 0.38 0.41 -0.51
-0.01 0.69 -0.87 0.41 0.42 -0.64
Ha
HA C N* NT
0.042 0.141 0.090 0.150 0.000 0.159
1.33 1.87 1.83 1.65 0.00 1.55
TABLE 3: Potential-Derived Charges for TS1 Formamide (au) center
qRHF g
RHF qaq
qMP2 g
MP2 qaq
H C N Ha Hb O
-0.05 0.90 -1.02 0.39 0.39 -0.61
-0.04 0.99 -1.07 0.41 0.41 -0.70
-0.06 0.82 -0.98 0.37 0.37 -0.52
-0.05 0.90 -1.03 0.39 0.39 -0.60
TABLE 4: Potential-Derived Charges for TS2 Formamide (au) center
qRHF g
RHF qaq
qMP2 g
MP2 qaq
H C N Ha Hb O
-0.10 0.89 -0.94 0.35 0.35 -0.55
-0.08 0.99 -1.03 0.39 0.39 -0.66
-0.13 0.82 -0.93 0.35 0.35 -0.46
-0.10 0.90 -1.01 0.39 0.39 -0.57
C Nb Nc H O a
H on C.
σ
CH N H O
0.105 0.170 0.000 0.210
3.75 3.25 0.00 2.96
(11)
and i and σi are the standard LJ energy and radius terms. Simulations were performed using the MP2/6-31G** optimized geometries and charges derived from both gas phase and continuum ab initio calculations Via fitting to the electrostatic potential using the scheme of Besler et al.32 and the CHELPG atomic radii.33 The charges from the gas phase RHF or MP2 and qMP2 density are denoted qRHF g g . Corresponding values from ab initio continuum calculations were obtained using the and PCM method with fixed atomic radii,34,35 giving qRHF aq qMP2 aq , the latter being obtained from an MP2 calculation based upon an RHF PCM reference wave function. Using these charges, it is possible to estimate the polarization free energy of the solute, ∆Gpol, Via FEP simulations in which the charges are mutated from their gas phase values to those obtained with the PCM wave function. The charges are listed in Tables 2, 3, and 4, while the LJ parameters used in these simulations are given in Table 5. The free energies obtained from the or qMP2 contain two distinct composimulations with qRHF aq aq nents; the free energy from the FEP simulations and a correction due to the change in internal energy of the solute required to produce the polarized charges. C. Semiempirical QM-MM Monte Carlo Calculations. In the hybrid QM-MM method the effective Hamiltonian for the solute-solvent system is defined as
ˆ°+H ˆ QM-MM H ˆ eff ) H
hybridized (GS). c sp3 hybridized (TS1, TS2).
{
}
qM zRqM ARM BRM H ˆ QM-MM ) -∑ + ∑ +∑ 12 iM riM RM RRM RM R R6RM RM
(13)
The first term on the right-hand side of eq 13 accounts for the interaction between solute electrons and the solvent center; riM is the distance between electron i and center M with partial charge qM. The second and third terms in eq 13 are the classical interaction of the solute nuclei (charge zR) with the solvent partial charge (RRM is the distance between the QM center R and the MM center M) and the solute-solvent LJ term. For this last term we use the CHARMM36 atomic parameters and combination rules, i.e. 12 6 , Bij ) 2ijRmin,ij Aij ) ijRmin,ij
(14)
ij ) x(ij), Rmin,ij ) Rmin,i + Rmin,j
(15)
with
where
Ai ) 4iσi12, Bi ) 4iσi6
sp2
where H ˆ ° is the QM Hamiltonian of the free solute and H ˆ QM-MM is the QM-MM solute-solvent interaction.
TABLE 5: OPLS Parameters Used To Model Formamide (E in kcal/mol, σ in Å) center
O b
(12)
where i is the LJ well depth parameter for atom i, Rmin,ij is the value of Rij where the LJ energy is a minimum (ij), and the LJ potential is zero at 2-1/6Rmin,ij. The CHARMM force field is used in preference to OPLS parameters and combination rules, as the five adjustable QM parameters in the semiempirical QMMM potential were optimized by Field et al.7 for the CHARMM parameter set and combination rules. The LJ parameters used in the QM-MM FEP simulations are listed in Table 6. QM ) in a QM-MM MC The solute-solvent energy (Esol-solv simulation is the expectation value of H ˆ eff minus the internal energy of the solute, QM ) 〈Φ|H ˆ eff|Φ〉 - 〈Φ°|H ˆ °|Φ°〉 Esol-solv
(16)
where Φ is the solvated solute wave function and Φ° the gas phase solute wave function. As a consequence, the solute internal energy effectively cancels and hence the quality of the QM wave function may not be paramount. Of more importance is the ability of the solute to polarize in the presence of the solvent and the value of the LJ parameters which model the solute-solvent QM exchange repulsion. The solute polarization free energy (∆Gpol) is obtained from the QM-MM ensemble as
∆Gpol ) RT ln〈e-[E°-E
QM]/RT
〉QM
(17)
where
E° ) 〈Φ°|H ˆ eff|Φ°〉 and
(18)
Hydration and the Barrier to Rotation in Formamide
EQM ) 〈Φ|H ˆ eff|Φ〉
J. Phys. Chem., Vol. 100, No. 15, 1996 6307
(19)
QM The calculation of Esol-solv from eq 16 implies a full QM calculation for each MC configuration. However, the computational effort required for the large number of energy evaluations can be reduced considerably. The evaluation of the solute one- and two-electron integrals can be avoided when the new MC configuration is obtained Via a solvent move with the solute geometry remaining unchanged. Furthermore, a full selfconsistent-field (SCF) procedure can be reduced to a single iteration if the solvent move occurs outside a certain cutoff radius from the solute. In this case the effect of that move on the charge distribution of the solute will be negligible, and QM can be evaluated with the existing density. A cutoff Esol-solv of 4-5 Å for this quick energy evaluation introduces a negligible error, but obviates the necessity of iterating in about 75% of the SCF calculations.37 The rate-determining step in a QMMM calculation with a large MM environment is now the calculation of H ˆ QM-MM at each step. Practical implementations of the MC algorithm only calculate the change in energy resulting from a move, in order to avoid recalculating the total energy of the system at each step. This concept has been extended to the calculation of the H ˆ QM-MM matrix. Since H ˆ QM-MM contains a contribution from each solvent molecule and consecutive configurations differ by the displacement of one solvent molecule, the new H ˆ QM-MM matrix is obtained from the current H ˆ QM-MM matrix and the change due to the solvent move. For example if a new configuration is generated by moving solvent molecule A from R to R′, then the new H ˆ QM-MM can be obtained Via
A A ˆ QM-MM(old) - H ˆ QM-MM (R) + H ˆ QM-MM (R′) H ˆ QM-MM ) H (20) A )H ˆ QM-MM(old) + δH ˆ QM-MM
(21)
TABLE 7: PCM-Optimized Geometries (Bond Lengths in Å, Angles in Deg, Gas Phase RHF Values in Parentheses) RCH RCN RCO RNH1 RNH2 θHCN θCNH1 θCNH2 θNCO δHCNH1 δHCNH2 δH2NCO
GS
TS1
TS2
1.09 (1.09) 1.33 (1.34) 1.21 (1.19) 0.99 (0.99) 0.99 (0.99) 113.6 (112.8) 121.3 (121.6) 119.6 (119.1) 124.5 (124.9) -0.1 (0.0) -179.9 (180.0) 0.1 (0.0)
1.09 (1.09) 1.42 (1.43) 1.19 (1.18) 1.01 (1.00) 1.01 (1.00) 114.2 (113.5) 108.3 (108.6) 108.3 (108.6) 124.8 (125.0) -122.9 (-122.7) 122.9 (122.7) -57.1 (-57.3)
1.09 (1.10) 1.41 (1.42) 1.19 (1.18) 1.01 (1.00) 1.01 (1.00) 117.1 (116.4) 108.3 (109.9) 108.3 (109.9) 122.5 (123.2) -56.8 (-58.6) 56.8 (58.6) -123.3 (-121.4)
free energies of solvation Via the cycle eq 5. To render the simulations tractable, energy evaluations were performed on a subset of the classical Monte Carlo ensemble. To generate this ensemble, classical MC simulations (using OPLS parameters and qRHF charges) were carried out on each structure. Followg ing equilibration, every 10th accepted configuration from a run of 1 × 106 configurations was taken, and its energy evaluated using an ab initio hybrid QM-MM approach employing a HF/ 6-31G** calculation for the solute and TIP3P charges for the water. The solute geometry remained fixed for all of those sample configurations. We have compared the average of the classical solute-solvent energy for the subset with that of the full MM ensemble as a check that the subset is representative of the whole and found insignificant differences. The independence of each configuration enables the QM energies to be computed in parallel and with great efficiency. Since the solute geometry is fixed, all two-electron integrals were retained in memory and each SCF procedure started with the gas phase density. The MM to QM free energy change is calculated from eq 22:
∆GMMfQM ) -RT ln〈e-[E
QM
sol-solv
- EMMsol-solv]/RT
〉MM
(22)
A H ˆ QM-MM
refers to the contribution solvent molecule A where ˆ QM-MM need makes to H ˆ QM-MM. Thus, only the change in H be calculated, each move resulting in a significant computational saving. With this procedure a small numerical error builds up in H ˆ QM-MM, and to minimize this, every time a solute or volume move is made (typically every 50-100 for solute or 10001500 for volume), the whole of the H ˆ QM-MM matrix is recomputed. A run of 2 × 106 moves for a QM-MM FEP MC calculation involves 1 × 105 integral evaluations, 1 × 106 full energy evaluations, and 5 × 106 quick energy evaluations and requires about 2 h using one processor of a Silicon Graphics Power Challenge computer. The QM-MM simulations were run for 1 × 106 configurations of equilibration, followed by 2 × 106 configurations for data collection. FEP calculations were performed in both the forward and backward directions using increments of 0.1 for the coupling parameter. The results quoted represent the average of the forward and backward runs. D. Ab Initio Hartree-Fock QM-MM Monte Carlo. A QM-MM MC simulation in which the solute is described using an ab initio Hartree-Fock or density functional theory method, rather than a semiempirical scheme, is naturally to be preferred. However, such a development presents formidable problems in terms of computational resources. We use here a scheme which employs a number of approximations designed to render it computationally feasible. We believe that the approximations are reasonable, although they can be rigorously justified only by comparison with a full ab initio QM MC simulation. Ab initio HF MC simulations have been performed to obtain the
QM is the QM solute-solvent energy (eq 16) and where Esol-solv MM Esol-solv is the corresponding MM energy. The free energy of polarization is calculated from
∆Gpol ) ∆GMMfQM + RT ln〈e-[E°,
QM
sol-solv
- EMMsol-solv]/RT
〉MM (23)
where °,QM ) 〈Φ°|H ˆ eff|Φ°〉 - 〈Φ°|H ˆ °|Φ°〉 Esol-solv
(24)
The QM calculations were carried out using a modified version of the Gaussian 9238 suite of programs on the Cray T3D massively parallel computer at the Edinburgh Parallel Computing Centre. III. Computational Results The objective of this study is to identify the effects responsible for the increase in barrier height upon hydration and to accurately predict their magnitude. We first discuss the results of the PCM calculations in terms of the minimum energy and transition state structures and charge distribution. In Table 7 we list the geometric parameters for all three PCM optimized stationary points, with the corresponding gas phase RHF values in parentheses. In both the gas phase and PCM calculations the C-N bond is considerably lengthened in TS1 and TS2 compared with the ground state; concurrently the C-O bond
6308 J. Phys. Chem., Vol. 100, No. 15, 1996
Craw et al.
TABLE 8: Ab Initio and Monte Carlo Calculations of ∆Gpol (kcal/mol)a method b
PCM AM1 QM-MM ab initio QM-MM FEP (RHF)c FEP (MP2)d
GS
TS1
TS2
-5.0 -2.2 ( 0.2 -5.1 ( 0.1 -5.0 ( 0.1 -3.2 ( 0.2
-2.6 -1.2 ( 0.1 -2.0 ( 0.1 -2.1 ( 0.1 -0.3 ( 0.1
-4.8 -1.9 ( 0.1 -2.9 ( 0.2 -3.9 ( 0.2 -2.0 ( 0.2
TABLE 10: Monte Carlo FEP Estimates of ∆Gsoln (kcal/ mol) Using MC FEP Values of ∆∆Gsolv and the Best Gas Phase Value of ∆Ggasa,b ∆Gsoln potential MM MM MM MM QM-MM QM-MM expt
a Error bars estimated from difference in forward and backward FEP. PCM calculation with cavity defined using fixed atomic radii. c Classical FEP simulation, perturbing the charges from RHF to PCM. d Classical FEP simulation, perturbing the charges from MP2 to MP2(PCM). b
TABLE 9: Ab Initio Estimates of ∆Gsoln (kcal/mol) Using PCM Values of ∆∆Gsolv and the Best Gas Phase Value of ∆Ggasa ∆Gsoln method
TS1
TS2
PCMb MP2(PCM)c PCMd expt
17.7 17.3 19.0
16.7 16.8 16.8
18.4 ( 2.2
a
Using values for ∆Ggas of 15.5 kcal/mol (TS1) and 16.9 kcal/mol (TS2). b PCM calculation with cavity defined using fixed atomic radii. c MP2 based upon a calculation with cavity defined using fixed atomic radii. d PCM calculation with cavity defined self-consistently using the 0.001e isodensity surface.
length in TS1 and TS2 is shorter than in GS. The PCM reduces the C-N bond length in all three structures by 0.01 Å and lengthens the C-O bond by between 0.01 and 0.02 Å. The PCM also predicts slightly longer N-H bond lengths in TS1 and TS2 compared with the gas phase RHF values. Turning now to the potential-derived charges (Tables 2, 3, and 4), there is clearly an increase in the polarity of each species upon solvation. The main effect is an increase in the oxygen and a reduction in the carbon electron density. This effect is somewhat greater in the planar structure compared with TS1 and TS2, and also greater in TS2 compared with TS1, since the latter is already polarized due to intramolecular O‚‚‚H interactions. These differing charge distributions are reflected in the differing polarization energies calculated from both the continuum and hybrid QM-MM models (Table 8). The values are significantly greater for the ground state and TS2 than for TS1, with the ground state values being somewhat greater than for TS2. The close similarity of these energies for all the treatments is encouraging, with the values for the MP2 treatment being somewhat smaller than the corresponding RHF values, as expected. The ab initio QM-MM calculation gives somewhat larger polarization energies than the semiempirical AM1 calculation, reflecting the more flexible atomic basis of the former. The values of ∆Gsoln obtained using eq 4, i.e., combining the best gas phase estimate of ∆Ggas with the various continuum calculations estimate of ∆∆Gsolv, are given in Table 9. Three different PCM calculations were used to obtain ∆Gsoln, firstly using a cavity shape obtained from overlapping atomic radii, secondly performing an MP2 calculation using this PCM charge density, and lastly using a cavity shape defined self-consistently as the 0.001e isodensity surface. All three continuum models predict the rotation to occur Via TS2 in aqueous solution, in contrast to the gas phase, where the lowest energy path is Via TS1. This reflects the competition between intramolecular and intermolecular hydrogen bonding, although these are not considered explicitly in the continuum models. In the gas phase
solute charges
TS1
TS2
18.0 ( 0.3 18.1 ( 0.1 19.2 ( 0.1 18.0 ( 0.2 17.3 ( 0.2 16.8 ( 0.3 20.9 ( 0.7 18.1 ( 0.3 18.2 ( 0.2 16.6 ( 0.2 20.4 ( 0.3 18.9 ( 0.1 18.4 ( 2.2
qRHF g RHF qaq MP2 qg MP2 qaq AM1 ab initio
a Error bars estimated from difference in forward and backward FEP. Using values for ∆Ggas of 15.5 kcal/mol (TS1) and 16.9 kcal/mol (TS2). Error bars represent the standard deviation of the forward and backward runs. b
TABLE 11: Dipole Moments in Solution (D) method
GS
TS1
TS2
PCM MP2(PCM) AM1 QM-MM ab initio QM-MM
5.2 4.8 4.9 5.2
2.1 1.6 1.5 1.9
5.2 4.6 4.6 5.1
the barrier Via TS1 is lower due to the favorable O‚‚‚H intramolecular interactions, while in aqueous solution TS2 is preferentially solvated compared with TS1. The additional solvation energy of TS2 (with respect to TS1) is enough to overcome the gas phase difference (which we estimate at only 1.4 kcal/mol), making the path Via TS2 preferred in aqueous solution. In Table 10 the results of the Monte Carlo free energy perturbation calculations (with various potentials) of ∆∆Gsolv have been combined with our best gas phase estimate of ∆Ggas to obtain ∆Gsoln Via eq 4. All of the classical MC FEP simulations are in excellent agreement with the experimental value of ∆Gsoln. It is noted here that due to the lack of solute polarization by the solvent, the use of gas phase charges underestimates the solvation energy of TS2 compared to TS1. The use of gas phase RHF charges (qRHF g ) predicts the rotation to occur Via TS1, although the two paths are practically isoenergetic. The other three charge sets all predict the path Via TS1 to be about 1 kcal/mol higher than that Via TS2. Turning to the QM-MM results, both the AM1-TIP3P QMMM results, obtained Via direct QM-MM FEP calculations, and the ab initio QM-MM results, calculated Via the thermodynamic cycle eq 5, are in accord in predicting the lowest energy path to be Via TS2 and are close to those from the PCM calculations (Table 9) and the Monte Carlo MM simulations. The similarity between the results of the ab initio PCM calculations and the results of the simulations with the QMMM potentials is also reflected in the aqueous phase dipole moments (Table 11). In Figure 2 we plot various radial distribution functions (rdf), gXY(R), where X ) C, N, or O of formamide and Y ) O of the water, from the AM1-TIP3P QM-MM simulations. In the ground state (Figure 2A), as expected, there is no evidence of direct hydrogen bonding to the nitrogen or carbon and approximately two hydrogen bonds to the carbonyl oxygen, slightly more than that reported by Jorgensen and Swenson.28 For TS1, Figure 2B, the hydrogen bonding to the carbonyl oxygen is considerably reduced compared with planar formamide, a reflection of increased intramolecular interaction. There is also evidence of hydrogen bonding to the (now) tetrahedral nitrogen, indicated by the small peak in the gNO(R) function at about 3 Å. For TS2, Figure 2C, there are again two hydrogen
Hydration and the Barrier to Rotation in Formamide
J. Phys. Chem., Vol. 100, No. 15, 1996 6309 Acknowledgment. We thank EPSRC for support of this research, Prof. W. L. Jorgensen for the use of the BOSS program, and Dr. A. J. Masters for helpful discussions. References and Notes
Figure 2. Solute-solvent radial distribution functions gXY(R) for formamide GS (A), TS1 (B), and TS2 (C).
bonds to the carbonyl oxygen and one hydrogen bond to the nitrogen center, and the rdf is overall very similar to that for the ground state, consistent with the ∆∆Gsolv and dipole data which indicate that the GS and TS2 are equally well hydrated.
IV. Discussion and Conclusion The calculations described here allow a comprehensive comparison between the various methods now currently employed to treat the chemically important problem of solvation by water. To our knowledge, such a comparison has not previously been reported. It is quite remarkable how quantitatively close the predictions from different methods are. Those that employ a quantum mechanical treatment to describe the solute, using either explicit water molecules or a continuum model, and those based upon a classical MC treatment all predict values within the error limits of the experiment. It is useful to note that our ab initio Hartree-Fock QM-MM Monte Carlo calculations were rendered computationally feasible by the use of the in-core SCF procedure. Thus, using the Cray T3D having 64 Mbytes of memory on each node, the method is restricted to a “problem” the size of formamide studied here. For substantially larger systems more rapid methods, such as semiempirical MO or empirical valence bond schemes, are viable alternatives. Turning to the prediction concerning the mechanism of the rotation in formamide, the calculations suggest that the transition state structure is likely to be different in the gas and aqueous phases. However, the two transition states are quite close in energy in water, and the level of accuracy of the present calculations cannot assign, with complete confidence, the most stable one. The important contributions of differential polarization of the solute by the solvent to the barrier in solution show the need to include these effects Via quantum mechanical models. Finally the consistent picture that has emerged from the calculations shows that the use of these methods to predict pathways of actual chemical reactions in solution now has real potential.
(1) Rinaldi, D.; Rivail, J. L.; Rguini, N. J. Comput. Chem. 1992, 13, 675. (2) Rinaldi, D. Comput. Chem. 1982, 6, 155. (3) Rivail, J. L.; Terryn, B. J. Chim. Phys., Phys. Chim. Biol. 1982, 79, 1. (4) Miertus, S.; Scrocco, E.; Tomasi, J. Chem. Phys. 1981, 55, 117. (5) Sun, Y.; Caldwell, J. K.; Kollman, P. A. J. Phys. Chem. 1995, 99, 10081. (6) Cao, M.; Teppen, B. J.; Miller, D. M.; Pranata, J.; Scha¨fer, L. J. Phys. Chem. 1994, 98, 11353. (7) Field, M. J.; Bash, P. A.; Karplus, M. J. Comput. Chem. 1990, 11, 700. (8) Wesolowski, T. A.; Warshel, A. J. Phys. Chem. 1994, 98, 5183. (9) Gao, J. J. Phys. Chem. 1992, 96, 537. (10) Sunners, B.; Piette, L. H.; Schnider, W. G. Can. J. Chem. 1960, 38, 681. (11) Kamei, H. Bull. Chem. Soc. Jpn. 1968, 41, 2269. (12) Drakenberg, T.; Forsen, S. J. Phys. Chem. 1970, 74, 1. (13) Hoffman, H. J.; Peinel, G. J. Mol. Struct. (THEOCHEM) 1983, 8, 289. (14) Tsuzuki, S.; Tanabe, K. J. Chem. Soc., Perkin Trans. 2 1991, 1255. (15) Wiberg, K. B.; Breneman, C. M. J. Am. Chem. Soc. 1992, 114, 831. (16) Wang, X. C.; Facelli, J. C.; Simons, J. Int. J. Quantum Chem. 1993, 45, 123. (17) Duben, A. J.; Miertus, S. Chem. Phys. Lett. 1982, 88, 395. (18) Burton, N. A.; Chiu, S. S.-L.; Davidson, M. M.; Green, D. V. S.; Hillier, I. H.; McDouall, J. J. W.; Vincent, M. A. J. Chem. Soc., Faraday Trans. 1993, 89, 2631. (19) Hall, R. J.; Davidson, M. M.; Burton, N. A.; Hillier, I. H. J. Phys. Chem. 1995, 99, 921. (20) Duffy, E. M.; Severance, D. L.; Jorgensen, W. L. J. Am. Chem. Soc. 1992, 114, 7535. (21) Morgantini, P.-Y.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 6057. (22) Gao, J. J. Am. Chem. Soc. 1993, 115, 2930. (23) Luque, F. J.; Orozco, M. J. Org. Chem. 1993, 58, 6397. (24) Orozco, M.; Bachs, M.; Luque, F. J. J. Comput. Chem. 1995, 16, 563. (25) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Cheeseman, J. R.; Robb, M. A.; Keith, T. A.; Petersson, G. A.; Montgomery, J. A.; Raghavachari, K.; Al-Laham, M. A.; Zakrzewski, V. G.; Ortiz, J. V.; Foresman, J. B.; Cioslowski, J.; Stefanov, B. B.; Nanayakkara, A.; Challacombe, M.; Peng, C. Y.; Ayala, P. Y.; Chen, W.; Wong, M. W.; Andre´s, J. L.; Replogle, E. S.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Binkley, J. S.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Head-Gordon, M.; Gonzalez, C.; Pople, J. A. Gaussian 94, Revision A.1; Gaussian, Inc.: Pittsburgh, PA, 1995. (26) Jorgensen, W. L., Boss Version 3.5; Yale University: New Haven, CT, 1994. (27) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (28) Jorgensen, W. L.; Swenson, C. J. J. Am. Chem. Soc. 1985, 107, 1489. (29) Jorgensen, W. L.; Swenson, C. J. J. Am. Chem. Soc. 1985, 107, 569. (30) Jorgensen, W. L.; Gao, J. J. Am. Chem. Soc. 1988, 110, 4212. (31) Jorgensen, W. L.; Tirado-Rives, J. J. Am. Chem. Soc. 1988, 110, 1657. (32) Besler, B. H.; Merz, K. M.; Kollman, P. A. J. Comput. Chem. 1990, 11, 431. (33) Breneman, C. M.; Wiberg, K. B. J. Comput. Chem. 1990, 11, 361. (34) Green, D. V. S. Ph.D. Thesis, University of Manchester, 1991. (35) Orozco, M.; Jorgensen, W. L.; Luque, F. J. J. Comput. Chem. 1993, 14, 1498. (36) Brooks, B. R.; Burrccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. J. Comput. Chem. 1983, 4, 187. (37) Gao, J. J. Am. Chem. Soc. 1995, 117, 8600. (38) Frisch, M. J.; Trucks, G. W.; Head-Gordon, M.; Gill, P. M. W.; Wong, M. W.; Foresman, J. B.; Johnson, B. G.; Schlegel, H. B.; Robb, M. A.; Replogle, E. S.; Gomperts, R.; Andre´s, J. L.; Raghavachari, K.; Binkley, J. S.; Gonzalez, C.; Martin, R. L.; Fox, D. J.; Defrees, D. J.; Baker, J.; Stewart, J. J. P.; Pople, J. A. Gaussian92; Gaussian, Inc.: Pittsburgh, PA, 1992.
JP952896D
