Molecular Mechanical Force Field

J. Phys. Chem. 1995,99, 4794-4804

4794

Hybrid Quantum MechanicaVMolecular Mechanical Force Field Development for Large Flexible Molecules: A Molecular Dynamics Study of 18-Crown-6 Mark A. Thompson Environmental Molecular Sciences Laboratory, Pacijk Northwest Laboratory, Richland, Washington 99352 Received: October 20, 1994; In Final Form: January 3, 1995@

We present a hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics study of 18-crown-6 ( 1 8 ~ 6 in ) both polar hydrogen bonding and nonpolar solvents. The QM/MM method described here is based on our previous QM/MM study of K+/18c6l which employed the semiempirical AM1 method for 18c6 and the SPC/e model for H20. We improve our previous QM/MM method by including MM torsion terms to better describe the energetics of OCCO and COCC rotamers. The torsion terms are parametrized The resulting torsion-modified against the ab inirio calculations of Jaffe et al. for 1,2-dimetho~yethane.~,~ AM1 (TAM1) Hamiltonian is used with our previous QM/MM parameters to describe the conformational energetics of free 18c6 in both H2O and CCb. The TAMl method describes optimized geometries and relative energetics for Ciand D3d configurations of 18c6 slightly better than does the unmodified AM1 method. TAMl simulation results of OCCO and COCC torsion angle distributions show 18c6 maintains a preorganized oxygenlined cavity in H20,in contrast to the simulation in CCh. The crown ether strongly binds two bridging H20 molecules that help maintain the 18c6 cavity in H20. These two waters, together, have a simulation average 18c6/H20stabilization energy of -23.8 kcaymol. These results contrast markedly with an 18c6/H20 simulation using the unmodified AM1 method. The total QM/MM 18c6/H20 interaction energy is -82.6 kcallmol and -63.7 kcallmol for the TAMl and AM1 simulations, respectively. The AM1 18c6/H20 simulation reveals that 18c6 lacks both a preorganized cavity and the two bound waters. Without correction for rotamer energetics, the AM1 method is not useful for simulations of flexible systems like 18c6. However, semiempirical QM/MM methods, with corrections for torsions, are perhaps preferable to both pairwise-additive and polarizable-empirical force field methods, since the QM treats electronic polarization as a natural part of the simulation.

I. Introduction Synthetic macrocycles have drawn much experimental and theoretical interest since Pederson first synthesized the crown ether 18-crown-6 ( 1 8 ~ 6 )in 1967.4$5 Crown ethers show a remarkable range of specificity for a wide variety of cations that depends, in part, on the size of the ether, the type of donor atoms (e.g., oxygen, nitrogen, sulfur), and the polarity of the solvent.6-8 Crown ethers are of particular interest to research efforts in environmental remediation. For example, at the Hanford nuclear facility 90Sr2+ and 137Csf are two major generators of heat which complicate the disposal of nuclear waste. A more thorough understanding of fundamental interactions of cationlcrown ether solution chemistry may provide the basis for rational design of new ligands useful in the separation of these and other radionuclides from radionuclidecontaining waste streams at hazardous waste storage facilities. Crown ethers have been extensively studied from both a t h e o r e t i ~ a l l , ~and - ~ ~e ~ p e r i m e n t a l ~perspective. ~ ~ ~ ~ ~ ~ ~18-~~ Crown-6 is the prototypical hostlguest system that provides a tractable, yet challenging, model system to test new methods. With 18 unique torsion terms, each involving four heavy atoms, there can be a large number of rotational isomers (rotamers) for the uncomplexed 18c6 macrocycle in solution. Crown ethers are thought to be preorganized for complexation in polar solvents but not in apolar solvent^.^^-^^.^^ A correct theoretical treatment of free energies of hostlguest complexation will require, at the least, an adequate treatment of the conformational energetics of 18c6. Previously, we presented a study, of 18c6 (interacting @

Abstract published in Advance ACS Abstrucrs, March 1, 1995.

with K+ in H20) that utilized QM directly in a MD simulation.’ In this communication, we continue that work by addressing issues of conformational energetics for uncomplexed 18c6 as described with the hybrid quantum mechanicallmolecular mechanical method. The question of macrocycle conformation in solution and the role it may play in the hostlguest complexation process has been the subject of several ~t~die~.l2,14,15,17-19,21,22,26,27,47,48 X-ray structures commonly show the uncomplexed 18c6 molecule to be in a Ci conformation, while K+/18c6 complexes are predominately D3d.34936-41Theoretical studies have confirmed that the Ci structure is lower in energy than D3d in the gas phase. 10*13,28 In polar hydrogen bonding solvents, both experiment and theory indicate that 18c6 (and other derivatives) maintains ‘‘D~like’’conformations (Le., with a cavity preorganized for cation binding). In a Monte Carlo (MC) study of rigid Ci,D3d, and C1 conformations of 18c6 in a H20 cluster, Ranghino et al. showed that the D3d conformation is the most stable due to conformation-dependent solute-solvent interactions.12 Sun and Kollman used free energy perturbation (FEP) calculations to obtain the relative free energies of various conformations of 18c6 in H20.21 They found the D3d conformation of uncomplexed 18c6 was 5 kcallmol more stable than the Ci conformation in HzO, consistent with the interpretation of experiment. Using FEP calculations, very long MD simulations, and symmetry arguments, Straatsma and McCammon demonstrated that 18c6 is highly flexible, adopting many configurations in solution, with D3d-like as the predominant c ~ n f i g u r a t i o n . ’ ~In , ~ ~a series of three papers, Ha and C h a k r a b ~ r t y ~systematically ~ , ~ ~ . ~ ~ described the issue of crown flexibility and preorganization of 18c6 in polar media. Their

QO22-3654/95/2099-4794$09.00/Q 0 1995 American Chemical Society

