J. Phys. Chem. 1994, 98, 13772-13779
13772
Polarization Effects on the Tautomeric Equilibria of 2- and 4-Hydroxypyridine in Aqueous and Organic Solution Jiali Gao* and Lei Shao Department of Chemistry, State University of New York at Buffalo, Buffalo, New York 14260 Received: September 14, 1994@
The role of electronic polarization effects on the tautomeric equilibria of 2- and 4-hydroxypyridine has been examined through Monte Carlo simulations using a combined quantum mechanical and molecular mechanical (QM/MM) potential. The results indicate that differential polarization in the hydroxy and oxo forms contributes significantly to the overall observed solvent effects on the tautomeric equilibria, amounting to 14%-68% of the total free energy of solvation in chloroform and in water. This is accompanied by concomitant changes in molecular dipole moment. Detailed insights into the structural and energetic nature of the solvent polarization effects are presented.
1. Introduction The tautomeric equilibria of heterocyclic compounds are of considerable importance in chemistry and biochemistry. Furthermore, solvent effects often play a major role in determining the position of eq~ilibrium.~,~ Undoubtedly, one of the most extensively studied cases is the hydroxypyridinetpyridone
H 2
1
OH
0
Q-0 I
H
3
4
In the gas phase, there is a preference of the hydroxy form in both 2- and 4-hydroxypyridine, while in nonpolar solvents such as cyclohexane and chloroform the two tautomers exist in comparable am0unts.43~However, the two equilibria are shifted entirely to the right in favor of the oxo form, 2- and 4-pyridone, respectively, in polar solvents as well as in the crsytalllne state.4~~ Similar solvent effects are found in other heterocyclic compounds including nucleotide bases.2J0 Consequently, quantitative prediction of the solvent dependence of the tautomeric equilibria in 2- and 4-substituted pyridine represents a critical test of computational methods in condensed phase simulations. Previous studies have focused on continuum, reaction field models in ab initio and semiempirical molecular orbital calculat i o n ~ . ~ *Although * ~ ~ ~ these methods can yield quantitative predictions and often provide valuable qualitative insights into electronic properties, specific solute-solvent interactions are excluded. In many cases, inclusion of one or several solvent molecules is necessary in order to reproduce the experimental data.'* Recently, use of combined quantum mechanical and molecular mechanical (QM/MM) methods has become popular thanks to advances in computer t e c h n ~ l o g y . ~ ~In- ~this ~ approach, the solute molecule is treated quantum mechanically in a similar fashion as that in the continuum self-consistent @Abstractpublished in Advance ACS Abstracts, November 15, 1994.
0022-365419412098-13772$04.5010
reaction field (SCRF) model. However, the surrounding solvent molecules are also treated explicitly, albeit using a molecular mechanics-type force field. Now, the combined QM/MM simulation method has the advantage of taking into account both the solute electronic properties including polarization, which is of particular interest in tautomerism of heterocyclic compounds, and computational efficiency of the MM force fields in statistical Monte Carlo and molecular dynamics simulations. In this paper, we report the results from analyses of the solute electronic polarization in the tautomeric equilibria of 2- and 4-hydroxypyridinelpyridone in water and in chloroform. Although this system has been extensively studied, the present investigation yields new insights into the molecular electronic polarization effects on tautomeric equilibria of heterocyclic compounds in general. Section 2 summarizes the theoretical model used in the present calculation, while the performance of the combined AMl/MM potential is verified in section 3 via comparison with ab initio 6-3 1G(d) results for bimolecular complexes with water. Findings from the fluid simulations are presented next, followed by conclusions.
2. Methods Intermolecular Potential Functions. In the present study, a combined QM/MM potential is employed by partitioning the condensed phase system into a quantum mechanical region and a molecular mechanical region.13-16 In the QM region, the solute molecule is treated by the semiempirical Austin Model 1 (AM1) theory," while the solvent molecules in the MM area are represented by the three-site TIP3P model for water or by the OPLS potential for chloroform.18 Thus, the total effective Hamiltonian of the system is given by
+ +
Heff= ifX HX, Hs, where Ei", is the Hamiltonian for the isolated QM solute, Hx, is the solute-solvent interaction Hamiltonian, and H,,represents the classical solvent-solvent interaction energy. Throughout this paper, the subscript X represents the solute molecule in the QM region, whereas s designates the classical solvent molecules that are treated by Fmpirical potentials. The QM/ MM interaction Hamiltonian Hx,can be further divided into an electrostatic, which is included in solving the Hartree-Fock equations, and a van der Waals term, 0 1994 American Chemical Society
Tautomeric Equilibria of 2- and 4-Hydroxypyridine
Consequently, the solute-solvent interaction energy is given as follows
J. Phys. Chem., Vol. 98, No. 51, 1994 13773 SCHEME 1: Thermodynamic Cycle Used in Computation of Free Energies Solvation between Pyridone (B) and Hydroxypyridine (A) MG8.3
A(HiL6Y)
-
B(M%6:w,
AGqmimm(A)
where Ysolis the wave function of the solute in solution, which minimizes the energy of the Hamiltonian of eq 1. The LennardJones terms account for dispersion interactions between the QM and MM regions, which contain empirical parameters (0and E) in the present s t ~ d y . ’ ~ ,These ’ ~ values have been derived previously by comparing with ab initio results for bimolecular complexes and are adopted here.15a Free Energy Calculations. The difference in free energy of solvation between the two tautomers in each equilibrium is determined using statistical perturbation theory in the context of Monte Carlo simulation^.'^ Although it is possible to design a direct proton-transfer path in the calculation, we choose to use a procedure described previously for free energy simulations using combined QM/MM potentials.20 In this procedure (Scheme l), the electrostatic term of the QM/MM interaction Hamiltonian (eq 2), is first annihilated, leaving only the Lennard-Jones terms, which are then gradually “mutated” from one tautomer into another. The latter calculations consist of, of course, empirical potentials, and standard methods are employed for the conversion.21 To investigate the solvent polarization effects, we compute the electronic polarization contribution to the total free energies of solvation via eq 4.
?ds,
free energy was computed with alternate forward and reverse perturbations via a double-wide sampling method.21b Ten simulations for each molecule were used in the “electrostatic decoupling” (vertical steps in Scheme 1),21awhile ten additional steps were executed to interconvert the “van der Waals tautomers”. For each simulation, 5 x lo5 to lo6configurations were taken in the equilibration, which were followed by 1.5 x lo6 configurations for data collection. During the fluid simulation, solute-solvent interaction energies are evaluated by singlepoint Hartree-Fock SCF calculations with the effective Hamiltonian of eq 1, and solvent molecules that are 10 8, (in water) or 13 8, (in chloroform) within any solute atoms are included. All calculations were carried out on a SUN sparc 10 and IBM RS6000/370 computers in our laboratory using the BOSS/ MCQUB program?3 in which the QM energies are determined by the MOPAC program.24 The required ab initio calculations for bimolecular complexes as well as gas phase energies were performed on a Cray C90 at Pittsburgh Supercomputing Center using Gaussian 90.25
3. Gas Phase Bimolecular Complexes
An important criterion on the validity of the AMlhIM model for liquid simulations is to be able to predict the structural and energetic properties of bimolecular complexes?6 Consequently, intermolecular geometry optimizations were performed for complexes of tautomers in the hydroxypyridine system with water. For each species, two or three orientations of the water where Yoand Ysolare the wave functions of the solute molecule hydrogen-bonded complex are considered with 6-3 1G(d) partial in the gas phase and in solution, respectively, and E(Yo)and geometry optimization (Figure 1). Since the 6-31G(d) basis E(Y,,I) are defined by eqs 5 and 6. set is known to provide good descriptions of hydrogen-bonding interaction^,^^ the results are used to verify the performance of E(@> = (yo IfiefiI @) (5) the combined AMl/TIP3P model. In the latter calculations, the heterocyclic compounds are treated by the AM1 theory, whereas E W S O I ) = (Yollfieffl~sol> (6) water is represented by the TIP3P model. Monomer structures are kept fixed during the partial geometry optimization at the Clearly, E(Ys01)and E(Yo)are energies of the polarized solute and “nonpolarized” solute in the same solvent c~nfiguration.’~~ corresponding level of theory, while the experimental configuration is adopted for water. The results are shown in Figure 1. The brackets ( )Y,, in eq 4 indicate that the ensemble average is computed using the fully relaxed solute wave function. The accord for the interaction energies between the combined AMl/TIP3P and ab initio 6-31G(d) results is shown in Figure It should be pointed out that the polarization free energy 2. The overall agreement resembles previous findings for a evaluated using eq 4 may be an overestimate of the effects because the solvent reorganization term, which results from variety of hydrogen-bonded complexes.26 For the present solvent configuration changes following the solute electronic system in particular, the root-mean-square (rms) deviation is 0.8 kcdmol for a total of 10 structures covering an energy range polarization, is not well averaged. In short, the rate of convergence of eq 4 might be slow. However, AGpol can be of about -3 to -8 kcdmol. The interaction between the accurately determined by computing solvation free energies of organic species and water may be divided into three energetic both the polarized and nonpolarized solute (using the fixed gas groups. Interactions involving the hydroxyl and carbonyl groups phase wave function). For the present analyses, the use of eq form the strongest hydrogen bonds with binding energies of 4 should be sufficient. 7-8 kcdmol. The next strong interactions are between hydroxypyridine nitrogen and water (4.5-6.5 kcaumol). The Monte Carlo Simulations. Statistical mechanical Monte weakest hydrogen bonds are between the hydroxyl oxygen and Carlo simulations have been performed employing Metropolis water hydrogen, which amount to about 3 kcaumol. The AM1/ sampling in the isothermal-isobaric (NPT) ensemble at 25 “C TIP3P model tends to underestimate the interaction energy by and 1 atm. A cubic box consisting of 506 H20 molecules (ca. about 0.5-1 kcdmol in comparison with ab initio 6-31G(d) 25 x 25 x 25 A3), or 260 CHCl3 molecules (ca. 31.5 x 31.5 x 3 1.5 A3),plus one solute was used. To facilitate the statistics results,26 except in the case of hydroxyl-acceptor bonds. Interaction distances from the AMl/TIP3P calculations are in near the solute, Owicki-Scheraga preferential sampling techgood agreement with the ab initio data for the carbonyl group; nique with l/($ C), where C = 150 A2, was adopted.22 The
+
Gao and Shao
13774 J. Phys. Chem., Vol. 98, No. 51, I994
E E
-8.28 1.7.451
-2.92 (-3.77)
I
2.13 (1.75)
j
i
(-8.96)
H
2-Hydroxypyridine-Water Complexes
I
E
-1.70 (-1.601
I
/"
4-Hydroxypprldine-Water Complexes
H
I
1.96 (1.93)
/ O
E I -8.01 (.6.86)
9
7
E = -8.03
I H
(-7.13)
2.02 '(1.68) ~
E
-7.49 (-6.48)
E
-7.22 (-7.01)
I
H H
2-Pyridona-Water Complexes
4-Pyridone-Water Complexes
Figure 1. Bimolecular complexes between hydroxypyridine/pyridone and water. 6-31G(d)results are first given, followed by AMl/TIP3P values
in parentheses. Interaction energies are given in kcdmol, distances in A, and angles in degrees.
-2.0 -4.0
I c ++
%
B
5
-6.0 -
9 -8.0
to the total interaction energy. In all, the combined AMlRIP3P potential appears to perform as well as molecular mechanics force fields in describing intermolecular interactions in the present system. There have been early studies of hydrogen-bonding interactions between pyridone and Recently, Del Bene reported the results of a high-level ab initio study of hydrogenbonded complexes of 2-pyridone with water at MEWcc-pVTZf/ /MP2/6-31+G(d,p) level.2g In that study, geometries are fully optimized and two equilibrium structures were found. In the first complex, a water molecule acts as both a hydrogen bond donor to the carbonyl and acceptor from the N-H unit of 2-pyridone with a binding energy of 12.0 kcal/mol. The second structure consists of a single hydrogen bond from water to the carbonyl group in 2-pyridone. The interaction energy was predicted to be 8.4 kcal/mol. The enthalpies at 298 K for these two complexes are 10.2 and 6.7 kcal/mol. In our calculation, the primary concern is to test the quality of the combined QM/ MM potential. Thus, only partial geometry optimizations were performed with fixed monomer geometries.
-
- 10.0-10.0
-8.0
-6.0 6-3 1G(d)
-4.0
-2.0
Figure 2. Comparison between ab initio 6-31G(d) and AMIITIP~P
interaction energies (kcal/mol). however, they are about 0.4 8, shorter in other pairs. Note that similar features are necessary in empirical force fields such as the OPLS potential, with which the hydrogen bond distances are typically 0.2-0.3 8, shorter than the 6-31G(d) counterpart.28 Although it seems that the short hydrogen bond distances predicted by the AMl/TIP3P model could affect the computed polarization energies, resulting in an exaggeration in this effect, it should be noted that the binding energies from this model are in accord with ab initio results. It might be expected that the contribution of the polarization effects be proportional
4. Tautomeric Equilibria in the Gas Phase and in Solution Tautomeric Equilibria. Perhaps the hydroxypyridine system is among the most intensively studied tautomeric equilibria. Numerous experimental and theoretical data exist in the Previous ab initio calculations demonstrated that use of large basis sets and inclusion of electron correlation are necessary to correctly predict tautomerization e n e r g i e ~ . ~ . ~ ~
J. Phys. Chem., Vol. 98, No. 51, 1994 13775
Tautomeric Equilibria of 2- and 4-Hydroxypyridine
TABLE 1: Thermodynamical Results for Tautomerization of 2- and 4-Hydroxypyridinein the Gas Phasd 2- or 4-hydroxypyridine 2-pyridone 4-pyridone 0.4 8.1 M(‘Qf1) -0.5 1.9 hEo(6-3l+G(d)) hEO(MP4SDTQ16-3l+G(pd)) 1.o 3.2 0.3 0.3 AEy 1.3 3.4 AH298 0.5 -0.7 AS 1.1 3.7 AQg8 0.3 7 Uexpt) 3 AS(expt) 0.8 AG(expt) “Units are kcdmol for energies and cd(mo1 K) for entropies. 6-31+G(d) geometries are used in ab initio calculations.
-.
,
3.0
0.0
,
I
1
30.0
90.0
60.0
120.0
150.0
6.0
TABLE 2: Computed Relative Solvation Free Energies (AAG,-d of the Tautomeric Equilibria of HydroxypyridineIPyridone in Water and in Chloroform at 25 “C (kdmol) 2- or 4-hydroxypyridine 2-pyridone 4-pyridone QM/MM expP QM/MM expta -1.7 f 0.1 -1.6 -2.8 f 0.1
