J. Phys. Chem. 1994, 98, 10465-10476
10465
The Nature of Ks/Crown Ether Interactions: A Hybrid Quantum Mechanical-Molecular Mechanical Study Mark A. Thompson,*Eric D. Glendening, and David Feller Molecular Science Research Center, Pacific Northwest Laboratory, Richland, Washington 99352 Received: May 16, 1994; In Final Form: July 18, 1994@
We present a hybrid quantum mechanicallmolecular mechanical (QM/MM) molecular dynamics study of dimethyl ether (DME) and 18-crown-6 ( 1 8 ~ 6 )interacting with K+. The QM/MM method employs the semiempirical AM1 method to describe the ethers, the MM parametrization of Dang for K+, and the MM SPC/e model for H2O. We parametrize the interaction Hamiltonian to the binding energies and optimized geometries for K+/DME using a b initio HF and MP2/6-3 1+G* results. The resulting Q W model describes the polarization response of both free DME and K+-complexed DME well. The QM/MM model gives good agreement with the experimental and ab initio structures for K+/18c6. We calculate gas-phase K+/18c6 binding energies of -70.2 and -72.0 kcallmol with the QM/MM and MP2/6-31+Gh (CP corrected) methods, respectively. Our simulation results for Kf/18c6 in H20 show that the most probable K+/18c6 center-ofmass displacement is 0.25 A, in marked contrast to previous molecular dynamics results of Dang and Kollman. Our result is consistent with K+ having an optimal “fit” for the cavity of 18c6. Still, we find that K+ retains significant solvent accessibility coordinating two H20 molecules, on average, in the K+/l 8c6 simulation. The simulation average polarization energy for 18c6 interacting with both K+ and the H20 solvent is -14.9 kcal/mol, which is 17% of the total electrostatic interaction energy. This result underscores the potential importance of QM in describing the solution chemistry of ion-macrocycle interactions. Our study is the first simulation of crown ethers that explicitly incorporates QM in the force field.
useful for specific ion recognition in the development of sensors. Crown ethers have drawn much experimental and theoretical There have been numerous theoretical studies of crown ethers interest since they were first described by Pederson in 1 9 6 7 . ’ ~ ~ that have provided new insights into crownkation interactions. These macrocyclic polyethers show a remarkable range of These include ab i n i t i ~ , ~ ~emiempirical,’~ -’~ and molecular specificity for a wide variety of cations that depends, in part, mechanics (MM).14-23 In a seminal study, Wipff el a1.l’ on the size of the ether, the type of donor atoms (e.g., oxygen, attributed the preferential binding of 18-crown-6 (18~6)for K+ nitrogen, sulfur), and the polarity of the ~ o l v e n t . ~This - ~ crown over Na+, in H20, to the greater difference in hydration energies ethedion specificity can also serve as a simple model for of the cations relative to their intrinisic complexation energies understanding host-guest recognition in enzymes. Crown with 18c6. This observation underscores the importance of ethers are of particular interest to research efforts in environsolvent effects in crown ether chemistry. Indeed, it is possible mental remediation. For instance, 90Sr and 13’Cs are the two to change the selectivity of 18c6 for the alkali metal cations by major generators of heat in nuclear waste which complicates changing the dielectric constant of the solvent, offering schemes disposal. A more thorough understanding of cationkrown ether for selectively loading and unloading crown-cation c ~ m p l e x e s . ~ ~ ~ solution chemistry may provide the basis for rational design of Though their theoretical results indicated that 18c6 exhibited new ligands useful in separation of these and other radionuclides selectivity for Rb+ over K+ (contrary to experiment), Wipff et from complex waste streams and hazardous waste storage al. did demonstrate the high degree of conformational flexibility facilities. One example of the use of crown ethers for of 18c6 in interactions with the alkali metal cations. Later, radionuclide separation is the strontium extraction (SREX) Howard et a1.I9 extended this work to include atom-centered process, described by Horowitz, Dietz, and Fisher, which uses polarizable dipoles and three-body exchange repulsion terms di-tert-butylcyclohexano- 18-crown-6for recovering ?3r2+from in a MM potential. They described comformational energetics acidic solution.6 There is also a growing interest in the use of and complexation enthalpies of crown ethers with H30+, N b + , crown ethers, cryptands, and other ligands for radio-immunoand (CH3)3NH+ with satisfactory agreement to experiment. Dang therapy treatment of carcinomas.’ The immunoconjugate and Kollmanl* reported good agreement with experiment for consists of a ligand-radionuclide-cation (e.g., 224Ra2+)comthe free energy of K+/18c6 binding in H20 by calculating the plex that is covalently attached to a monoclonal antibody, potential of mean force (PMF). They observed a minimum in specific for tumor antigens. The immunoconjugate selectively the PMF at 1.75 8, K+-crown center of mass distance, binds to the tumor, delivering therapeutic doses of radionuclide significantly different from what is observed in the crystal to the tumor site, where the subsequent decay of the radionuclide structure, which they attributed to favorable interaction of K+ inactivates tumor cells. Fluoroionophores, consisting of a with solvent. Marrone et aLzoreported a PMF calculation of fluorophore linked to an ionophore (e.g., crown ether), represent 18c6 and K+ in methanol with an interesting modification of another interesting use for crown ethers and related macrocyclic the MM atom-centered charges. During the simulation, they ligands.* Measurable changes in the photophysical properties periodically perform a semiempirical AM124Hartree-Fock (HF) of the fluorophore upon ion binding by the ionophore may be calculation of 18c6 in the field of the solvent, represented as _____ point charges, and fit a new set of atom-centered charges to the Abstract published in Advance ACS Abstracts, September 1, 1994.
I. Introduction
@
0022-3654/94/2098-10465$04.50/0
0 1994 American Chemical Society
10466 J. Phys. Chem., Vol. 98, No. 41, 1994 electrostatic potential of the HF wave function. These new MM atomic charges are then used to continue the simulation. With this scheme, Marrone et al. are able to describe polarization of the crown ether as it binds K+. More recently, Glendening et al.” have utilized ab initio methods to compute binding enthalpies of 18c6 for the alkali metal cations Li+, Na+, K+, Rb+, and Cs+. Their study utilized large basis sets and results obtained at the second-order Moller-Plesset (MP2) perturbation level. Their calculations clearly show that solvation effects strongly influence cation selectivity. Their gas-phase calculations show 18c6 preferentially binds Li+, not Kf, as found in aqueous environments. However, they demonstrated that K+ selectivity is recovered when even a few waters of hydration are considered. The question of macrocycle conformation in solution and the role it may play in the complexation process has been the subject of several s t ~ d i e s . ~ ~X-ray - * ~ structures commonly show the uncomplexed 18c6 molecule to be in a Ci conformation, while K+/18c6 is predominately D3d.28-34 Theoretical studies have confirmed that the C, structure is lower in energy than D3d in the gas p h a ~ e . ’ ~ $ ’Sun ~ J ~and KollmanZ7 used free energy perturbation (FEP) calculations to obtain the relative free energies of various conformers of 18c6 in H20. They found that the D3d conformation of uncomplexed 18c6 was 5 kcaY mol more stable than the Ci conformation in H20, consistent with interpretation of experiment. Using FEP calculations, very long MD simulations, and symmetry arguments, Straatsma and M ~ C a m m o n calculate ~ ~ ~ * ~ the free energy of an ensemble of 18c6 molecules with D3d symmetry as 0.62 kcaYmol higher in energy than an unconstrained population in aqueous solution. They demonstrate that 18c6 is highly flexible, adopting many configurations in solution, with D 3 d as the predominant configuration. In a series of three papers, Ha and Chakraborty10s21s22 systematically describe the issue of crown flexibility and preorganization of 18c6 in polar media. They extensively parametrized a MM force field (including nonadditive terms) for the interaction of N h + with 18c6 using more than 150 data points obtained with density functional theory (DFT). To date, the MD and Monte Carlo studies of crown ethers have all employed MM-type potentials using both pair-wise additive as well as many-body effects.16J8,20-23~25.27 Recently, methods that hybridize quantum mechanics (QM) with MM have begun to be utilized with good success.35-40 These hybrid QMMM methods treat part of the system with QM (typically the solute(s) of interest or substrate and residues of an enzyme active site) and the remainder of the system with MM (such as the solvent). The two motifs are then coupled by a QM/MM interaction Hamiltonian. Hybrid QMMM methods have the potential advantage of being able to describe the role of the electrons (e.g., molecular polarization and charge transfer) as a natural part of the simulation, as well as treat chemical reactions. One of our goals is to suggest new or modified crown ethers and related ligands useful for separation of radionuclides from nuclear waste. We anticipate that as we test new ligand structures, a QM-based method will not require the extensive reparametrization of typical MM-based force fields. 18Crown-6 is considered a “hard donor”, i.e., less polarizable relative to “soft donors” such as aza and thia crown ethers.5J4341-43 Howard et al.19 found polarization effects to be significant in oxygen-containing crown ethers interacting with cations. Marrone et aL20 found that the AM l electrostatic-potential-fitted charges, for the oxygens of 18c6, changed significantly during the complexation of K+ by 18c6 in methanol, suggesting the importance of consideration for polarization effects. In this and future communications, we will describe the extent to which
Thompson et al. QM polarization is important in the theoretical description of the crown ethers in solution. We anticipate that such polarization effects will be even more important in the aza and thia crown ethers as well as crown ethers functionalized with charged or highly polarizable groups. Also, we will extend our use of the QMRMM methods to treat reactive systems, including macrocyclic enzyme mimics43 as well as enzyme active sites. Before we study reactive systems, it is important to benchmark these methods on large flexible molecules to determine whether they are capable of accurately calculating free energy quantities for crown ethers and related macrocycles in solution. Gains in computer hardware performance as well as newer more efficient software implementations of hybrid Hamiltonian methods (vide infra) make QM/MM simulations of condensed phase crown ether systems tractable. In this communication, we describe our initial implementation of a QM/MM method and benchmark an initial parametrization useful for study of K+/18c6 in aqueous systems using MD. Due to the size of the QM solute (18c6), we initially focus our efforts on the use of semiempirical Hamiltonians to represent the QM motif. Our treatment of hybrid Hamiltonian systems is easily extensible to include more accurate ab initio and DFT methods. In section 11, we present our method, review the QM/MM technique, and describe the implementation we have developed. In section HI,we present the parametrization scheme and results of MD simulations. In section IV, we present conclusions and discuss limitations and future use of the QM/MM method. 11. Methods
Hybrid QMMM methods have been described several times in the literat~re.~~-~O We simply review the basic formalism here for completeness. For a system composed of a QM solute in a MM solution, we can write an effective Hamiltonian
where HQMrepresents the quqtum mechanical region, Hm the molecular mechanical regjon, H Q the potential ~ that couples these two regions, and Hboundw the boundary region outside the system. In our study, we treat the boundary region by using periodic boundary conditions with the minimum image conventiomU The complete solute/solvent interaction is given by HQM/MM which we write as
(2) where indices i and a represent QM electrons and nuclei, respectively, and index s represents MM atoms. The first two terms on the right-hand side of eq 2 represent the interaction between the MM atoms and the QM electrons and nuclei, respectively, and the third term is the van der Waals interaction between QM nuclei and MM atoms. We use the combining rules for the Lennard-Jones parameters: E , = (Eacs)1’2 and u, = ‘/2(ua us). In our study, the ether molecule is treated by the Hartree-Fock molecular orbital theory, while the solvent and KC are described with MM. Specifically, we use the AM1 parametri~ation~~ of the semiempirical neglect of diatomic differential overlap (NDDO) model H a m i l t ~ n i a nto~ represent ~ the QM solute while H20 and K+ are treated with SPC/e threesite model for H20& and the parameters of respectively. We treat K+ non-quantum mechanically because there are no reliable semiempirical parameters for this cation. Wave function polarization effects in the QM motif arise from inclusion of the MM-atom/QM-electron term added to the one-electron
+
J. Phys. Chem., Vol. 98, No. 41, I994 10467
Nature of K+/Crown Ether Interactions matrix used in the Fock equations. We treat the AM1 implementation of the electrostatic terms in eq 2 in a fashion similar to Field et ~ 1 . ~ ~ Ab initio restricted Hartree-Fock (RHF) and second-order Moller-Plesset perturbation (MP2) calculations were performed with the GAUSSIAN 9248 and GAMESS49 programs. For calculations involving K+/DME/H20, the 6-3 1+G* basis was employed. For 18c6 calculations, a hybrid basis set was used. This basis set consisted of the standard 6-31+G* set for oxygen and hydrogen and the 6-31G* set for carbon. Carbon diffuse functions were judged to be relatively unimportant for describing ether-cation interactions and were neglected. The resulting basis set for 18c6 consisted of 342 contracted functions. For potassium, we used a (5s5p)/[3s2p] contraction of Hay and Wadt's valence basis sets with effective core potentials (ECP).50 The outermost (n - 1) shell of core electrons was treated explicitly, the field of the remainder being described by the ECP. Full geometry optimization of each structure was performed using the "tight" gradient convergence threshold of GAUSSIAN 92 or a threshold of 0.000 03 au for GAMESS. Calculations typically required more than a month of CPU time on an IBM Model 340. Gas-phase binding energies for the M+/18c6 and M+/18c6/H20 complexes were calculated as previously described." The full counterpoise (CP) correction of Boys and Bernardisl was applied to each ab initio binding energy reported here. Semiempirical and QM/MM calculations, MD simulations, and analyses of MD trajectories used the Argus computer program version 3.0.52 Classical MD simulations were carried out in the NVT ensemble in a periodic system using the minimum image convention.44 In all cases, the combination of periodic box size and nonbonded cutoffs were such that the QM solute did not interact directly with its image in neighboring boxes. Periodic box sizes were 20 A (K+/H20, DME/H20 simulations) and 25 A (K+/18c6/H20). SHAKE constraint^^^ were used for rigid SPC/e H20, DME (C-H bonds), and 18c6 (C-H, C-C, and C-0 bonds). For the QM molecules, constrained bond lengths were the AM1 optimized gas-phase values from the QM/MM K+/18c6 optimized geometry: C-H (1.121 A), C-C (1.524 A), and C-0 (1.426 A). Velocities were scaled by the method of Berendsen to simulate contact with an infinite thermal bath.54 Nonbonded interactions employed group cutoffs. Typical nonbonded cutoff distances were 9-10 A for MM/MM interactions and 10-12 A for QM/MM interactions. Time steps of 1.0 fs (K+/H20, DME/H20) and 1.5 fs (K+/18c6/H20) were used. The general protocol for preparing the system and collecting simulation data was (1) minimize the energy of the entire system with steepest descents for 10-50 steps, (2) anneal to 50 K for 1-2 ps using a temperature coupling constant of z = 0.020 ps, (3) equilibrate system to 300 K for 50- 100 ps, z = 0.150 ps, and (4) simulate at 300 K, z = 0.150 ps. Configurations were periodically saved to disk for later analysis. Simulations were carried out on the NERSC Cray C90, HP 735, and DEC AXP 3000/500 workstations. Initial configurations for solutes embedded in a box of H20 were obtained using InsightII 2.2 (Biosym Technologies). These files were then converted to input files for subsequent simulations with Argus. Our implementation of the QM/MM code52deserves further discussion. To date, most QM/MM implementations have been achieved by merging large mature codes, often with highly dissimilar data models, which can hinder the extensibility of the application. Typically, these implementations are driven by the simulation code, including necessary drivers from the QM codes to supply forces on QM atoms necessary for the MD.
Legend A: AM1 B: HF/6-31+G* hybrid C: Expt.
\ //T
7" 8 A: 111.7' B 113.6O
A:4316A B:4.792A c:4.267 A
k
89
I\
\
Dihedral 1-2-34 A: 186.2' B: 192.9' c:184.50
Figure 1. Ci conformation of 18c6. Selected geometric parameters are shown for (A) AM1, (B) HF/6-31+G*(hybrid basis)," and (C) experimental X-ray structure.28
Our approach was to develop a new data model that was, from its inception, designed for a multi-Hamiltonian description of systems ranging from single isolated molecules in the gas phase through condensed phase descriptions of solutes as large as proteins. Our data model and software were designed with techniques borrowed from object-oriented design principle^,^^ and the program was implemented in ANSI standard C.56 While our approach required us to design and write an entire electronic structure and empirical force field MD code, we find that the flexibility to readily extend the methods within a highly modular data model makes our approach more useful for our long-term goals.
III. Results and Discussion
A. Parametrization of the QM/MM Force Field. We began by assessing the quality of the AM1 method on energy minimized structures for dimethyl ether (DME) and 18-crown-6 (18c6). AM 1 results for DME have been previously reported57 and agree well with experiment for structure58 and dipole moment.59 Our implementation of the AM1 method has been extensively tested and in all cases reproduces the published results. AM1 gives the following results for DME (experimental values in parentheses): C-0-C angle 112.9' (1 11.3'), C-0 bond length 1.415 8,(1.410 A), and dipole moment 1.45 D (1.30 D) . The X-ray structure for uncomplexed 18c6 has been reported by Dunitz et al. to possess Ci symmetry.28 Beginning with the X-ray structure, we obtained optimized geometries at the AM1 and ab initio HF/6-31+G* level. The results are summarized in Figure 1 for selected geometrical elements of 18c6. Both AM1 and ab initio methods retained the Ci symmetry of the X-ray structure. Comparison of both ab initio and AM1 results with experiment is generally good. The AM1 bond lengths are generally longer than experiment. The C-C and C-0 differences are well within the average unsigned error in bond lengths reported in the literature for AM1 (0.050 while the error in C-H bond length is somewhat larger at 0.083 A. Both the COC and CCO bond angles show good agreement with experiment and are within the average unsigned error for AM1 bond angles of 3.3°.57 The CCOC dihedral angle shown is in excellent agreement with experiment. The transannular 0-0 distance of 4.316 A for AM1 is in reasonable agreement with
Thompson et al.
10468 J. Phys. Chem., Vol. 98, No. 41, 1994 Legend A: QM/MM (AMl$PC/e, Dang (K+) B: MW6-31&* hybrid C: Experiment
A 1.425 A
B: 1393 A B: 1 0 9 . 8 O
R Az2.735A
R A:2.5466j
B: 2.643 di
A: 5 74 A B: 5.802 A
B: 2 . 0 di
I
Dihedral 1-2-3-4 A: M.2" B: 75.40
Figure 2. D3d conformation of 18c6. Selected geometric parameters shown are for (A) AM1 and (B) HF/6-31+G*(hybrid basis)."
the experiment value of 4.267 A, considering the sensitivity of this geometry element to bond lengths and angles around the ring. Next, we compare the AM1 and ab initio calculated D3d structure of 18c6, which is summarized in Figure 2. The C-C and C-0 bond lengths generally agree with similar accuracy observed for Ci 18c6 (Figure 1). The angles COC and CCO are also in good agreement. The OCCO dihedral angle is calculated to be 88.2" for AM1 and 75.4' for the ab initio result. The AM1 dihedral angle lies between the ab initio result and a recently reported DFT structure by Ha and Chakraborty, who give a value of 97'.1° We have no ready explanation for the difference between the three results, especially in light of the good agreement between ab initio and AM1 values for the OCCO angle of gauche-1,Zdimethoxyethane (vide infra). The transannular 0-0 distance for AM1 is 5.974 A, which is longer than the ab initio value of 5.802 A, a result consistent with the larger OCCO dihedral angle for AM1. The 18c6 Ci structure is 4.4 kcaVmol more stable than D3d for the ab initio result at HF/6-31+G* (5.4 kcaYmol at the MP2 level using the HF geometry). Using AM1, we calculate the Ci structure to be more stable than D3d by 4.2 kcaVmo1, in good agreement with the ab initio HF result. This compares with Howard et al.,l9 who calculated an energy difference of 1.1 kcaY mol using a polarizable variant of the AMBER force field, and Sun and K ~ l l m a n who, , ~ ~ more recently, reported an energy difference of 2.0 kcaYmol using a nonpolarizable AMBER force field. Rotation barriers for AM 1 are generally underestimated relative to both experiment and ab initio results.60 The smallest "rotamer subunit" of 18c6 can be modeled with 1,2-dimethoxyethane in which there are three torsion angles involving heavy atoms: COCC, OCCO, and CCOC (the OCCO angle is the central angle in 1,Zdimethoxyethaneon which we focus). We can label each of these torsions with the labels s (syn), g (gauche), and a (anti). Using AM1,we calculate the syn-anti barrier (asa-aaa) to be 5.0 kcaVmol compared to a value of 8.9 kcaYmol reported by Jaffe et al. (obtained at the ab initio MP2/ D95+(2df,p) level).61 Our underestimation of this barrier is caused, in part, by the lack of proper description of repulsion between nonbonded centers inherent in all zero differential overlap (ZDC; methods. In the present study, we will attempt to assess the magnitude ci this effect in benchmark aqueous
Binding Energy A: -17.9 kcaVmol B: -18.8 kcaVmol C: -17.9 kcaVmol (AH)
Binding Energy A: -17.5 kcaVmol B: -18.9 kcaVmol
t Q
A:2.849A B: 2.845 A
Binding Energy A: -5.4 kcaVmol B: -5.4 kcaVmol
Figure 3. Model systems used to parametrize the Q M M M method: K+/H20, K+/DME, and DME/H20. Experimental binding enthalpy for K+/H20 is from Dzidic and Kebarle.*'
simulations of K+/18c6 that may suggest improvements in the model Hamiltonian methods. The SPC/e model has been shown to give good agreement with experiment for radial distribution functions, density, dielectric constant, and the diffusion constant of H ~ 0 . 4The ~ K+ parameters, from Dang>7 are designed to reproduce thermodynamic and structural properties of K+ in H20.62*63 To date, all reported uses of QM/MM methods have begun with the assumption that the QM and MM methods separately describe their respective motifs reasonably well and hence are not modified in the subsequent QM/MM parametri~ation.~~-~~ On the basis of the overall quality of the above results, and with the above mentioned caveat on rotation barriers, we proceed in a similar fashion. The solute-solvent interaction, described entirely by the H Q ~ term M (eq 2), contains all degrees of freedom which we use to fit parameters. The free parameters in HQM/MM that we vary are the LJ parameters E and 0 for the QM atoms (eq 2, term 3) and the AM1 parameters Qo and M a used in the two-electron coulomb integral and core-core repulsions, respectively, for the MM atoms (eq 2, terms 1 and 2). The fundamental interactions which we examine in the parametrization scheme are the gas-phase binding energies and minimized structures for K+/DME, DME/H20, and K+/H20 shown in Figure 3. In our approach, the K+/H20 interaction is already fixed by the MM parameters for these two We then utilize ab initio (MP2/6-3l+G*/CP corrected) calculations to parametrize HQW to the geometry and binding energy for K+/DME. Finally, we compare to the ab initio geometry and gas-phase binding energy for DME/H20 using the QM/ MM parameters fitted for K+/DME. All fits to the ab initio results were straightforward and done by hand. For the QM atoms, we began with the LJ parameters given by G ~ owhich , ~ ~ themselves are based on the OPLS parameters of Jorgensen et ~ 1 We . ~then explored various modifications to the LJ
J. Phys. Chem., Vol. 98, No. 41, 1994 10469
Nature of K+/Crown Ether Interactions parameters for DME to reproduce the ab initio K+/DME binding energy and geometry given in Figure 3. We observe that fitting the U parameters for the DME oxygen (DME(0)) is sufficient to achieve reasonable agreement with the ab initio result. A more systematic approach to this parametrization, including all the alkali metal cations, as well as other ether molecules, will be described in a forthcoming communication. The eo and aM values of 0.0 8, and 5.0 respectively, used by Field et ~ 1 gave the best fits to the ab initio data and were unmodified in our work. The K+/DME fitted parameters were then checked against the ab initio DME/H20 binding energy and geometry (Figure 3). Note that the use of the SPC/e water model (polarized for the bulk phase) will give a DME/H;?Obinding energy that is incorrect, relative to experiment, mainly by (1) the self-energy of polarizing the H20 molecule (i.e., the work required to generate polarized H20 from its gas-phase state) and (2) the incorrect polarization for H20 in the cluster as described by SPC/e. Berendsen et al. reported the self-polarization energy for the SPC/e model to be 1.25 kcal/mol.46 Thus, DME and H20, at infinite separation, will be too high in energy by this amount. For the DME/H;?Ocluster, natural energy decomposition analysis (NEDA)65results for our ab initio structure show an induced dipole moment of 0.18 D for H20, relative to infinite separation. This represents a “true” cluster polarization for H20 that is about half of the induced dipole for SPCle water, which is 0.5 D.46 Thus, to first order, DME/H20, with H20 having the “true” cluster polarization, should be lower in energy (more stable) than our DME/H20(SPC/e) result. Our binding energy for DME/H20 (Figure 3) should be decreased by the self-energy of SPC/e H20 but also increased by consideration of the “proper” cluster polarization for H20; these two sources of error tend to cancel. In addition, the error from the self-energy for SPC/e H20 should be larger than the error from H20 cluster polarization since the difference in the SPC/e H20 dipole relative to gas-phase H20 (0.5 D) is larger than the difference between SPC/e and our ab initio H20 induced dipole in the cluster (0.32 D). Thus, an upper limit on the error to our binding energy will be the 1.25 kcaymol self-polarization energy for the SPC/e H20. Hence, our true binding energy for DME/H20 will have a lower limit of 4.2 kcal/mol. We note the use of pairwiseadditive MM models in a method that explicitly includes polarization effects in only the QM solute can inhibit a consistent description of both structure and energetics (which is one possible source of error in all QM/MM studies to date). A better description of the cluster results should be achieved with a QM/ MM method that includes polarization in the MM motif. The K+ cation presents a large electric field perturbation on the DME molecule. To measure the response of the AM1 wave function to K+, we calculate the induced dipole moment (Ap) of DME to be 1.0 D in the QMlMM optimized K+/DME structure (Figure 3). For comparison, NEDA65of the ab initio (HF/6-31+G*) minimized structure for K+/DME gives Ap = 1.26 D for DME. While we observe the semiempirical method slightly underestimates total polarization relative to ab initio results (an observation consistent with the valence-only nature of the basis set expansions used in semiempirical Hamiltonians), we find the Ap results to be in generally good agreement. This result can be contrasted with the MM study of Wipff et a1.l’ They concluded that multiple charge schemes for the ether oxygens were necessary to reproduce experimental results for DME gas-phase dipole moment and binding enthalpies for K+-H20 and Na+-H20. For example, they found an oxygen charge of -0.3 gave a dipole moment (p) for DME of 1.3 D, in agreement with experiment and was the charge they used
TABLE 1: Core-Repulsion and Lennard-Jones Parameters Used in H Q ~ atom a.A-1 On. A 0.A E, kcaVmol C 0
H
.
~
Dimethyl Ether AM1 3.5w AM1 2.700 2.000 AM1
AM 1“ AM1 AM1
~
0.080 0.190 0.0 10
Water
5.w 5.0 5 .O
0 H K+
0.0 0.0 0.0
3.166
0.1554
0.0
0.0
3.332
0.100
Carbon and hydrogen values from G ~ o ; ~ ~ a Standard AM1 values.” oxygen values modified from G ~ o .Field ~ ~ et u L . ~ ~
-20
1
I
I
I
I
1
2
3
4
5
K+:Oxygen (DME)Distance (A)
Figure 4. Binding energy curve for K+/dimethyl ether (DME) comparing MP2/6-3 l+G*(CP) and QM/MM methods. Geometry for DME is from the MP2/6-31+G* optimized K+/DME structure. All other points obtained by varying R(K+-DME(0)).
for free, uncomplexed crown ethers. However, an oxygen charge of -0.6 (which gave a DME dipole moment of 2.48 D) was found to be optimal for cation-bound 18c6. We note that the QM/MM method describes both of these interactions well (DME p = 1.45D and 2.45 D for free and K+-complexed DME, respectively) within the same scheme and underscores the importance of including molecular polarization, especially for systems that include ion-neutral interactions. The final parameters we use for the HQW term are summarized in Table 1. The differences in K+/DME(O) distance is -0.05 8, relative to the ab initio result (Figure 3). This error is within the average unsigned error for bond lengths given with the AM1 method.57 The QM/MM K+/DME binding energy is less stable than the ab initio result by 1.4 kcaymol. We attribute the largest source of error in the binding energy to lack of (1) charge transfer effects, (2) polarization in the K+ cation, and (3) the quality of molecular polarization given by the minimum basis set expansions used in valence-only semiempirical methods. Indeed, NEDA65 of the ab initio K+/DME structure gives -2.2 kcdmol contribution to the binding energy by charge transfer effects. Nevertheless, we deem our parametrization scheme to be within acceptable ranges given the known limitations of the methods used and uncertainties in experiments. To further test the QM/MM parameters, we compare the QM/MM and ab initio binding energy curve for K+/DME in Figure 4. For both curves, the geometry for DME is obtained from the ab initio optimized geometry for K+/DME. We observe good agreement in most regions of the binding
Thompson et al.
10470 J. Phys. Chem., Vol. 98, No. 41, 1994
DME(O)M20(0) rdf (Figure 6) shows the first solvation shell centered -2.9 8, and integrates to 1.6 H20 0 : O pairs (Figure 6). Both of these results are consistent with -1.5 H20 molecules H-bonding with the oxygen of DME. The slight hump at -3.2 8, in Figure 5 represents the second hydrogen for the H-bonding water. This feature is mostly obscured by the larger second solvation shell centered at -4.5 8,. The DME(H)/H2O(O) rdf has essentially no structure imposed on the solvent at the H end of DME (data not shown). Our structural results are somewhat different from those reported for OPLS simulations by Jorgensen and co-workers6'@ which showed that ethers have 1-1.2 H20's H-bonding to the oxygen atom. To better understand the solute-solvent energetics, we describe the total interaction energy between a QM solute and the MM solvent in a fashion similar to Gao and Xia.36 The total solute/solvent interaction consists of electrostat$ and van der Waals interactions. The electrostatic piece of Heff (eq 1) can be written
2
1.5
2 z '
0.5
0
Figure 5. Dimethyl ether oxygen-water hydrogen radial distribution function from QM/MM MD simulation. Analysis from configurations saved during a 200 ps simulation (NVT), 300 K.
(3) where & is the zeroth-order gas-phase Hamiltonian ( ~ Q Mfrom eq 1) and is the electrostatic part of A Q ~ (the M first two terms on the right hand side of eq 2). The total system and the gas-phase wave functions, Y and YO,are defined by k$lY) = and fiolY0) = EolYo), respectively. We define the electrostatic solute- solvent interaction energy as
*dm
el") Ef;,
=