Molecular Dynamics Study of 18-Crown-6 results also indicate that D3d-like configurations predominate in H20 and more Ci-like configurations predominate in CC4. In a recent MD study, Kowal and Geiger employ a modified MM2 force field to study 18c6 in H20.27 Their results showed two H2O molecules interacting strongly with the crown ether to maintain D ~ l i k econformations. In a 2 ns simulation of 18c6 in acetonitrile, Troxler and WipfP8 showed that 18c6 favors preorganized &-like conformations. Recently, Ha and Chakraborty have employed a reference interaction site model (XRISM) to study the solvation thermodynamics of 18c6 and dibenzo-18-crown-6 in H20, CC4, and a ~ e t o n i t r i l e . As ~ ~ with Troxler and Wipff, they also observe favorable solvation for the D3d conformation of 18c6 in acetonitrile. To date, MD and MC studies of crown ethers have all employed MM-type potentials, using both pairwise additive and many-body effect^.^^,^^,^^.^^,^^,^^,^^-^^,^^ In recent years, methods that hybridize quantum mechanics (QM) with MM have been utilized successfully.1~49-54These hybrid QMMM methods treat part of the system (typically a solute of interest) with QM, and the remainder of the system (typically solvent) with MM. The two motifs are then coupled by a QMMM interaction Hamiltonian. QM/MM methods have the potential advantage of being able to describe the role of the electrons as a natural part of the simulation. 18-Crown-6 is considered a “hard donor”, Le., less polarizable relative to “soft donors” such as aza and thia crown However, using a MM method that included atom-centeredpolarizable dipoles, Howard et al. found polarization effects to be significant in oxygencontaining crown-ethers interacting with cations.13 In another MM study, Marrone et al. found that AM1 electrostaticpotential-fitted charges, for the oxygens of 18c6, changed significantly during the complexation of K+ by 18c6 in methanol, suggesting the importance of polarization effects.26 We are interested in exploring the extent to which QM polarization is important in the theoretical description of crown ethers and related macrocycles in solution and developing QMbased simulation methods useful for large flexible systems. We anticipate polarization effects will be important in the aza and thia crown ethers, crown-ethers functionalized with charged or highly polarizable groups, and macrocycles interacting with alkaline earth dications. We are also interested in extending these methods to treat reactive macrocycles that can act as enzyme mimic^.^^.^^ Recently, we developed a QM/MM force field and utilized it in a MD study of K+/18c6 in H20.l We demonstrated that QM-based simulations for systems as large as 18c6 are feasible. We furthermore demonstrated a significant contribution from QM polarization in 18c6, in response to both the cation and solvent. However, we observed that the unmodified semiempirical AM1 method described the crown ether as too flexible, possibly lacking an adequate description of conformational energetics. While this problem was only minimal in the K+complexed 18c6 (due to the cation maintaining the macrocycle in &-like conformations), it could adversely affect the conformational statistics of the uncomplexed crown ether and prevent an adequate description of the free-energy of cation complexation. In this communication, we describe modifications to our previous QMMM Hamiltonian which improve the treatment of rotamer energetics in 18c6. We characterize this Hamiltonian with MD calculations of 18c6 in H 2 0 and CC4. In section 11, we present the method and describe the implementation we have developed. In section 111, we present the parametrization scheme and results of MD simulations. In section IV, we present conclusions.

J. Phys. Chem., Vol. 99, No. 13, 1995 4795

11. Methods Hybrid QM/MM methods have been described several times in the literature'^^^-^^*^^ (to which the reader is referred for a more exhaustive discussion). As in our previous study,’ the solute is treated by the Hartree-Fock (HF) molecular orbital theory while the solvent is 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 n ~ ~ to represent 18c6 while H20 is treated with SPC/e three-site model for H z O ? ~and C c 4 uses the united atom model of Rebertus et al. 6o We use this CC4 model primarily to facilitate comparison with the 18c6/CC14 MC study of Ha and Chakraborty.’* Wave-function polarization effects arise in the QM motif due to the addition of the MM-atom/QM-electron QM/MM interaction term added to the one-electron matrix used in the Fock equations.’ We treat the AM1 implementation of the QM/MM electrostatic terms in a fashion similar to Field et al.49 In this work, we augment the semiempirical HF Hamiltonian by including MM torsion terms for the heavy atom torsions about C - 0 and C-C bonds in 18c6 (vide infra). Semiempirical and QM/MM calculations, MD simulations, and analyses of MD trajectories used the Argus computer program version 3.0.1361-65 Classical MD simulations were carried out in the NVT ensemble in a periodic system using the minimum image convention.66 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 25 A and contained the 18c6 solute and 498 H20 molecules or 87 CC4 molecules. SHAKE constraint^^^ were used for rigid SPC/e H20 and 18c6 (C-H, C-C, and C - 0 bonds). For the QM solute, bond lengths were constrained to the AM1 optimized gas-phase values for 18c6: 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.68 Nonbonded interactions employed neutral group cutoffs. Both MMMM and QMMM nonbonded lists were constructed with 12.5 A cutoffs, while interactions calculated from these lists were truncated at 10 A. The time step was 1.5 fs The protocol for preparing the system and collecting simulation data was (1) minimize the energy of the entire system using steepest descents for 10-50 steps, (2) anneal to 50 K for 1-2 ps using a temperature coupling constant t = 0.020 ps, (3) equilibrate system in six equal phases of 6 ps each, with target temperatures from 50 to 300 K and corresponding t = 0.7 ps to 0.2 ps (velocities were rethermalized to new target temperature at the beginning of each phase), and (4) simulate at 300 K, t = 0.20 ps Configurations were saved to disk every 45 fs for later analysis. Simulations were performed on an IBM Model 590 workstation. The simulations typically accumulated -25 ps of simulation time/ day on an IBM 590 workstation. This performance makes potential of mean force free-energy calculations tractable for systems the size of 18c6, where total simulation times of 1-1.5 ns would require < 2 months of total computing time. Ongoing tuning of the code for performance and parallelization is in progress. 111. Results and Discussion

A. Parameterization of the QMlMM Force Field. To date, all reported uses of QM/MM methods have begun with the assumption that the QM and MM methods separately describe their respective systems reasonably well and hence are not modified in the subsequent Q m M parametrization. We will show that this approach is inadequate for flexible molecules 1349-54

4796 J. Phys. Chem., Vol. 99, No. 13, 1995

0

I

20

t

40

1

60

1

80

1

100

1

120

1

140

Thompson

1

160

I80

OCCO dihedral angle (deg)

0

COCC dihedral angle (deg)

