6432
J. Phys. Chem. 1992,96, 6432-6439
Comparison of the Hybrid AMI/TIP3P and the OPLS Functions through Monte Carlo Simulations of Acetic Acid in Water Jiali Gao Department of Chemistry, State University of New York at Buffalo, Buffalo, New York 14214 (Received: March 24, 1992)
The hybrid quantum mechanical (QM) AMl/TIP3P and the empirical OPLS potential functions have been examined through statistical mechanical Monte Carlo simulations of acetic acid in water at 25 OC and 1 atm. It was found that the results obtained from these two approaches are in good accord, supporting the validity of the combined QM and molecular mechanical ( M M )approximation in condensed-phase simulations. In the AMI /TIP3P treatment, acetic acid was characterized by the valence electrons and the nucleus cores with the AM1 theory, while the TIP3P model was adopted for the solvent molecules. It was noted that the predicted hydrogen-bonding energies for bimolecular complexes using the AMl/TIP3P potential are in good accord with the ab initio 6-31G(d) results; however, hydrogen-bond distances in the QM hydrogen-bond donor complexes are about 0.2-0.4 A shorter than those predicted by the empirical function and ab initio calculations. This was further characterized by energy distribution and radical distribution functions obtained from the fluid simulations. In addition, Mulliken population analyses in the Monte Carlo simulations revealed details of the polarization effects on solute charge distribution, leading to large induced dipole moments in aqueous solution.
Introduction Computer simulations can provide valuable insights into the dynamics and structure of condensed-phase matters.’ A key element underlying the success of these calculations is the need for a set of intermolecular potential functions that can accurately describe molecular interactions in the system considered.2 Much effort has gone into the development of empirical potentials, and a number of excellent alternatives have been available for study of organic and biomolecular systemsa3 However, the fitting procedure to obtain good parameters is rather time-consuming and often difficult due to the lack of experimental data for structural and energetic comparisons. Furthermore, pairwise empirical potential functions are not appropriate for describing detailed polarization effects and chemical procases involving bond forming and breaking.Ib To overcome these problems, a hybrid quantum mechanical (QM) and molecular mechanical (MM) In this method, part of the system, approach has been e.g., the solute molecule in a dilute solution or the active site region in an enzyme, is explicitly treated quantum mechanically, while the “environmental” solvent molecules are approximated by empirical MM potentials. Thus, it is no longer necessary to parametrize the potential functions for the solute molecules, aside from the approximations involved in the quantum mechanical theory. Although hybrid QM/MM potentials have been extensively explored in energy minimization for determination of enzymatic reaction paths: application and testing in condensed-phase simulations are To provide a detailed analysis of the QM/MM potential function for liquid simulations, results from Monte Carlo simulations of the hydration of acetic acid using the hybrid AM1/TIP3Pg,bs and the OPLS9 models are presented and compared here. CH3COOH is a strong organic acid with a pK, value of 4.7. Thus, the protonated form only exists at low pH in aqueous solution. Carboxylate groups are important nucleophiles, and the basicities of the syn and anti lone pairs have an important implication in enzymatic reaction mechanisms.I0 In the gas phase, syn-acetic acid is about 6 kcal/mol more stable than the anti conformer from both ab initio MP3/6-311+G(pd) and AM1 calculations,*.l1while the difference is diminished in aqueous solution with a free energy of only 1 kcal/mol in favor of the syn.8JodJ2Both syn and anti structures are stable monomers in the gas phase and in aqueous solution, separated by a rotational barrier of 10.6 and 8.3 kcal/mol at the AMI/TIP3P level in the two phases, respectively.8 Wiberg’s ab initio calculations yielded a barrier height of 12.6 kcal/mol in the gas phase.’lb Experimentally, the barrier height for acetic acid was determined to be 10 kcal/mol from an IR study.I3 The anti/syn energy difference 0022-3654/92/2096-6432S03.00/0
does not appear to be available in the literature, whereas a Raman spectroscopic measurement predicted a value of 3.90 kcal/mol for formic acid.14
Computational Details Intermolecular Potential Fmctiow. The combined QM/MM potential function for fluid simulations is partitioned into a QM region consisting of the solute molecule and an MM region of solvent molmles.48.’ In order to compute the energies of the QM particles throughout the condensed-phase simulations, a computationally efficient method must be used. In the present study, the semiempirical Austin model 1 (AM1) theory15 developed by Dewar and co-workers was employed for the QM region.”* Extension to ab initio methods should be straightforward. The AM1 method can yield excellent results for many classes of systems that have been parametri~ed.’~ In many cases, its performance is comparable to that of the 3-21G and 4-31G level in ab initio c a l c u l a t i ~ n s . ~Thus, ~ J ~ for the quantum mechanical region, the solute molecule, acetic acid, is represented by the valence electrons and nucleus cores, while the restricted Hartree-Fock (HF) wave function for these particles is described by a single Slater determinant of all doubly occupied molecular orbitals (MO), ($JI6 = IJ/I.(1)$18(2)...$,B(2N))
(1) where UV is the total number of electrons, ci and 8 are the electron spin eigenfunctions, and the MO’s are linear combinations of a minimum basis set for each atom, (4J. In the MM area, the solvent molecules are approximated by interaction sites located on each nucleus. Consequently, the total potential energy of the system is given by the equation4gv7 (fr
Em, E M M + EQM + EQM/MM (2) Here, the terms on the right-hand side of eq 2 are respectively the total MM pair interaction energy for the solvent molecules, the H F energy of the QM region, and the interaction energy between the QM and MM particles. As usual, the intermolecular energy of a water dimer in the MM region is given as the sum of Coulomb interactions between all atomic pairs plus a Lennard-Jonesterm between the two oxygen atoms (eq 3).3 The threesite TIP3P model was adopted for water ona onb hE,b
=
i
J
(qtqj/rij) + 4%0[(a00/r00)’z- (‘J00/r00)61
(3) with experimental geometry held fixed throughout the simulation~.~’ The TIP3P model has been shown to be appropriate for 0 1992 American Chemical Society
The Journal of Physical Chemistry, Vol. 96, No. 15, 1992 6433
Monte Carlo Simulations of Acetic Acid in Water TABLE I: Lennrrd-Joaesand Core-Repulsion Parameters Used in the AMl/TIPJP Model OPLS AM 1 /TIP3P atom u, A c, kcal/mol u, A e, kcal/mol a PO Acetic Acid 0 3 .OO 0.170 2.10 0.150 C 3.75 0.105 3.50 0.080 =o 2.96 0.210 2.95 0.200 C(CH3) 3.91 0.160 3.50 0.080 e HO 0.0 0.0 0.80 0.100 HC 2.00 0.010 ~
0 H
3.1506 0.0
0.1521 0.0
Water 3.1506 0.0
0.1521 0.0
5.0" 5.0"
0.0" 0.0"
Reference 4.
studying hydration properties.18 On the other hand, the QM energies are computed via the HF-SCF procedure using the AM1 Hamiltonian augmented by a solute-solvent interaction term, HQMIMM, and are given as f o l l ~ w s ~ g , ~ EQM (@lfiQM(e)l@) - (@"(fiQM(e=o)l@") (4) and EQMIMM =
(5) where e is the C-O torsional angle, fiQM(6) 1s the electronic Hamiltonian for the QM acetic acid molecule, HQM/MM(e) is the QM/MM interaction Hamiltonian, and @" and @ are respectively the wave functions of the QM atoms in the gas phase and in aqueous solution. Note that eq 4 can further be dm-mposed into an elptron distortion energy term ( E & ~= (@IH(O)QMl@) ( @ o l ~ ( d ) Q M l @ o ) ) and a gas-phase torsional energy term (Etor= (@"IH(~)QMl@") - ( @ " ~ H ( ~ = O ) Q M ~ @ " ) ) . The SUllllllatiOn in eq 5 features the van der Waals interactions between the QM and MM nuclei, whose need has been discussed previously.4g It should be pointed out, however, that the Lennard-Jones parameters are the only adjustable variables in the hybrid AM1 /TIP3P model, and their values have been found to be not sensitive in energy and geometry computations4g except when the QM atoms are hydrogen-bond donors. In the latter case, optimization of the Lennard-Jones parameters for the heteroatom and hydrogen is necesary. The resulting parameters along with the a and po values for water, which were taken from ref 4g, are given in Table I. Monte Carlo Simulations. Statistical mechanical Monte Carlo calculations were carried out for systems of one acetic acid and 216 water molecules in a cubic cell with periodic boundary conditions at 25 "C and 1 atm. Standard procedures were adopted including Metropolis sampling and the isothermal-isobaric ensemble (NPT).2 To facilitate the statistics near the solute molecule, Owicki-Scheraga preferential sampling technique was used TABLE II: Computed Complexation Energies (kcal/mol)" AM 1 /TIP3P comDlex Rnu dl -AE
I I1 I11 IV anti-AcOH-H20 I I1 I11 IV V
1.92 1.94 1.84 1.53
128 117 139 140
5.46 5.69 3.58 7.1 1
with l / ( 9 + c) weighting, where c = 150 Two systems were considered in this study with fixed anti (e = 180") and syn (e = 0")conformations for the C-O rotation,8 while the C-C torsion was sampled via dihedral variations during the Monte Carlo calculations.20 In addition, the solute geometry was optimized whenever there is a solute move in thc simulation, though the solutmlvent interaction Hamiltonian (HQMIMM) was not included in the optimization. So, the solute structure corresponds to that in the gas phase. Geometry optimizations were needed because relaxation of the bond lengths and angles in acetic acid should be allowed in the dihedral variation.8 The solute-solvent interaction energy was computed by a single-point SCF calculation with the inclusion of H Q M I M M . ~ The intermolecular interactions were feathered to zero at a spherical cutoff distance between 8.5 and 9.0 A based roughly on the center-of-mass separation. New configurations were generated by randomly selecting a molecule, translating it in all three Cartesian directions, rotating it along a randomly chosen axis, and varying the internal rotation. For the dilute aqueous solution of CH3COOH an acceptance rate of ca. 40% was maintained by using ranges of f0.13 A and f13" for the solvent motions and fO.10 A and *loo for CH3COOH. Volume changes were restricted to within f l O O A3 on every 1250 configurations. For comparison, similar computations were repeated using the empirical OPLS potential function for CH3COOH9and the TIP3P model for water.17 These functions have been derived through gas-phase ab initio computations of hydrogen-bonding complexes and pure liquid simulation^.^ Here, acetic acid is described by five interaction sites, treating the methyl group as a single unit. In each case 1 X lo6 configurations were discarded for equilibrating the system followed by an additional 1.5 X lo6 configurations for averaging of the reported properties. All computations were carried out with the BOSS programZomodified to include the QM/MM energy terms,' which were evaluated using the MOPAC package.2i The QM/MM calculations required 2.5 days, while the pure molecular mechanical potential needed 20 h on a SUN SPARCstation 2 computer in our laboratory.
Results and Discussion Bimolecular Complexes. To assess the validity of the QMAMl/TIP3P model for simulations of the hydration of acetic acid in water, hydrogen-bondingcomplexes of CH3COOH with water were investigated. Intermolecular geometry optimizations were executed using a program developed in our laboratory, treating CH3COOH quantum mechanically with the AM1 method and water molecular mechanically. The results are summarized in Figure 1 and Table I1 and are compared with the results predicted by the OPLS and ab initio 6-31G(d) calculation^.^ The CH3COOH structure was fixed in the equilibrium geometry at the corresponding level of theory, while experimental configuration was used for water. The computed results from the AMl/TIP3P model in Table I1 are in reasonable accord with the data from the 6-31G(d) calculations. The agreement in hydrogen-bonding energy is
Rnu
1.73 1.81 1.84 1.66
OPLS/TIP3P dl -AE 103 132 144 150
7.01 5.39 3.18 9.63
2.06 1.94 1.91 1.56
98 118 105 156
5.38 5.73 6.19 7.72
1.81 1.78 1.92 1.88
137 137 99 162
5.77 6.44 6.47 9.18
1.54
169
7.68
1.71
164
11.07
Rnu
6-31G(d) dl
2.01 2.05 2.20 1.91b
122 120 142 1 16b
1.88
111
5.52 5.40 2.38 7.36b 8.45
2.1 1 2.03 2.35 1 .96b 1.95 1.94
121 120 98 14ob 136 175
5.18 5.85 4.86 7.84b 8.46 8.01
#Distance and angle are given in angstroms and degrees. bComputed using the 6-31+G(d) basis set and geometry.
-AE
6434 The Journal of Physical Chemistry, Vol. 96, No. 15, 1992
1’
+
Q
P
A Figure 1. Geometrical arrangements for the CH,COOH-water complexes. particularly good for the CH3COOH acceptor complexes (I and 11), whereas when CH3COOHacts as a hydrogen-bond donor (IV) the hydrogen-bond distances have to be shortened in order to obtain good energetic results. This difficulty was readily attributed to the much smaller partial charge on the hydroxyl hydrogen according to Mulliken population analysis in the AM 1 model than that obtained using the 6-31G(d) basis set and that optimized in the OPLS function (Table 111). Hydrogen-bonding interactions in water are somewhat strengthened in the AMl/TIP3P model due to electric polarization of the solute (see below). This trend also exists in other QM hydrogen-bond donors and perhaps is a major shortcoming of the AMl/TIP3P potential. In our previous study of the relative aqueous basicities of the carboxylate lone pairs, the OPLS parameters were used for the Lennard-Jones terms in eq 5 for acetic acid and the interaction energy for structure IV was predicted to be ca. 2-3 kcal/mol smaller than the 6-31G(d) value for both syn and anti conformations,8while hydrogen-bond distances were in good accord with the ab initio calculations. In that study, the major concem was whether there is a net bias of hydration of either syn or anti conformer, and the results indicate that the AMl/TIP3P model with the OPLS parameters for the solute can provide reasonable description of the relative interaction energies for syn- and anti-acetic acid with waters8 We have now optimized the Lennard-Jones parameters
Gao in eq 5 such that hydrogen-bonding energies for both QM hydrogen-bond donor and acceptor complexes can be correctly predicted, but at the sacrifice of some geometrical accuracy for donor complexes. It should be mentioned that the Lennard-Jones parameters do not affect signifhntly the hydrogen-bonding energy and geometry if the QM atoms are hydrogen-bond acceptors4g (see structure 111); they only become sensitive when the QM particles are hydrogen-bond donors. Structure IV was predicted by all three models to have a stronger interaction energy in either syn or anti cbnfiguration than hydrogen-bond acceptor complexes. Using the OPLS function, structure V was predicted to be the lowest energy form (-1 1.1 kcal/mol) with a bent hydrogen bond due to an unfavorable van der Waals contact between the methyl group and water oxygen. This represents an overestimate of ca. 3 kcal/mol compared to the ab initio number. Note that the COH angle in the OPLS function was kept the same in both anti and syn conformers. Increasing the COH angle in anti-acetic acid by 10’ gives the same complexation energy, but with a roughly linear hydrogen bond. Structural features predicted by the AMl/TIP3P model nicely mirror the trends from the 6-31G(d) and OPLS optimizations (Table 11) except for structure IV. The computed intermolecular distances for acceptor complexes at the AMl/TIP3P level are shorter than the corresponding ab initio values by 0.1-0.15 A but are ca. 0.1 A longer than the OPLS estimates. The shorter hydrogen-bond distances in the OPLS function are needed in order to produce good agreement between the calculated and experimental liquid pr~perties.~ The intermediate interaction distance in the AMl/TIP3P model is also reasonable since the HF-SCF calculations have only included the polarization effects into the solute molecule through the solute (QM)-solvent (MM) interaction Hamiltonian. A complete description of both the gas-phase and condensed-phase properties requires a full consideration of the solute and solvent polarizations.22 Atomic Charges. A major issue in developing parameters for empirical potential functions is to obtain the atomic charge^.^ Typically, they are taken from Mulliken population analysis using large basis set and adjusted to give reasonable results for gas-phase complexes and for the pure liquid^.^ Alternatively, they can be fitted to electrmtatic potential surface generated from HF wave functions.23 The first approach provides the best results for condensed-phase simulations but requires a laborious iterative optimization.” The latter approximation, on the other hand, yields atomic charges that can reproduce the gas-phase electrostatic properties. However, electron clouds are polarized in solution owing to interactions with the surrounding solvent molecules. A direct consequence of the polarization effects on the atomic charges in solution is an increase in dipole moment by ca. 10-2095 over the corresponding gas-phase value^.^^" Since the solvent effect is directly incorporated into the H F theory in the present AM 1/TIP3P approach,+’ analysis of the electron population calculations throughout the Monte Carlo simulations could provide valuable information concerning the reorganization of the solute electron clouds in water. It should be noted that in our current approach induced solvent polarization is not explicitly included, although it has been implicitly averaged into the effective potential function for water.” Table I11 recorded the average atomic charges for syn- and anti-CH3COOH in water along with the OPLS partial charges and those obtained from AM1 and 6-31G(d) calculations (gas phase). Computation of aqueous atomic charges through QM/MM simulations has been reported by Bash et al. in their study of an SN2 reaction in water,6 although previous investigations in the context of self-consistent reaction field (SCRF) theory treating the solvent as a dielectric continuum medium could also give similar results.24 Obviously, the latter cannot provide any detailed information on molecular interactions in solution. It is first interesting to compare the partial charges of acetic acid in the gas phase and in water obtained from Mulliken pop ulation analyses using the AMl/TIP3P model. Clearly, there is a substantial increase in the magnitude of the net charges on all atoms. The largest variation from the gas-phase value occurs
The Journal of Physical Chemistry, Vol. 96, No. 15, 1992 6435
Monte Carlo Simulations of Acetic Acid in Water TABLE III: Computed Partial Atomic Charges (e) atom AMl/TIP3P (as) AM1 (gas) 6-31G(d) (gas) OPLS syn-Acetic Acid 4 -0.437 i 0.004 -0.361 -0.551 -0.50 0 -0.380 i 0.004 -0.321 -0.689 -0.58 C 0.349 i 0.002 0.306 0.724 0.55 C(CH3) -0.234 i 0.001 -0.217 -0.571 0.08 HO 0.307 i 0.002 0.243 0.465 0.45 HC' 0.132 f 0.002 0.117 0.207 A D 2.21 i 0.04 1.89 1.79 1.51 4 0 C C(CH3) HO HC" P, D
anti-Acetic Acid -0.426 i 0.002 -0.296 -0.364 i 0.005 -0.293 0.356 f 0.001 0.294 -0.270 i 0.001 -0.268 0.306 i 0.003 0.221 0.132 i 0.001 0.1 14 5.77 f 0.02 4.39
a Averaged
-0.519 -0.690 0.755 -0.625 0.463 0.206 4.73
-0.50 -0.58 0.55 0.08 0.45
syn-Acetic Acid
5.05
over all three methyl hydrogen atoms.
TABLE IV: Computed Energies at 25 OC and 1 atm (kcal/mol) species syn-CH3COOH (OPLS) anti-CH3COOH (OPLS) SYWCH~COOH (AMl/TIP3P) anri-CH$OOH (AMl/TIP3P)
Ecdi,,
4, E,or -22.2 i 0.5 0.0 -34.1 i 0.4 5.0 1.8 i 0.1 -1.8 i 0.1 -22.1 i 0.7 0.0 Epot
3.8 i 0.1 -4.2 i 0.1 -31.8 i 0.8 6.0
= (@I& + AQM/MM(@)- (@'I& + I ? Q M / M M ) @ ' ) , where a' a n i k r e wave functions of the solute in the gas phase and in aqueous solution.
at the carbonyl oxygen with deviations of 0.13 and 0.08 e for the anti and syn conformers. The overall charge polarization is greater for anti-CH3COOH than the syn monomer in water, due to a stronger solvation effect (Table IV). This is also indicated by the required to polarize the greater electric distortion energy E\,, electron clouds in anti-CH,COOH (Table IV). To illustrate the polarization effects, Figure 2 displays the electron density difference (EDD)plots for acetic acid in water at an instantaneous configuration during the simulation relative to that in the gas phase. (The surrounding solvent configurations are shown in Figure 13.) Gains and losses in electron density are indicated by the solid and dashed curves.25 Note that the charge distributions in Table I11 are essentially the same for the two conformers in water, while differences up to 0.05 e exist in the gas phase from both AM1 and 6-31G(d) computations. The average dipole moments in water are 5.77 f 0.02 D for anti-CH3COOHand 2.21 f 0.04 D for syn-CH3COOH, which represent increases of 1.38 and 0.32 D over the gas-phase data. These may be compared to values of 5.05 and 1.51 D based on the OPLS charge^.^ The experimentaldipole moment is 1.74 D for syn-acetic acid in the gas which was reproduced by both AM1 and 6-31G(d) calculations (Table 111). Finally, a few comments should be made. The ab initio 63 1G(d) charges are much larger than both the ones optimized in the OPLS function and those computed by the AMl/TIP3P simulations, whereas the gas-phase AM1 charges are too small to be useful in empirical potentials for liquid simulations. Of course, the Mulliken population only provides a rough estimate of the atomic charges and certainly depend strongly on the basis set used and on the optimzation procedure adopted by semiempirical quantum mechanical methods. Although the computed AMl/TIP3P charges for acetic acid in water are still much smaller than the OPLS data, the net charges are significantlygreater than the gas-phase values. Furthermore, it is not unreasonable to observe good results for bimolecular complexes and aqueous solution using the hybrid AMl/TIP3P method because the electron clouds are delocalized and are "closer* to the solvent molecules than the static atomic charges projected on the nuclei from the Mulliken population analysi~.~'The results presented in Table I11 provide a good assessment of the polarization effects on the
anti-Acetic Acid
F w 2. Electron density difference (EDD) plots for acetic acid in water and in the gas phase. Dotted contours represent regions where there is a depletion of electron density, and solid curves indicate areas where there is a gain in electron density, on transferring the solute from the gas phase into aqueous solution. The surrounding solvent molecules are shown in Figure- 13. Bonding Energy Distributions
-50
-42
-34
-26
-18
-10
Bonding Energy (kcal/mol)
Figure 3. Distributions of total solutesolvent energies obtained using the OPLS function (solid line) and the AMl/TIP3P model (dashed line). Units for the ordinate are mole percent per kcal/mol.
solute atomic charges in solution and a quantitative measure of the induced dipole moment. Note that the present study employed a semiempirical MO method in the condensed-phase simulations due to the limitation of computer speed;it is convenient to extent the QM/MM approach to ab initio MO calculations in the future. Energy Distribution Functions. Figures 3-5 revealed some details of the solutesolvent interactions regarding the total solute-water bonding energy (Figure 3) and the individual acetic acid-water interaction energy (Figures 4 and 5). The averages obtained from the bonding energy distributions in Figure 3 are
6436 The Journal of Physical Chemistry, Vol. 96, No. 15, 1992
Gao
Energy Pair Distributions
3 .O
1.5
I
O(syn)-H Radial Distribution Functions
O(anti)-H Radial Distribution Functions
-15
-10 -5 0 Interaction Energy
5
h
E5 80
Figure 4. Distributions of individual solute-solvent interaction energies (kcal/mol) for syn-acetic acid in water. Units for the ordinate are number of solvent molecules per kcal/mol. Energy Pair Distributions
3.0
;
2
1
3
7
8
Figure 6. Carbonyl oxygen-water hydrogen radial distribution functions for syn-acetic acid (top) and anti-acetic acid (bottom) obtained using the OPLS function (solid l i e ) and the AMl/TIP3P model (dashed line). This convention is adopted throughout.
I
I
-10 -5 0 Interaction Energy
6
5
R (A)
I
-15
4
I
O(syn)-0 Radial Distribution Functions
5
Figure 5. Distributions of individual solutesolvent interaction energies (kcal/mol) for anti-acetic acid in water.
the E, values listed in Table IV, which is the total solutesolvent interaction energy. In Figure 3, there is a continuous distribution of the bonding energy in each case. The computed ranges for syn-acetic acid from the hybrid AMl/TIF'3P model and the OPLS function are nearly identical with averages of -22.1 f 0.7 and -22.2 f 0.5 kcal/mol, providing strong support to the AM1/ TIP3P approach. However, the range for anti-acetic acid obtained from the OPLS simulations is stronger than the corresponding results from the AM1 /TIP3P computations. The average total bonding energies are -31.8 f 0.8 and -34.1 f 0.4 kcal/mol using the AMl/TIP3P and OPLS functions, respectively (Table IV). This difference is consistent with predictions for bimolecular complexes (Table 11) where the hydrogen-bonding energy for structure V in anti-acetic acid was somewhat overestimated using the OPLS functions. The energy pair distribution functions listed in Figures 4 and 5 provide details of solutewater interaction energies. The lowenergy bands are from the hydrogen-bonding interactions between acetic acid and water, while the central peaks centered at 0 kcal/mol represent the more distant water molecules. For both the OPLS and the AMl/TIP3P potentials, there are clearly identifiable peaks in both systems at the low-energy region, corresponding to donor and acceptor hydrogen bonds. The first band can be assigned to the hydrogen-bonding complex with donor CH,COOH. Note that the best interaction energy for the anti structure is lower than the syn conformation (also see Table 11), while the bands near -5 kcal/mol are due to hydrogen bonds to oxygen atoms in acetic acid, which showed much less structure using the OPLS potential than the AMl/TIP3P model. The OPLS distribution of bimolecular interaction energies only showed a small shoulder in these plots. The best interaction energies recorded in water are -1 1.0 and -9.3 kcal/mol for anti and syn conformations, respectively, using the AM 1/TIP3P potential, whereas they are -11.8 and -10.2 kcal/mol from the OPLS function. The energy pair distributions obtained using the OPLS function are consistent with the bimolecular interactions listed in Table 11, while the enhanced hydrogen-bonding interactions
1'
0.0
-
1.5
O(anti)-0 Radial Distribution Functions !I
I
0.0
r
1
2
I
3
4
6
5
8
7
R (A)
Figure 7. Carbonyl oxygen-water oxygen rdfs. 1.5
Oh(syn)-H Radial Distribution Functions
1.0
-
0.5
-
h
E5 00
0.0
1
I II
Oh(anti)-HRadial Distribution Functions
h
B M
1
2
3
5
4
6
I
8
R (A)
Figure 8. Hydroxyl oxygen-water hydrogen rdfs.
observed with the AMl/TIP3P model reflect the polarization effects on the solute (see above).
The Journal of Physical Chemistry, Vol. 96, No. 15, 1992 6437
Monte Carlo Simulations of Acetic Acid in Water
1.5
Oh(syn)-0 Radial Distribution Functions
H(syn)-H Radial Distribution Functions
1.5 r
1.0
-
0.5
-
EM 5
,-a_\
1.o
\
I
1.5
17
!I
I
, ,
II
0.0
1
H(anti)-H Radial Distribution Functions
Oh(anti)-0 Radial Distribution Functions
-
2
1
3
4
5
6
R (A) Figure 11. Hydroxyl hydrogen-water hydrogen rdfs.
Figure 9. Hydroxyl oxygen-water oxygen rdfs. 3.0
I
i
2.0
0.0
0.0
I, 1
//I
'
C(syn)-0 Radial Distribution Functions
H(syn)-0 Radial Distribution Functions
\-y I
1
I
I
I
1
2
3
4
5
6
I
8
R (A)
Figure 10. Hydroxyl hydrogen-water oxygen rdfs.
Radial Distribution Functions. Seven radial distribution functions (rdfs) computed with both the AMl/TIP3P and OPLS functions are compared in Figures 6-12 for syn- and anti-acetic acid. In these figures the first atom for an xy distribution, g&), refers to an atom of the solute and the second atom is either the oxygen or the hydrogen of water. For the solute atoms, 0 and Oh are the carbonyl and hydroxyl oxygens, and H is the hydroxyl hydrogen. The error range in these computations can be approximated by one half of the bin size (0.04-0.05 A) in rdf computations. Hydrogen-bonding interactions are clearly indicated by the strong first peaks in the 0-H rdfs in Figure 6. The maxima are at 1.94 and 1.86 A for the syn and anti conformers using the OPLS function, which may be compared to a value of 2.00 A for both structures in AMl/TIP3P models. Integration to 2.5 A yields a coordination number of 1.5 and 1.6 for syn- and anti-acetic acid in the empirical approximation, while the results are 1.9 and 2.1 waters, respectively, in the QM/MM approach. Hydrogenbonding interactions appeared to be more favorable at the carbonyl group at the AMlfTIP3P level, which was echoed in the 0-0 rdfs with higher and more structured first peaks than the OPLS results (Figure 7). This follows the greater polarization effect and atomic charge developed at the carbonyl oxygen in water. In Figures 8 and 9, the Oh-H and Oh-O rdfs in the AM1/ TIP3P model are 0.1-0.2 A ahead of those predicted by the OPLS potential. In contrast to the previous 0-H pair, the interaction
-
I
C(anti)-0 Radial Distribution Functions
I
8
6438 The Journal of Physical Chemistry, Vol. 96, No. 15. 1992
Gao
BB Figure 13. Stereo plots of configurations from the simulations of syn-acetic acid (top) and anti-acetic acid (bottom) in water using the AMl/TIP3P
function. Water molecules that are more than 5 A from any solute atoms have been removed for clarity.
shown in Figure 13. For clarity, water molecules that are more than 5 A from any atoms of the solute have been deleted, leaving the viewing windows slightly larger than the first solvation layer. The ubiquitous hydrogen bonding to acetic acid is apparent for both structures. In both cases, the plots reinforce the results noted above. At the carbonyl center, there are typically two hydrogen bonds between the carbonyl oxygen and the hydrogens of water, while the two hydrogen bonds at the hydroxyl group are divided into one donor and one acceptor hydrogen bond. Qualitatively, the hydrogen-bonding patterns and polarization effects are nicely mirrored in the EDD plots discussed earlier (Figure 2). Conclusion The present study provides a direct comparison of the hybrid QM-AM 1/TIP3P and the empirical OPLS potential functions through Monte Carlo simulations of the dilute aqueous solution of acetic acid. In the AMl/TIP3P approach, acetic acid was treated quantum mechanically using the AM1 theory throughout the simulations,while the solvent was approximated by the TIP3P model. The computed structural and energetic results were found to be in good agreement with the OPLS function. For acetic acid-water dimer complexes, the hydrogen-bond donating ability in CH3COOH predicted by the AMl/TIP3P model was about 2-3 kcal/mol weaker than the ab initio 6-31G(d) results if the OPLS Lennard-Jones parameters were used for acetic acid, while the ab initio results for CH3COOHacceptor complexes were well reproduced by the AMl/TIP3P and the OPLS calculations.8 Reoptimization of these parameters was made, and good agreement between the computed interaction energies (AM1 /TIP3P vs 6-31G(d)) for the bimolecular complexes was obtained, however, hydrogen-bond distances between the carboxyl hydrogen and water oxygen were ca. 0.4 A shorter in the AMl/TIP3P model than the ab initio predictions. Further investigation of the Mulliken population analyses in the Monte Carlo simulations revealed the polarization effects on solute charge distributions,and a significant increase in the net atomic charges was observed in aqueous solution, leading to large induced dipole moments. Structural features were also in reasonable a w r d between the AM1 /TIP3P and OPLS potentials as shown by radial distribution functions.
Since the OPLS function has been well tested and applied to many system^,^ the agreement in the computed results reported here indicates that the AM 1/TIP3P approach is a reliable procedure for studying condensed-phase matters. Acknowledgment. Gratitude is expressed to the donors of the Petroleum Research Fund, administered by the American Chemical Society, for partial support and to the Julian Park Publication Fund for covering page charges. Registry No. Acetic acid, 64-19-7; water, 7732-18-5. References and Notes (1) (a) Karplus, M.; Petsko, G. A. Nature 1990, 347, 631. (b) McCammon, J. A.; Harvey, S.C. Dynamics of Proteins and Nucleic Acids; Cambridge University Press: Cambridge, 1987. (2) Allen, M. P.; Tidesley, D. J. Computer Simulations of Liquids; Oxford University Press: London, 1987. (3) (a) Jorgensen, W. L.; Tirado-Rives, J. J . Am. Chem. Soc. 1988,110, 1657. (bl Weiner. S. J.; Kollman, P. A.; Case. D. A.; Sinah. U. C.: Ghio. C.; Alagona, G.; Profeta, S.Jr.; Weiner, P. J. Am. Chem. Szc. 1984, 106, 765. (c) Brooks,B. R.; Bruccoleri, R. E.; Olafson, E. D.; States, D. J.; Swaminathan, S.; Karplus, M. J . Comput. Chem. 1983, 4, 187. (4) (a) Warshel, A.; Karplus, M. J. Am. Chem. SOC. 1972,94,5612. (b) Wang, I. S.;Karplus, J. J . Am. Chem. SOC. 1973, 95, 8160. (c) McCreery, J. H.; Christoffersen, R. E.; Hall, G. C. J . Am. Chem. Soc. 1976, 98, 7191. (d) Thole, B. T.; van Duijnen, P. Th. Chem. Phys. 1982, 71,211. (e) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7,718. (f) Tapia, 0.;Colonna, F.; Angyan, J. G. J . Chim. Phys. 1990.87, 875. (9) Field, M. J.; Bash, P. A.; Karplus, M. J . Comput. Chem. 1990, 11, 700. (5) (a) Alagona, G.; Desmeules, P.; Ghio, C.;Kollman, P. A. J. Am. Chem. Soc. 1984, 106, 3632. (b) Weiner, S.J.; Seibel, G. L.; Kollman, P. A. Proc. Narl. Acad. Sci. U.S.A. 1986,83,649.(c) Hwang, J.-K.; King, G.; Greighton, S.;Warshel, A. J . Am. Chem. Soc. 1988, 110,5297. (d) Rullmann, J. A. C.; Bellido, M. N.; van Duijnen, P. T. J . Mol. Biol. 1989, 206, 101. (e) Waszkowycz, E.; Hillier, I. H.; Gensmantel. N.; Payling, D. W. J. Chem. Soc., Perkin Trans. 2 1991, 225. ( f ) Bash, P. A.; Field, M. J.; Davenport, R. C.; Petsko, G. A.; Ringe, D.; Karplus, M. Biochemistry 1991, 30, 5826. (6) Bash, P. A.; Field, M. J.; Karplus, M. J . Am. Chem. Soc. 1987, 109, 8192. (7) Gao, J. J . Phys. Chem. 1992, 96, 537. (8) Gao, J.; Pavelites, J. J. J. Am. Chcm. Soc. 1992, 114, 1912. (9) (a) Briggs, J. M.; Nguyen, T. E.; Jorgensen, W. L.J . Phys. Chem. 1991,95,3315. (b) Jorgensen, W. L.; Briggs, J. M.; Contreras, M. L. J. Phys. Chem. 1990, 94, 1683. (10) (a) Gandour, R. D. Bioorg. Chem. 1981,10, 169. (b) Kirby, A. J. In The Anomeric Effect and Related Stereoelectronic Effects at Oxygen;
J. Phys. Chem. 1992, 96,6439-6441 Springer-Verlag: Berlin, 1983; p 135. (c) Zimmerman, S. C.; Korthals, J. S.;Cramer, K. D. Tetrahedron 1991,47,2649. (d) Tadayoni, B. M.; Rebek, J. Jr. Bioorg. Med. Chem. Lett. 1991. I , 13. (11) (a) Wiberg, K. B.; Laidig, K. E. J. Am. Chem. Soc. 1987,109, 5935. (b) Wiberg, K. B.; Laidig, K. E. J . Am. Chem. Soc. 1988,110, 1872. (c) Li, Y.; Houk, K. N. J. Am. Chem. Soc. 1989,111,4505. (12) Allen, F. H.; Kirby, A. J. J. Am. Chem. Soc. 1991, 113, 8829. (13) Miyazawa, T. Bull. Chem. Soc. Jpn. 1961, 34,691. (14) Bertie, J. E.; Michaelian, K. H. J . Chem. Phys. 1982, 76, 886. (15) (a) Dewar, M. J. S.;Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J . Am. Chem. SOC.1985, 107, 3902. (b) Dewar, M. J. S.;Healy, E. F.; Holder, A. J.; Yuan, Y. C. J. Compur. Chem. 1990, 1 1 , 541. (c) Stewart, J. J. P. J . Compur. Chem. 1990, I I , 543. (16) Pople, J. A.; Beveridge, D. L. Approximate Molecular Orbital Theory: McGraw-Hill: New York. 1970. (17) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (18) For example; Gao. J.; Kuczera, K.; Tidor, B.; Karplus, M. Science 1989, 244, 1069. . (19) (a) Owicki, J. C.; Scheraga, H. A. Chem. Phys. Lett. 1977,47,600. (b) Owicki, J. C. ACS Symp. Ser. 1978, No. 86, 159. (c) Jorgensen,W. L.; Bigot, B.; Chandrasekhar, J. J. Am. Chem. Soc. 1982, 104, 4584. (20) Jorgensen, W. L. BOSS,Version 3.0; Department of Chemistry, Yale University, 1991. (21) Stewart, J. J. P. MOPAC, Version 5; Quantum Chemistry Program Exchange 455, 1986; Vol. 6, No. 391. . I
6439
(22) Caldwell, J.; Dang, L. X.;Kollman, P. A. J . Am. Chem. Soc. 1990, 112,9144. (23) (a) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1984.5, 129. (b) Cieplak, P.; Kollman, P. J. Comput. Chem. 1991, 12, 1232. (24) Two recent papers have compiled extensive lists of references of previous dielectric continuum SCRF computations. (a) Wong, M. W.; Frisch, M. J.; Wiberg, K. B. J . Am. Chem. SOC.1991, 113,4776. (b) Cramer, C. J.; Truhlar, D. G. J. Am. Chem. Soc. 1991,113,8305. I thank one referee for bringing us attention of the work described in ref 4f using a similar approach as the present QM/MM model. In that study, atomic charges and dipole moment of water were computed with various basis sets for 59 structures of one QM water in a bath of water molecules represented by polarizable dipoles. (25) (a) Jorgensen, W. L.; Severance, D. A modified version Psi/88; Department of Chemistry, Yale University. I thank Dan Severance for advice on the EDD plots. (b) Morrison, H.; Jorgensen, W. L.; Bigot, B.; Severance. D.; Munoz-Sola, Y.; Strommen, R.; Pandey. B. J . Chem. Educ. 1985,62,298. The following is quoted from (b): Because these are contours of varying density, no quantitative conclusions can be drawn from a delta plot; only a qualitative picture of gains and losses is available without a more extensive "population analysis". An attempt to use contour size as an indication of absolute magnitude is akin to attempting to estimate the relative masses of a mushroom and a lead ball by inspection. (26) Weast, R. C., Ed. CRC Handbook of Chemistry and Physics; CRC Press: Boca Raton, FL, 1988; pp E55-58. (27) Mulliken, R. S.J . Chem. Phys. 1955, 23, 1833.
Effect of Temperature on the Structural Properties of the Poly(ethy1ene oxlde)/Poly(propylene oxlde)/Poly( ethylene oxide) Triblock Copolymer Studled wlth Time-Correlated Fluorescence Problng Techniques Panagiotis Lianos' School of Engineering, Physics Section, University of Patras, P.O. Box 5054, 26004 Patras, Greece
and Wyn Brown Department of Physical Chemistry, University of Uppsala. Box 532, 75121 Uppsala, Sweden (Received: November 22, 1991; In Final Form: March 16, 1992)
The effect of temperature between 10 and 50 OC on the structural properties of a poly(ethy1ene oxide)/poly(propylene oxide)/poly(ethylene oxide) triblock copolymer has been studied with time-correlated fluorescence probing techniques using pyrene. Two different models have teen used for the analysis of the fluorescencedecay profiles, one probing micellar structures and the other probing percolating structures. Between 15 and 25 "Cwe have detected a pronounced transition from a micellar to a solidlike gel phase, while with a further increase of temperature our results indicate a possible second phase transition.
Introduction Much interest has been recently shown in the aggregation properties of block copolymers in water.14 In the present work we use excimer pyrene fluorescence probing to study the effect of temperature, in the range 10-50 O C , on the structural properties of the triblock copolymer poly(ethy1ene oxide)/poly(propylene oxide)/poly(ethylene oxide), abbreviated here as Pluronics P-85. The monomeric molecular weight is 4500, of which 50%belongs to the PEO (poly(ethy1ene oxide)) moiety.' P-85 forms micelles in aqueous solutions' which, in analogy to other micelleforming block copolymer^,^ are presumably composed of a core of water-insoluble PPO (poly(propy1ene oxide)) blocks with a swollen shell of PEO end blocks. It has been previously found' that temperature affects micelle formation and, generally, the aggregation properties of P-85. It is this phenomenon that we are scrutinizing now by employing time-correlated (and assisting steady-state) fluorescence probing techniques, which have been proven to be a powerful tool in the investigation of the structural properties of organized a s s e m b l i e ~ . ~When ' ~ pyrene is used as fluorescent probe, as in the present case, the monomer pyrene fluorescence decay profiles in the presence of excimers in the dispersed phase are usually analyzed, offering information on the
structure and the dynamics of the organized phase itself. Supporting information is also obtained by analysis of the fluorescence spectra of both monomer and excimer pyrene.'"-I6 This work also consists of analyses of pyrene fluorescence spectra and pyrene fluorescence decay profiles. Experimental Section The triblock copolymer P-85 was obtained from Serva AG, Heidelberg, Germany, and used without further purification. It has a declared molecular weight of 4500, 2200 for the PEO component and 2300 for PPO (E,P E,, where x = 25, y = 40). In the present work we have used sofutions containing 18.7 g/cL of P-85 in highly purified water. The above aqueous solution easily solubilizes pyrene at all concentrations used in the present work. In fact the probe concentration varied from 10 FM up to 1 mM. Pyrene was purchased from Fluka (puriss.) and was used without further purification. Fluorescence decay measurements were recorded with the time-correlated photon-counting technique using a specially constructed hydrogen flash, Ortec and Schlumberger-Enertec Electronics, and a nucleus multichannel scaler card with an IBM-PC. An interference ftlter was used for excitation (337 nm),
0022-3654/92/2096-6439$03.00/00 1992 American Chemical Society