Figure 1. Torsion potential energy for the central OCCO torsion of 1,2-dimethoxyethane(the two COCC torsions are maintained anti for all points in these curves). Energies calculated with (A) unmodified AMI, (B) AM1 with the OCCO angle augmented with a single MMtorsion term having K+ = 1.O kcdmol and functional form E$ = &[ 1 cos(3#)], and (C) AM1 augmented with a single MM-torsion term having Km = 1.4225 kcdmol.

+

Figure 2. Torsion potential energy for a terminal COCC torsion of 1,2-dimethoxyethane(the other COCC and central OCCO torsions are maintained anti for all points in these curves). Energies calculated with (A) unmodified AM1, (B) AM1 with the COCC angle augmented with a single MM-torsion term having K+ = 1.0 kcaYmol and functional form E+ = K+[1 + COS(^^)].

such as 18c6 when the current generation of semiempirical methods are used to treat the QM. Rotation barriers for AM1 are generally underestimated relative to both experiment and ab initio re~ults.6~ The smallest "rotamer-subunit" of 18c6 can be modeled with 1,Zdimethoxyethane (12DME) in which there are three torsion angles involving heavy atoms: COCC, OCCO, and CCOC. We can designate each of these torsions with the labels s (syn), g (gauche), and a (anti). For example, the all-anti conformation of 12DME is designated aaa. Using AM1, we calculate the central C-C syn-anti (asa-aaa) energy difference to be 5.0 kcaV mol, compared to a value of 8.9 kcdmol reported by Jaffe et al. (ab initio MP2/D95+(2df,~)).*,~Interestingly, the MM potential used by both Kowal and Geiger2' (modified MM2) and Straatsma and M ~ C a m m o n ' (GROMOS) ~*~~ give the asaaaa energy difference as -5 kcaYmo1, in good agreement with AM1, but in poor agreement with the ab initio results of Jaffee et al. AM1 underestimates this barrier, in part, because zero differential overlap (ZDO) methods lack a proper description of the repulsion between nonbonded centers. In Figures 1 and 2 we show calculated torsion potentials for the OCCO and COCC dihedral angles, respectively, in I2DME. The configurations for OCCO and COCC are generated by first minimizing the geometry for the asa and saa conformations, respectively; then increasing the selected torsion angle from 0" (syn) to 180" (anti) in increments of lo", while maintaining all other geometric elements fixed. The resulting single-point energies for these configurations describe the potential energy curves in Figures 1 and 2. Since all other degrees of freedom are frozen in these curves, they do not strictly represent minimum energy paths (which should differ by only -1 kcaY mol). For the purposes of our study, this should not be a serious issue, since we wish to demonstrate how addition of the MM torsion terms to the AM1 Hamiltonian can have dramatic effects on both the qualitative and quantitative rotamer statistics for 18c6 in solution. In Figure 1 we show the OCCO curve for AM1 and for two different potentials in which we augment the AM1 Hamiltonian with a single MM torsion term. The torsion term has the functional form: E# = K#[1 COS(~$J)]. For the

OCCO dihedral term, we show two MM-augmented curves which use Kb = 1.0 and 1.4225 kcdmol, respectively. Most evident in the AM1 curve is an asa-aaa energy difference of 4.3 kcaYmo1 and the complete lack of a banier between the aga and aaa structures. Addition of MM-torsion terms to AMI dramatically alters the description of this curve. With K4 = 1.4225 kcaVmol (curve C) we see the asa-aaa energy difference has increased to 7.2 kcaVmo1 and the aga-aaa barrier is now 2 kcal/mol. By comparison, Jaffe et al. report an asa-aaa energy difference for the OCCO angle of 8.90 kcaYmo1, and the agaaaa barrier as 2.31 kcaYm01.~.~ In Figure 2 we show a similar analysis applied to the COCC torsion angle (with the other two torsions in 12DME maintained as anti). The AM1 curve gives a gaa-aaa barrier of -0.25 k c d mol. Addition of a MM torsion term (same functional form as OCCO), with K# = 1.0 kcaYmo1, increases the gaa-aaa barrier to 1.25 kcaYmol and increases the saa-aaa energy difference from 15.5 kcdmol to 17.5 kcaYmo1. Jaffe et al. report a value of 2.04 kcaumol for the gaa-aaa To further test the MM-augmented AM1 Hamiltonian, we employ K# = 1.4225 kcaVmol and K4 = 1.O kcdmol for the OCCO and COCC torsion terms, respectively. We energy minimized the asa, aga, and aaa configurations for the central OCCO angle. The structural and energetic results are presented in Table 1. The OCCO angle in the aga-optimized structure is 70.8", which compares to 73.6" reported by Jaffe et al. We calculate the aga-aaa energy difference as 1.42 kcdmol compared to 0.15 kcaumol reported by Jaffe et al. The above observations have implications on the use of QM/ MM methods to treat large flexible systems that presently, and for the foreseeable future, remain intractable to ab initio and density functional theory (DFT) methods. In this regard, we stress the need for new model Hamiltonians that better describe intermolecular interactions and geometries away from minima. Such work is being pursued in our laboratory and elsewhere.70 To proceed with the present study, we chose to augment the AM1 method by including dihedral angle potential terms for all heavy atom torsions (COCC and OCCO) in 1 8 ~ 6 .One advantage of our approach over polarizable-MM methods is that

+

J. Phys. Chem., Vol. 99, No. 13, 1995 4797

Molecular Dynamics Study of 18-Crown-6 TABLE 1: Structural and Energetic Results for Minimized Conformations of 12DME' Using the TAMlb Hamiltonian OCCO angle (deg) COCC angle (deg) COC (deg) OCC (deg)

180.0 180.0 112.2 105.9 1.5251 1.4266

Relative Energy (kcaYmo1)

0.0

cc (4 oc (interior) (A)

70.8 179.6 112.2 108.3 1S233 1.4250 1.42

Legend A: AM1 B: TAMl C: HF/6-31+G* hybrid

A: 1.523 A R: 1.523 A A: 111.3"

c: 1.511 A

0.0

180.0 112.2 108.8 1S268 1.4217 7.74

1,ZDimethoxyethane. The MM torsion-augmented AM 1 form of the QM/MM Hamiltonian. & = 1.4225 and 1.O kcdmol for the OCCO and COCC torsion angles, respectively.

A: 1.425 A R: 1.424 A C: 1.393 A

a

Legend A: AM1 B: TAMl C: HF/6-31+G* hybrid D: Experiment

B: 5.745 A

11

Dihedral OCCO A: 88.2" R: 75.9' c: 75.40

I

\\

//

COCC A: 171.7" R: 178.6" C: 176.7'

Q Figure 4. D3d conformation of 18c6. Selected geometric parameters shown are for (A) AM1, (B) TAMl (see text), and ( C ) HF/6-31+G*(hybrid basis).28

Dihedral 1-2-3-4 A: 186.2' R: 184.8' C: 192.9' D: 184.5"

Figure 3. Ciconformation of 18c6. Selected geometric elements are shown for: (A) AM1, (B) TAMl (see text), ( C ) HF/6-31+G*(hybrid basis),28and (D) experimental X-ray structure.34

QM-based methods should treat electronic effects more realistically. We view the semiempirical QMMM methods (even MM-augmented QMMM methods) as a possible alternative to polarizable-MM methods in this regard. For convenience only, we hereafter refer to the MM-torsion-augmented AM1 method as TAM1 (torsion-augmented AM1). We chose not to pursue an exhaustive parameter search for the torsion parameters nor utilize more sophisticated functional forms of torsion potentials but rather to employ values used in the above calculations (Figures 1 and 2), that generally give good agreement with the results of Jaffe et al. and show clear improvement over the unmodified AM1 results. Thus, we use K4 = 1.4225 and 1.0 kcal/mol for the OCCO and COCC (CCOC) torsions, respectively. All other QM/MM parameters are the same as those previously described. The X-ray structure for uncomplexed 18c6 has been reported by Dunitz et al. to exhibit Ci symmetry.34 Beginning with this X-ray structure, we obtained optimized geometries with TAM1. The results are presented in Figure 3 along with previously reported AM1 and ab initio HF/6-31+G* and the experimental X-ray results of Dunitz et al.34 Previously, we reported generally good agreement of AM1 with the experimental Ci structure.' Here, we observe a slight improvement in the description of the transannular 0. -0distance from 4.3 16 A for AM1 to 4.237 A for TAM1, which is in better agreement with the experimental value of 4.267 A.34 Also, there is a small improvement in the COCC dihedral angle from 186.2" for AM1

to 184.8" for TAM1, the latter of which is also in better agreement with the experimental value of 184.5°.34 In Figure 4 we show the optimized structure of the uncomplexed D 3 d configuration of 18c6 and compare it with our previously reported AM1 and ab initio HF/6-3 1+G* results. Generally, the results of TAM1 are similar to AM1, and both are in good agreement with the ab initio results. However, the OCCO dihedral angle does show significant improvement from an AM1 value of 88.2" to 75.9" for TAM1, which compares to 75.4" for the ab initio result. The transannular 0 -0distance is also improved slightly from 5.974 8, for AM1 to 5.745 A for TAM1, compared to the ab initio value of 5.802 A, a result consistent with the smaller OCCO dihedral angle for TAM1. The 18c6 Ci structure is 4.4 kcaVmol more stable than D3d for the ab initio result at HF/6-31+G* (5.4 kcaVmo1 at the MP2 level using the HF Using AM1, we previously calculated the Ci structure to be more stable than D3d by 4.2 kcaVmol in good agreement with the ab initio HF result.' The TAM1 result shows the Ci structure to be more stable than D3d by 6.6 kcal/mol, in good agreement with our previous MP2 This compares with Howard et al.,13 who calculated an energy difference of 1.1 kcaVmo1 using a polarizable variant of the AMBER force field. Sun and Kollman2*report an energy difference of 2.0 kcal/mol using a nonpolarizableAMBER force field. Using a modified MM2 potential, Kowal and Geiger recently reported D3d to be nearly isoenergetic with the Ci configuration, with Ci more stable by only 0.3 kcaVm01.~~ B. Molecular Dynamics Results for 18c6 in H20. We performed an MD simulation of 18c6 in a 25 x 25 x 25 A3 box of H20, using periodic boundary conditions (section 11). The QM/MM Hamiltonian used the TAMl QM Hamiltonian for 18c6 and the SPC/e model for H20. The periodic box contained 498 H20 molecules. The initial configuration chosen for 18c6 was the TAM1 optimized D3d structure (Figure 4). The simulationphase was run for 300 ps, using 1.5 fs time steps, with positions and velocities saved every 45 fs All analyses were performed on the saved configurations. The simulation averaged distribution of OCCO and COCC dihedral angles are shown in Figures 5 and 6, respectively. The OCCO distribution has a distinct peak centered near 65" with a smaller peak near 180". Integration of this distribution shows

Thompson

4798 J. Phys. Chem., Vol. 99, No. 13, 1995 0.2

0.1

0.1s 0.15

$

i g

Bp

i

I

I

0.1

B

0.1

0.05

0.05

5 1 : OCCO dlhedral angle (deg)

6

20

40

60

80

100

120

140

160

180

OCCO dihedral angle (deg)

Figure 5. Unsigned distribution of the OCCO dihedral angle value along with the running total averaged for all six OCCO angles in 18~6, of OCCO angles under the curve. Results shown are for the TAMl Hamiltonian. The analysis is from configurations saved during a 300 ps simulation of 18c6 in HzO.

Figure 7. Unsigned distribution of the OCCO dihedral angle value averaged for all six OCCO angles in 18c6, along with the running total of OCCO angles under the curve. Results shown are for the unmodified AM1 Hamiltonian. The analysis is from configurations saved during a 200 ps simulation of 18c6 in HzO.

0.6-

0.5

0.5

-I

g

M

0

+

II

-

3a

0.4-

8

0.3-

I

0.4-

8

U 0.3U

0.2-

8

o.lj

0.24 I

0.1

1

0 0- n

20

40

60

80

1 100

120

140

160

180 0

COCC dihedral angle (de@

Figure 6. Unsigned distribution of the COCC dihedral angle value averaged for all 12 COCC angles in 18c6, along with the running total of COCC angles under the curve. Results shown are for the TAMl Hamiltonian. The analysis is from configurations saved during a 300 ps simulation of 18c6 in HzO.

that, on average, five OCCO torsions are accounted for in the 65" peak and one in the 180" peak. The COCC distribution shows, on average, 9 of 12 COCC torsions centered about 180°, with three COCC torsions having a value centered near 70". For reference, the 18c6 D3d optimized geometry has all OCCO and COCC angles at 75.9" and 178.6", respectively (Figure 4). Visualization of the simulation shows a distinct oxygen lined cavity being maintained by 1 8 ~ 6with , one oxygen occasionally being pointed up and slightly away from the cavity (not shown). For comparison, we ran a 200 ps simulation of 18c6 under the same conditions as above, but using the unmodified AM1 method in the QM/MM Hamiltonian. We show the resulting OCCO and COCC dihedral angle distributions in Figures 7 and 8, respectively. The difference between the TAMl and AM1

20

40

60

80

100

120

140

160

180

COCC dihedral angle (deg)

Figure 8. Unsigned distribution of the COCC dihedral angle value averaged for all 12 COCC angles in 18c6, along with the running total of COCC angles under the curve. Results shown are for the unmodified AM1 Hamiltonian. The analysis is from configurations saved during a 200 ps simulation of 18c6 in H20.

results is marked. The latter shows a much lower, broad OCCO distribution, mainly centered around 85", but with significant occupations in all values up to 180". The AM1 COCC result is also decidedly non-&-like, with significant occupancy in angles about 80", and a much reduced peak near 180", relative to the TAMl results. Visualization of the AM1 simulation results shows a highly flexible 18c6 macrocycle having a largely disordered structure (not shown), with no discernible oxygenlined cavity as observed in the TAM1 result. Our TAMl result is very similar to the observations of Ha and Chakraborty,I8 who observed similar OCCO and COCC torsion distributions for 18c6 in H20. The TAMl result also agrees with the experimental interpretation that crown ethers are preorganized for cation binding in polar protic solvent^.^^^^^ We next demonstrate how solvent effects help to maintain the 18c6 cavity in the TAMl result. In Figure 9, we show the

Molecular Dynamics Study of 18-Crown-6

R (A)

Figure 9. Radial distribution function for the separation of the 18c6 center-of-mass (COM) and the H20 (COM) along with the running coordination number. Results shown are for the TAMl Hamiltonian. The analysis is from configurations saved during a 300 ps simulation.

radial distribution function (rdf) for the 18c6 center-of-mass (COM) and H?;O(COM)from the TAMl simulation. There is a distinct peak, centered near 1.5-1.75 A, which integrates to two waters (out to a 3.0 8, separation). This peak is the first solvation shell for the 18c6 cavity; in essence, the bound waters that act as guest molecules in the absence of metal cations. The sharpness of the peak and the distinct minimum between the first and second peaks in the rdf indicate that these waters are tightly bound to the 18c6 cavity. The smaller peak near 4.04.5 8, is the second solvation shell and integrates to -6 waters. The peak near 7.0 8, is solvent in van der Waals contact with the periphery of the 18c6 macrocycle. Our results are very similar to those of Kowal and Geiger, who also observe two waters strongly bound in the cavity of 1 8 ~ 6 . ~Their ' 18c6/H20 rdfs look very similar to ours. In addition to the sharp peak centered near 1.5-1.75 8,, they observe the second solvation shell near 4.0 A, as well as the peak near 7.0 8,. Our simulation results show one water interacting on either side of the crown macrocycle plane. In Figure 10 we show a typical configuration from the TAMl simulation showing the two bound waters. Ranghino et al. observed conformation-dependent solvation in a MC study of 18c6 in a cluster of 100 waters.12 They utilized a rigid crownether, constrained in Ci,C1, and D3d conformations. For the perfect D3d configuration, they observed a total of four waters having significant interactions with the crown cavity in which (on each side) one water H-bonds to two crown oxygens and the second water H-bonds with a third crown oxygen and the oxygen of the first water. We note that Ranghino's result is from a rigidly constrained D3d crown, without thermal fluctuations in the macrocycle (which are included in both our and Kowal's simulations). For comparison, we show in Figure 11, the 18c6(COM)/H20(COM)rdf using the unmodified AM1 QM/ MM Hamiltonian. The difference between the TAMl and AM1 method is again marked. Both the first and second solvation peaks are missing in the AM1 result. Integrating the AM1 rdf out to 3.0 8, gives -0.4 waters. This observation is consistent with the lack of a well-defined cavity for the AM1 result. To better understand the average orientation of the tightly bound waters, we measured the orientation of the H20 dipole

J. Phys. Chem., Vol. 99, No. 13, 1995 4799

vector with a vector connecting the H20 oxygen and the COM of 18c6. The H20 dipole vector is defined as the vector originating on the H20 oxygen and bisecting the H-0-H angle. For reference, if an H20 dipole is oriented directly at the 18c6 COM for D3d crown ether, the measured angle is exactly zero. In Figure 12, we show the distribution of this angle for the strongly bound waters (Le., those in the first peak of the g(R) in Figure 9). The distribution in Figure 12 is centered near 40°, which indicates that the two bound waters generally have their hydrogens oriented toward the 18c6 cavity, in position to maximize H-bonding with the 18c6 oxygens. For comparison, Kowal and Geiger observe this s q e distribution to be strongly oriented near 0" (Le., both hydrogens pointed almost straight at the COM of the crown, on average).27 It seems apparent from their study that they also observe the crown maintaining virtually pure &-like configurations. Our two results are in qualitative agreement. We attribute the differences, in part, to the slightly different 18c6 conformational statistics (our crown shows occasional movement of an oxygen up and away from the cavity, which would displace the 18c6 COM slightly). The 18c6(0):H20(H) and 18c6(0):H20(0) rdfs from the TAMl simulation are shown in Figures 13 and 14, respectively. These distributions are averaged over all six oxygens of 18c6 and thus will reflect the correlated effects of the bound waters interacting with multiple crown oxygens. The first peak in the 18c6(0):H20(H) rdf is centered near 2.0 8, and integrates to two O:H pairs. The first peak for the 18c6(0):H20(0) is centered near 2.85 8, and integrates to two 0 : O pairs. These distributions show increased H20 occupation in the first solvation peaks, relative to our previous QM/h4M simulation of free dimethyl ether in H20,' a result consistent with the two additional strongly bound waters in 18c6. Our results differ somewhat from the MC simulations of Ha and Chakrab~rty,~' who report no discernible structure in H20 near the crown oxygens (their Figure 4), despite their observation of &-like conformations in H20. In a more recent study, Ha and Chakraborty report XRISM calculations of the solvation thermodynamics of 18c6 in in which they report solvent structuring about the crown-ether oxygens very similar to our results. To better understand the role of the bound waters, we report the results of a QM/MM interaction energy analysis. The details of this type of analysis have been previously described;'.50we proceed here in an identical fashion as our previous work.' To summarize, the total QM-solute/MM-solvent interaction energy consists of electrostatic and van der Waals interactions, which we write

The terms on the right-hand-side of eq 1 are: E l , the interaction of the solute's gas-phase wave function with the solvent (unpolarized solute-solvent interaction energy); the energy penalty for reorganizing the solute's gas-phase wave function to the solvated phase (solute cost to solvate); Epl/smb, the interaction of the distorted wavefunction with the solvent (polarization stabilization); and the van der Waals interaction energy between the QM and MM atoms. The total QM polarization response is given by (Epol)= (Edist Epo~smb). All of these interaction energy terms can be directly computed without resorting to perturbation theory.' The results we report are simulations averages over the saved configurations. The results of the above energy analysis for 18c6 in H20 are given in Table 2 for both the TAMl and AM1 simulations. The total interaction energy (Eint)between 18c6 and H20 is -82.6

+

Thompson

4800 J. Phys. Chem., Vol. 99, No. 13, 1995

Q

d

Figure 10. Stereoview of a sample configuration showing 18c6 and two bound waters. Results shown are for the TAMl Hamiltonian from a 300 ps simulation of 18c6 in H20.

0

i

i

i

i

io

0

i2

R (A)

Angle (de@

Figure 11. Radial distribution function for the separation of the 18c6 center-of-mass (COM) and the H20 (COM) and the running coordination number. Results shown are for the unmodified AM1 Hamiltonian. The analysis is from configurations saved during a 200 ps simulation.

Figure 12. Distribution of the angle between the H20 dipole vector and the vector connecting the H20 oxygen with the 18c6 (COM) (see text). Results shown are for the TAM1 Hamiltonian from a 300 ps simulation of 18c6 in H2O.

kcaVmol and -63.7 kcdmol for the TAMl and AM1 simulations, respectively. The TAMl -AM1 difference in Eht is -18.9 kcaumol (TAM1 more stable than AMl), and is roughly accounted for by: AE1 = -14.1 kcaVmo1; AEpol= -2.6 k c d : E A = -2.1 kcaVmo1 (difference of columns 1 mol; and & and 2 in Table 2). In Table 3, we show an interaction energy analysis applied to only those waters whose COM is within 3.0 from the COM of 18c6 (Le., the first peak in the 18c6(COM):H20(COM) rdf, Figures 9 and 11). For the TAMl simulation, the first solvation shell waters have an average Eint of -23.8 kcdmol with the crown ether. In Figure 15 we show the distribution of Eint for the bound waters. There are three distinct regions in the distribution, centered near -12, -25, and -35 kcdmol, respectively. Futher examination of our simulation data reveals that these regions represent configurations with, respectively, 1, 2, and 3 waters bound to the crown cavity (data not shown), with the majority of configurations having two bound waters. The simulation average of Eint per water in the first solvation shell is -12.3 kcaVmol (not shown). By comparison, Kowal and Geiger report the total interaction of the two bound waters (“type 1” in their work) with 18c6 as -9.93 kcaVmol per bound H20.27 Using the AM1 QM/MM Hamiltonian, we previously obtained an energy minimized structure for dimethyl ether

(DME)/H20 showing a single hydrogen bond to the DME oxygen with a binding energy of -5.4 kcaVmo1’ (this result is identical for both TAMl and AM1). Thus, the Eint value of -12.3 kcdmoVH20 is consistent with -2 hydrogen bonds between each bound H20 and 18c6. Recall, this same solvation peak is completely missing in the AM1 18c6(COM):H;?O(COM) rdf (Figure 11). Not surprisingly, the TAM1-AM1 difference in Eint for all of the solvent waters (Table 2, - 18.9 kcdmol) is fairly close to Eint for the first solvation peak in the TAM1 results from Table 3 (-23.8 kcaumol). Thus, we observe that the two bound waters account for nearly the entire difference between the TAMl and the unmodified AM1 simulations both in solute/solvent interaction energy as well as the difference in the conformation distribution of the macrocycle. The TAMl results are consistent with 18c6 being preorganized for cation binding, with two bound guest molecules (H20) that need to be desolvated during cation complexation. The simulation average Mulliken charge on the oxygens of 18c6 is -0.37 for TAM1, compared to the gas-phase optimized D3d geometry which gives a value of -0.28. In agreement with the polarization energy results (Table 2), we see significant QM repolarization upon solvating 18c6. The more polarizable the crown ether or the stronger the cation (e.g., Sr2+), the larger will be the change in polarization between the complexed and

Molecular Dynamics Study of 18-Crown-6

J. Phys. Chem., Vol. 99, No. 13, 1995 4801 TABLE 3: Interaction Energies for 18-Crown-6/H 0 from MD Simulation for Waters with COM within 3.0 of 18-Crown-6 COM"

2

I

energies (kcaymol) EinP

TAMlb -23.8

E i

-20.4

Edist Epovstab

Eple

Evdw QW

no. of H20f

1.6 -3.3 -1.7 -1.8 1.94

AM1' -3.3 -2.7 0.2 0.5 -0.2 -0.4 0.45

"Results are simulation averages from an analysis applied to configurations saved every 45 fs during the simulation. The MM torsion-augmented AM1 form of the QM/MM Hamiltonian (see text). 'The unmodified AM1-based QM/MM Hamiltonian used in our previous study.' Interaction energy terms defined in eq 1 (see text). e (E,!) = (&is[ -k EPustab). f Simulation average number of waters with COM within 3.0 %, of 18-crown-6 COM. R (A)

Figure 13. Radial distribution function for the separation of 18c6 oxygen and HzO hydrogens along with running coordination number. Results shown are for the TAMl Hamiltonian from a 300 ps simulation of 18c6 in H20.

ti Eint (keaUmol)

0

2

4

6

8

10

I2

R (A)

Figure 14. Radial distribution function for the separation of 18c6 oxygen and H20 oxygen along with running coordination number. Results shown are for the TAMl Hamiltonian from a 300 ps simulation of 18c6 in H20. TABLE 2: Interaction Energies for 18-Crown-6/H20 from MD Simulation" energies (kcaymol) TAMlb AMI' Eintd -82.6 -63.7 E1 -50.8 -36.7 Edist 8.4 5.8 .&ustab -16.9 -11.6 Epof -8.5 -5.9 Evdw -23.3 -21.2 QW

"Results are simulation averages from an analysis applied to configurations saved every 45 fs during the simulation. bThe MM torsion-augmented AM1 form of the QM/MM Hamiltonian (see text). 'The unmodified AM1-based QM/MM Hamiltonian used in our previous study.' Interaction energy terms defined in eq. 1 (see Text). e

= (Edist

+ Epoustab).

uncomplexed crown ether. Painvise additive MM models can only describe this effect in an average fashion, while polarizable MM models can, at best, mimic QM-based methods where polarization is a natural part of the calculation. C. Molecular Dynamics Results for 18c6 in CC4. We performed an MD simulation of 18c6 in a 25 x 25 x 25 A3

Figure 15. Distribution of total QM/MM interaction energies (Elnt; see text, eq 1) for only the 18c6-bound waters (first peak in Figure 9). Results shown are for the TAMl Hamiltonian from a 300 ps simulation of 18c6 in H20.

box of CC4, using periodic boundary conditions (section 11). The QM/MM Hamiltonian used the TAMl QM Hamiltonian for 18c6. The initial configuration chosen for 18c6 was the TAMl optimized D3d structure (Figure 4). The simulationphase was run for 200 ps, using 1.5 fs time steps, with positions and velocities saved every 45 fs All analyses were performed on the saved configurations. The simulation averaged distribution of OCCO and COCC dihedral angles are shown in Figures 16 and 17, respectively. The OCCO distribution is markedly different from that obtained in the H20 simulation (Figure 5). In the CC4 result, there is significant occupation for OCCO angles in the vicinity of 180", with smaller probability near 70". Integration of this curve, shows, on average, four OCCO torsions in the 180" and two in the 70" regime. The 18c6 COCC distribution in CCl4 (Figure 17) is also distinctly different from that obtained in H20 (Figure 6). Here, .we observe a dramatically reduced occupation in the near-180" regime, and significant increase in the 70" regime, relative to H20. Integration of the curve in Figure 17 shows, on average, eight COCC torsions in the 70" regime and four in the 180" regime. Taken together, these results describe 18c6 in a nonpolar solvent as not maintaining a preorganized cavity as in H20. Indeed, visualization of the simulation results shows the crown has some of its -CH2- groups turned inward,

Thompson

4802 J. Phys. Chem., Vol. 99, No. 13, 1995

i OCCO dihedral angle (deg)

Figure 16. Unsigned distribution of the OCCO dihedral angle, averaged over all six OCCO angles in 18c6, along with the running total of OCCO angles under the curve. Results shown are for the TAMl Hamiltonian. The analysis is from configurations saved during a 200 ps simulation of 18c6 in CC4 12

0.7

0.6

IV. Conclusions

10

0.5

8

0.4

8 U

U U

6

8

1P

0.3 4

:

0.2

2

0.1

0 ;

I

I

I

are being explored. Our results do not show distinctly Ci-like statistics, but we also observe many conformations being explored. We note that crown-ether conformational transitions involve concerted motions of several atoms. In this regard, it may be difficult for MC simulations to adequately explore phase space for these motions (unless this is specifically built into the MC sampling algorithm). Molecular dynamics simulations can describe concerted motions important to such conformational changes. Finally, it is possible that Ha's MC and our MD results are not exploring the same configurationalphase space for these finite-length simulations. Indeed Straatsma and McCammon's MD study of 18c6 in H20 indicate that even after 1.5 ns, the crown is still exploring new configuration^.'^ A more detailed comparison of the convergence properties of MC and MD simulations, for a flexible molecule such as 18c6, would be interesting. We performed an interaction energy analysis of 18c6 in CC4 (eq 1) Since the united-atom CC4 model is uncharged, all terms vdw will be exactly zero (i.e., Eht = EQmM). We except E:& calculate the simulation average interaction energy (EhJ for 18c6 with CC4 to be -28.1 kcal/mol, much less than the H20 result. The relative values of Eintfor 18c6 in H20 (-82.6 kcallmol) and in CC4 (-28.1 kcdmol) clearly demonstrate the larger degree of 18c6/solvent interaction in polar protic solvents relative to nonpolar solvents.

0

Figure 17. Unsigned distribution of the COCC dihedral angle, averaged for all 12 COCC angles in 18c6 along with the running total of COCC angles under the curve. Results shown are for the TAMl Hamiltonian. The analysis is from configurations saved during a 200 ps simulation of 18c6 in CCL.

essentially filling the cavity of 18c6, with several of the oxygens pointing up and/or away from the center of the macrocycle (not shown). The enhanced occupation near 180" in the OCCO distribution of Figure 16 is consistent with the electronegative oxygens being as far apart as possible in a nonpolar medium. Ha and Chakraborty, employing the same solvent model for CC4, also observed the crown-ether to be non-&~like,relative to the H20 results.'* However, their torsion distributions are somewhat different than ours. Ha's 18c6/CC4 COCC distribution is more broadly distributed over the range 60-180°, without distinct peaks as seen in our data (Figure 17). Their OCCO distribution does show distinct populations near 70" and 180°, as ours, but their distributions are more evenly weighted. They conclude that their C c 4 results indicate that Ci-like configurations predominate, though they observe that many conformations

We have presented molecular dynamics simulation results for 18c6 in H20 and in CCk, to better understand the role of solvent in the conformational energetics of 18c6 described by the QM/MM method. We employed a MM torsion-augmented AM1 Hamiltonian to describe the QM solute ( 1 8 ~ 6 )to obtain agreement for rotamer energetics of 1,2-dimethoxyethane with the ab initio results of Jaffe et al.*s3 Simulations using this torsion-augmented AM1 (TAM1) Hamiltonian show that the crown ether maintains a preorganized cavity and strongly binds two H20 molecules, one on either side of the macrocycle plane. These two H20 molecules interact with multiple oxygens from the crown and together have an average interaction energy with 18c6 of -23.8 kcaymol. This result is consistent with both waters maintaining -2 H-bonds to the crown ether oxygens. The total interaction energy of the crown with all the HzO solvent is -82.6 kcallmol. This is very different from simulations using unmodified AM1. In the latter case, the crown does not maintain a preorganized cavity and is completely missing the two bound waters. What is interesting about this result is that both TAMl and AM1 treat QM electrostatics and vdW terms identically, the only difference being the torsion augmentation in TAM1. Apparently, only a small barrier ( - 2 kcaU mol in free 1,2-dimethoxyethane)separating the gauche and anti conformations of the OCCO angles is sufficient to allow the crown oxygens to interact favorably with H20 relative to competition from thermal fluxuations that would drive the predominately gauche OCCO angles to their trans conformation, thereby breaking the multiple H-bonding pattern seen in the bound waters. The total solute/solvent interaction energy for the unmodified AM1 result is -63.7 kcallmol. The difference between the TAMl and AM1 18c6/H@ total interaction energies is close to the total interaction energy of TAMl 18c6 with just the two strongly bound waters. The TAMl results are consistent with both e ~ p e r i m e n t ~and ~ , ~previous * theoretical studies12.14,1*,21,27 which indicate that 18c6 is preorganized for cation binding in polar protic solvents. Total quantum mechanical polarization contributes -8.5 kcaY mol stabilization to the TAMl 18c6/H20 simulation, which

Molecular Dynamics Study of 18-Crown-6 represents -10% of the total solute/solvent interaction energy. The simulation average Mulliken charges on the 18c6 oxygens is -0.37 relative to a value of -0.28 for the gas-phase optimized D3d structure; again consistent with the polarization of 18c6 in H20. A TAMl simulation of 18c6 in CCb shows the crown distinctly lacking a preformed oxygen-lined cavity. Rather, we observe the crown with some of its -CH2- groups pointing inward and occupying the cavity. The OCCO torsion distribution for 18c6 shows -4 of these torsions maintaining angles near 180°, which is consistent with keeping the electronegative oxygens as far apart as possible in a nonpolar solvent. Our results for 18c6/HzO are generally consistent with the observations of Kowall and Geiger, who also report two strongly bound waters.27 Furthermore, our result for the 18c6 conformational differences between H20 and CCl4 simulations is qualitatively consistent with those previously reported by Ha and Chakraborty.’* What is interesting, and satisfying, about this comparison is that our potential is very different in origin and in its treatment of electrostatics than the MM potentials of both Kowall and Ha. Proper treatment of flexible systems has not been extensively discussed in the Q W M literature. It is self-evident that correct treatment of rotamer energetics is required to properly describe free energy quantities for flexible macrocycles with multiple low-lying conformations. In this regard, we conclude that the unmodified AM1 method is not generally useful for these types of QMhlM calculations. However, we note that the TAMl method does show promise of being able to calculate free energies of crown ether systems. Furthermore, we postulate that methods like TAMl may even be preferable to treating the solute with polarizable-MM methods, since QM-based methods should treat electronic polarization more realistically. This latter characteristic may be important to a correct and consistent treatment of more polarizable macrocycles or for treating strong cations (e.g., Sr2+). The TAMl simulation of 18c6/H20 clearly shows that the crown is preorganized for cation complexation and is strongly solvated. While preorganization favors binding of cations, strongly bound waters may compete with the cation. Both of these characteristics will be important to describing free energies of cation complexation of 18c6 in H20, which we are currently pursuing. Acknowledgment. This work was performed under the auspices of the Division of Chemical Sciences, Office of Basic Energy Research, U.S. Department of Energy, under Contract DE-AC06-76lUO 1830 with Battelle Memorial Institute, which operates the Pacific Northwest Laboratory, a multiprogram national laboratory. We gratefully acknowledge Eric Glendening for a critical reading of the manuscript. We would also like to acknowledge Thom Dunning and Bruce Garrett for encouragement. References and Notes (1) Thompson, M. A,; Glendening, E. D.; Feller, D. F. J. Phys. Chem. 1994, 98, 10465. (2) Smith, G. D.; Jaffe, R. L.; Yoon, D. Y. J. Phys. Chem. 1993, 97, 12752. (3) Jaffe, R. L.; Smith, G. D.; Yoon, D. Y. J. Phys. Chem. 1993, 97, 12745. (4) Pedersen, C. J. J. Am. Chem. SOC. 1967, 89, 7017. ( 5 ) Pedersen, C. J. J. Am. Chem. SOC. 1967, 89, 2495. (6) Parker, D. Adv. Inorg. Chem. Radiochem. 1983, 27, 1. (7) Samec, Z.; Papoff, P. Anal. Chem. 1990, 62, 1010. (8) De Jong, F.; Reinhoudt, D. N . Adv. Phys. Org. Chem. 1980, 17, 279.

J. Phys. Chem., Vol. 99, No. 13, 1995 4803 (9) Yamabe, T.; Hori, K.; Akagi, K.; Fukui, K. Tetrahedron 1979, 35, 1065. (10) Wipff, G.; Weiner, P.; Kollman, P. J. Am. Chem. SOC.1982, 104, 3249. (11) Hon, K.; Yamada, H.; Yamabe, T. Tetrahedron 1983, 39, 67. (12) Ranghino, G.; Romano, S.; Lehn, J. M.; Wipff, G. J. Am. Chem. SOC.1985, 107, 7873. (13) Howard, A. E.; Singh, U. C.; Billeter, M.; Kollman, P. A. J. Am. Chem. SOC. 1988, I1 0, 6984. (14) Straatsma, T. P.; McCammon, J. A. J. Phys. Chem. 1989,91, 3631. (15) Mazor, M. H.; McCammon, J. A,; Lybrand, T. P. J. Am. Chem. SOC.1990, 112, 4411. (16) Hancock, R. D. Acc. Chem. Res. 1990, 23, 253. (17) Dang, L. X.; Kollman, P. A. J. Am. Chem. SOC. 1990, 112, 5716. (18) Ha, Y. L.; Chakraborty, A. K. J. Phys. Chem. 1991, 95, 10781. (19) Ha, Y. L.; Chakraborty, A. K. J. Phys. Chem. 1992, 96, 6410. (20) Seidl, E. T.; Schaeffer ID,H. F. J. Phys. Chem. 1991, 95, 3589. (21) Sun, Y.; Kollman, P. A. J. Chem. Phys. 1992, 97, 5108. (22) Ha, Y. L.; Chakraborty, A. K. J. Phys. Chem. 1993, 97, 11291. (23) Rencsok, R.; Kaplan, T. A.; Harrison, J. F. J. Phys. Chem. 1993, 92, 9758. (24) Hay, B. P.; Rustad, J. R.; Hostetler. C. J. J. Am. Chem. SOC.1993, 115, 11158. (25) Leuwerink. F. T. H.: Harkema. S.: Briels. W. J.: Feil. D. J. Comout. Chem.’ 1993, 14, 899. (26) Marrone. T. J.: Hartsough, D. S.; M e n Jr., K. M. J. Phys. Chem. 1994,98, 1341. (27) Kowall, T.; Geiger, A. J. Phys. Chem. 1994, 98, 6216. (28) Glendening, E. E.; Feller, D.iThompson, M. A. J. Am. Chem. SOC. 1994, 116, 10657. (29) Ha, Y. L.; Chakraborty, A. K. J. Phys. Chem. 1994, 98, 11193. (30) Live, D.; Chan, S. I. J. Am. Chem. SOC.1975, 98, 3769. (31) Lehn, J.-M. Angew. Chem., Inr. Ed. Engl. 1988, 27, 89. (32) Cram, D. J. Angew. Chem., Int. Ed. Engl. 1988, 27, 1009. (33) Bourson, J.; Pouget, J.; Valeur, B. J. Phys. Chem. 1993, 97,4552. (34) Dunitz, J. D.; Seiler, P. Acta. Crystallogr. 1974, 830, 2739. (35) Horwitz, E. P.; Dietz, M. L.; Fisher, D. E. Solvent Extract. Ion Exchange 1991, 9, 1. (36) Seiler, P.; Dobler, M.; Dunitz, J. D. Acra Crystallogr. 1974, 830, 2744. (37) Dunitz, J. D.; Seiler, P. Acta Crystallogr. 1974, 830, 2750. (38) Dunitz, J. D.; Dobler, M.; Seiler, P.; Phizackerley, R. P. Acta. Crystallogr. 1974, B30, 2733. (39) Dobler, M.; Dunitz, J. D.; Seiler, P. Acta Crystallogr. 1974, B30, 2741. (40) Dobler, M.; Phizackerley, R. P. Acta Crystallogr. 1974, B30,2746. (41) Dobler, M.; Phizackerley, R. P. Acta Crysrallogr. 1974, B30,2748. (42) Doxsee. K. M.: Wieman. H. R.: Weaklev. T. J. J. Am. Chem. Soc. i99’2,1i4, (43) Srivanavit. C.: Zink. J. I.; Kechter. J. J. J. Am. Chem. SOC. 1977, 99,‘5$76. (44) Izatt, R. M.; Bradshaw, J. S.; Nielsen, S. A,; Lamb, J. D.; Christensen, J. J. Chem. Rev. 1985, 85, 271. (45) Gokel, G. Crown Ethers and Cryprands; The Royal Society of Chemistry: Cambridge, 1991; Vol. 3. (46) Li, G.; Still, W. C. J. Am. Chem. SOC.1993, 115, 3804. (47) Straatsma, T. P.; McCammon, J. A. J. Chem. Phys. 1989,90,3300. (48) Troxler, L.; Wipff, G. J. Am. Chem. SOC.1994, 116, 1468. (49) Field, M. J.; Bash, P. A,; Karplus, M. J. Comput. Chem. 1990, 11, 700. (50) Gao, J.; Xia, X. Science 1992, 258, 631. (51) Gao, J. J. Phys. Chem. 1992, 96, 6432. (52) Gao, J.; Chou, L. W.; Auerbach, A. Biophys. J. 1993, 65, 43. (53) Singh, U. C.; Kollman, P. A. J. Comput. Chem. 1986, 7, 718. (54) Staiton, R. V.; Hartsough, D. S.; Merz Jr, K. M. J. Phys. Chem. 1993, 97, 11868. ( 5 5 ) Coooer. S. R. Acc. Chem. Res. 1988. 21, 141. (56) Warihel, A.; Levitt, M. J. Mol. Biol. 1976, 103, 227. (57) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. J. Am. Chem. SOC.1985, 107, 3902. (58) Pople, J. A.; Santry, D. P.; Segal, G. A. J. Chem. Phys. 1965,43, S129. (59) Berendsen, H. J. C.; Grigera, J. R.; Straatsma, T. P. J. Phys. Chem. 1987, 91, 6269. (60) Rebertus, D. W.; Beme, B. J.; Chandler, D. J. Chem. Phys. 1979, 70, 3395. (61) Argus 3.0; Thompson, M. A.; 1994, Pacific Northwest Laboratory, Richland, WA 99352, [email protected] (62) Thompson, M. A,; Zemer, M. C.; Fajer, J. J. Phys. Chem. 1991, 95, 5693. (63) Thompson, M. A,; Zemer, M. C. J. Am. Chem. SOC.1991, 113, 8210. (64) Thompson, M. A.; Fajer, J. J. Phys. Chem. 1992, 96, 2933. (65) Barkigia, K. M.; Miura, M.; Thompson, M. A,; Fajer, J. Inorg. Chem. 1991, 30, 2233. I

Thompson

4804 J. Phys. Chem., Vol. 99, No. 13, 1995 (66) Allen, M. P.; Tildesley, D. J. Computer Simulations of Liquids; Oxford University Press: New York, 1987. (67) Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327. (68) Berendsen, H. I. C.; Postma, J. P. M.; van Gunsteren,W. F.; DiNola, A,; Haak, J. R. J. Chem. Phys. 1984, 81, 3684.

(69) Dixon, D. A.; Matsuzawa, N.; Walker, S . C. J. Phys. Chem. 1992, 96, 10740.

(70) Thiel, W.,personal communication. (71) Ha, Y. L.; Chakraborty, A. K. J. Phys. Chem. 1991, 95, 10781. JP942853H