Aggregation and Solvation of Steroid Molecules in Different Solvents

Two novel self-assemblies based on π-π and belly-to-belly interactions. Martha C. Mayorquín-Torres , Rafael Arcos-Ramos , Marcos Flores-Álamo , Ma...
0 downloads 0 Views 73KB Size
Download PDF
Aggregation and Solvation of Steroid Molecules in Different Solvents Juha T.

Pa¨iva¨rinta,*,†

Antti T.

Poso,‡

Matti

Hotokka,†

and Esa

Muttonen§

Department of Physical Chemistry, Åbo Akademi University, Porthansgatan 3-5, FIN-20500 Åbo, Finland, Department of Pharmaceutical Chemistry, University of Kuopio, P.O. Box 1627, 70211 Kuopio, Finland, and Orion Co., P.O. Box 65, FIN-02101 ESPOO, Finland

CRYSTAL GROWTH & DESIGN 2002 VOL. 2, NO. 2 121-126

Received November 16, 2001

ABSTRACT: Ab initio quantum chemical studies at the Hartree-Fock (HF) level with the 6-31G** basis set and molecular dynamics calculations were performed on two pharmaceutical compounds, namely, the budesonide and beclomethasone dipropionate molecules. The possible water binding positions in the budesonide molecule were analyzed using the PM3 Hamiltonian, which gives fairly good geometries for budesonide-water complexes but underestimates dramatically the stabilization energies of the complexes as compared to Hartree-Fock calculations. The most stable complex according to our Hartree-Fock calculations is one where water is bound to the keto group. Its stabilization energy is -30.2 kJ/mol. In the molecular dynamics simulations, clear differences in aggregation of steroid molecules were observed in different solvents. The conformations of steroid molecules from molecular dynamics simulations differ from structures predicted from quantum chemical calculations showing the importance of evaluation of structure in different solvents. Introduction A detailed knowledge of the microscopic properties of drug substances is indispensable for a deep understanding of many chemical and biological processes. The properties of drug molecules are influenced by solvent molecules. Thus, it is vital to investigate the behavior of molecules in different solvents. Crystal structures are determined by intermolecular interactions, and so the study and the understanding of interactions are of considerable importance in the development of crystal engineering. Weak interctions can exert an influence during crystallization that is out of all proportion to their final contribution to overall crystal stabilization energy. The roles played by some of the weaker interactions during crystallization could be of a specific nature.1 The adoption of the particular crystal structure depends on a number of factors, such as enthalpic, entropic, and kinetic. In the budesonide crystal where ordinary hydrogen bonded network has proved to be important for crystallization, molecular dynamics is useful tool for understanding the formation of dimers or bigger agglomerates and for studying the usual conformational changes in different solvents. The antiinflammatory corticosteroid budesonide, 11β,21-dihydroxy-16R,17R-propylmethylenedioxy-1,4-pregnadiene-3,20-dione, has been investigated with X-ray crystallography,2 molecular mechanics calculations,3 NMR analysis,4 mass spectrometry,5 and chromatography methods.6 The molecular mechanics calculations referred to the gas phase. The other steroid structure * To whom correspondence should be addressed. Current address: Åbo Akademi University, Department of Physical Chemistry, Porthansgatan 3-5, 20500 Åbo, Finland. E-mail: [email protected]. Phone: +358-2-215 4951. Fax: +358-2-215 4706. † Åbo Akademi University. ‡ University of Kuopio. § Orion Co.

studied is a salt of beclomethasone, beclomethasone dipropionate, 9-chloro-11β,17,21-trihydroxy-16β-methylpregna-1,4-diene-3,20-dione 17,21-dipropionate. Most chemical processes, in particular biochemical reactions, and industrial processes take place in condensed media.7 The influence of solvent on different molecules has been calculated by other authors using for instance different self-consistent reaction field (SCRF) models,8-10 the conductor-like screening model (COSMO) approach,11 and the explicit solvent method.12,13 The easiest way to handle the solvent environment is to use electrostatic screening, but this does not take into account for instance hydrophobic interactions. In the earliest classical continuum models,14,15 the solute and the solvent are treated by separated nonoverlapping wave functions. Studies of solvent effects involving the interaction of the solute with a proton donating solvent usually lead to hydrogen bonds. Therefore, it is very difficult to justify that the overlap of the electron densities between the two subsystems is neglected.16 Such models also neglect the possible exchange effects between the solute and the solvent. It has also been suggested that even the newest continuum models cannot treat hydrogen bonds, especially solvation energy.17-19 These problems can be overcome by explicitly including solvents. In this so-called “supermolecule approach”,20 the interaction between species A and B is calculated as the difference between the energies of the complex AB and the isolated species. The present study concentrates on interactions between budesonide and different solvents, namely, water, methanol, and ethanol, and between bechlometasone dipropionate and water, methanol, and acetone. Quantum chemical calculations are used in identifying the most stable complexes of both epimers (22R and 22S) of budesonide and water molecule using the “supermolecule” approach. We have compared these results with molecular dynamics simulations using the GROMACS

10.1021/cg010032d CCC: $22.00 © 2002 American Chemical Society Published on Web 02/06/2002

122

Crystal Growth & Design, Vol. 2, No. 2, 2002 Scheme 1

united atom force field21 and explicit water molecules. Molecular dynamics was also used to simulate aggregation of both steroid structures, budesonide and beclomethasone dipropionate, in different solvents. It has been shown that the C-H groups are able to form hydrogen bonds in crystal structure,22 but they are not considered here as groups capable of forming hydrogen bonds for the reasons given by Mills et al.23 The C-H hydrogen bonds may occur in crystals when hydrogen-bond donor groups are not present. They are so weak that their geometries are often distorted as compared to ordinary hydrogen bonds. Quantum Chemical Calculations. The ab initio calculations were done at the restricted Hartree-Fock level of approximation using the 6-31G** basis set24-27 and the Gaussian-94 program.27 This basis set has proven to be suitable for the complex formation incorporating hydrogen bond formation.28,29 The geometry of budesonide molecule was fully optimized using the PM3 Hamiltonian, starting from its X-ray structure.2 In the ensuing calculations on the budesonide-water complexes, this structure was then kept fixed. A grid was generated around the budesonide molecule. One water molecule was systemically placed at each grid point where its position was optimized. In this way, all the hydrophilic sites of budesonide were detected. All the detected budesonide-water complexes were afterward optimized at restricted Hartree-Fock level of approximation by using standard split-valence 6-31G** basis set allowing all geometrical parameters to change. Polarization functions were added to the basis set because it has been demonstrated that these are needed in hydrogen bond calculations.30 In the ab initio calculations, only the essential portions of the budesonide molecule were included to make the calculations feasible. This is not expected to affect the stabilization energies of water-budesonide complexes. In the optimization of the water-O26, waterO27, and water-O28-O31 complexes rings A and B, B and C, D and E with side chains, respectively, were included. The atom numbering used in the budesonide molecule in this paper is shown in Scheme 1. We assumed that stabilization energies for different epimers differ only in positions O28 and O29. Molecular Dynamics Simulations. The same structure of budesonide-S epimer was used in the beginning

Pa¨iva¨rinta et al.

of the molecular dynamics simulation calculations as in the PM3 calculations. Simulations were carried out for the following systems: (i) 10 budesonide/beclomethasone molecules with 1018 water molecules in an initial box size of 3.34 × 3.34 × 3.34 nm3; (ii) 10 budesonide/beclomethasone molecules with 1018 methanol molecules in an initial box size of 4.44 × 4.44 × 4.44 nm3; (iii) 10 budesonide/beclomethasone molecules with 1018 ethanol molecules in an initial box size of 5.13 × 5.13 × 5.13 nm3. For all simulations, the initial positions of budesonide molecules were first randomly placed, and after that the systems were solvated. All solvent molecules that approached budesonide molecules closer than the sum of their van der Waals radii were excluded. At the beginning of the simulation, the system was minimized 2000 cycles with steepest descent algorithm to remove all close contacts. The same principle was also employed in the simulations where a single budesonide (iv) and beclomethasone (v) molecule was simulated in different solvents. Molecular dynamics simulations were run using GROMACS 2.0.21 In the GROMACS force field, a united atom approach is used so that the hydrogen atoms of, e.g., a CHn group, are not considered explicitly. The only exception in our system is the A-ring where hydrogens are needed for a better description of the double-bonded carbon atoms. The Lennard-Jones 12-6 type and Coulomb terms describe the interactions between sites. The ring structure and the correct conformation of the epimer were kept at the desired orientation applying the improper torsion potentials that prevent transitions by the other possible conformations. In the hydroxyl groups and in the A-ring, hydrogens were maintained explicitly using charges recommended in GROMACS force field. Several of the parameters needed in our structure were not available in the force field. The OS atom type was used for oxygens 28 and 29. A proper description of electrostatic interaction in the force field is essential for valid calculations. To get a consistent description of the system, we have used the 6-31G* basis set in the assigning charges for oxygen atoms 28 and 29. The same basis set has been used in the development of the GROMACS force field. The validity of used intramolecular parameters was tested from the crystal structure.2 In the following text and tables, the CHn atom type means the carbon atom bonded to n hydrogen atom(s). In cases where a bare carbon (atom type CB or C) was bonded to a sp3 hybridized carbon atom CH1, parameters were adopted from the case when CH2 was used instead. All missing parameters were created using the same principle. The important missing parameters are listed in Table 1. The partial charge of the chloride atom, 0.1 e, was calculated with the same method as was used in hydrogen bond evaluation. The force constants for Cl bonding were adopted from the GROMOS-96 force field31 implemented in GROMACS, using modifications as instructed in GROMACS32 to transform GROMOS type potentials to harmonic ones. The geometry of acetone molecules was maintained fixed during MD simulations and parameters were adopted from the work by Jorgensen et al.33 and tested for instance using Monte Carlo simulations by Freitas et al.34 and molec-

Aggregation and Solvation of Steroid Molecules

Crystal Growth & Design, Vol. 2, No. 2, 2002 123

Table 1. Missing Parameters in GROMACS Force Field and Used Values for Reference Values and Force Constants GROMACS parametera CH1-CB CB-OS CB-CB CH2-CH2-CS1 CH2-CS1-OS OS-CS1-OS CB-CH1-CH2 OS-CB-C CH3-CB-CH2

used parameter CH2-CB C-OS CB-CBb CH2-CH2-CH1 CH2-CH1-OA

ref value (nm, angle) and force constant

example in structure

0.15300, 334720. 0.14350, 251040. 0.15300, 334720. 111.000, 460.240 111.000, 460.240 109.5000 111.000, 460.240 109.500, 460.240 111.000, 460.200

C16-C17 C17-O29 C13-C17 C22-C23-C24 C23-C22-O28 O28-C22-O29 C17-C16-C15 O29-C17-O20 C18-C13-C12

a Atom types explained in original paper and reference manual for GROMACS force field. b Reference value changed to describe better the bond found in budesonide.

ular dynamics by Ferrario et al.35 in work where acetone-water mixtures were investigated. A twin range cut off 0.9 nm was used for van der Waals interactions, and the particle-mesh Ewald summation method36,37 for electrostatic interactions, where the cut off was also 0.9 nm. The time step was 2 fs, using LINCS38 to constrain bonds with hydrogen atoms. We used NPT conditions in the simulations. A constant pressure, 1 bar, independently in all three directions was used with a coupling constant 1.0 ps. Budesonide and solvent were coupled separately to a temperature bath at 300 K, using a coupling constant 0.1 ps. Water was modeled as simple point charges (SPC) because it has been shown that the SPC model produces a better thermodynamic potential than flexible water models, say SPC/E39 and TIP3P.40 After energy minimization, the system was allowed to equilibrate. Analysis from simulations were carried out using 1 ns simulation time on 30 processors at the Center of Scientific Calculations, Finland. Analysis was performed using the facilities within GROMACS.

Figure 1. Optimized bond lengths (pm) and angles (degrees) for budesonide-water complex in the O26 position by the (a) AM1; (b) PM3; and (c) HF methods.

Results and Discussion Possible water binding positions and complex energies have been considered using the semiempirical PM3 and AM1 Hamiltonians and ab initio calculations at the Hartree-Fock 6-31G** level of approximation. The structures optimized using the PM3 and ab initio methods differ at the C22-C25 side chain. The changes in these optimized structures are not expected to affect the stabilization energies, because the same computational method was used throughout. The optimized water positions calculated with different methods are shown in Figures 1-6. The AM1 optimized structures of budesonide-water complexes differ dramatically from those obtained using the PM3 Hamiltonian and Hartree-Fock approximation. The AM1 optimizations lead to unreasonable bifurcated structures, which the PM3 and ab initio calculations do not predict, as is seen in Figures 1 and 4-6. The PM3 and ab initio structures resemble each other. The AM1 method predicts the longest and the PM3 method predicts the shortest hydrogen bond distances. In the PM3 optimized budesonide-water complex (1b), water tends to move toward a hydrogen atom of the A-ring. The distance between the nearest hydrogen atom in the A-ring and the oxygen atom of water is 233, 280, or 264 pm calculated with AM1 Hamiltonian, PM3 Hamiltonian, and RHF method, respectively. The hydrogen bonds in the epimers 22S (2) and 22 R (3) do not differ markedly. The one in (2) is slightly

Figure 2. Optimized relevant bond distances (pm) and angles (degrees) in the budesonide (22S)-water complex using (a) the PM3 method and (b) Hartree-Fock.

Figure 3. Optimized relevant bond distances (pm) and angles (degrees) in the budesonide (22R)-water complex using (a) the PM3 method and (b) Hartree-Fock approximation.

weaker than that in (3) because the side chain of budesonide distorts the structure, so that the angle O28-Hwater-Owater, 170 degrees, differs slightly from the optimal hydrogen bond angle. The AM1 optimized complexes (4a) and (5a) are similar to the previous ones. The hydrogen bonds in the PM3 and HF optimized structures (4b) and (4c) also are not typical because of steric effects caused by the Me18 group and the hydrogen atom at C16 atom. In (5),

124

Crystal Growth & Design, Vol. 2, No. 2, 2002

Pa¨iva¨rinta et al. Table 2. Calculated Stabilization Energies (kJ/mol) of Budesonide-Water Complexes Using the AM1 Hamiltonian, the PM3 Hamiltonian, and Hartree-Fock Approximation ∆E (AM1)

∆E (PM3)

∆E (HF)

-22.5 -22.5 -17.0 -18.9 -18.6 -17.7 -20.7 -21.8 -24.7

-16.1 -10.8 -9.2 -11.9 -11.1 -11.9 -12.3 -11.7 -20.6

-30.2 -21.9 -17.9 -20.4 -19.3 -19.3 -26.1 -27.7 -28.8

O26 O27 O28 (S) O28 (R) O29 (S) O29 (R) O30 O31 O30-O31

Table 3. Mulliken Atomic Charges Computed for the Isolated Budesonide and Budesonide-Water Complex (1) Where Relevant Charges Are Shown

Figure 4. Optimized relevant bond distances (pm) and angles (degrees) in the budesonide-water complex using the (a) AM1; (b) PM3; and (c) HF methods.

C1 C2 C3 C4 C5 H33a H32 O26 O (water) H (water) H (water)b

HF

PM3

AM1

HF

PM3

AM1

-0.06 -0.23 0.52 -0.25 0.08 0.20 0.16 -0.62 -0.71 0.32 0.37

-0.18 -0.30 0.39 -0.32 -0.04 0.24 0.19 -0.35 -0.44 0.18 0.23

-0.17 -0.29 0.32 -0.29 -0.04 0.25 0.19 -0.35 -0.47 0.24 0.24

-0.06 -0.21 0.50 -0.25 0.08 0.17 0.15 -0.58

-0.19 -0.29 0.37 -0.32 -0.06 0.22 0.19 -0.32

-0.17 -0.28 0.31 -0.29 -0.05 0.22 0.18 -0.31

a Weak interaction with water oxygen. b Forms a hydrogen bond in the PM3 and HF calculations.

Figure 5. Optimized relevant bond distances (pm) and angles (degrees) in the budesonide-water complex using the (a) AM1; (b) PM3; and (c) HF methods.

Figure 6. Optimized bond lengths (pm) and angles (degrees) for the budesonide-water complex in the O30-31 position by the (a) AM1; (b) PM3; and (c) HF methods.

water forms a clear hydrogen bond when calculated with PM3 Hamiltonian (5b), but in HF (5c) calculations the hydrogen bond angle is distorted. The most stable complex according to HF optimization, 30.2 kJ/mol, was found when the water forms complex in A-ring (1c), see Table 2. According to the AM1 and PM3 calculations, the most stable complex was found adjacent to the keto and hydroxyl groups at C17, where the stabilization energies are 24.7 and 20.6 kJ/ mol for (6a) and (6b), respectively. In HF calculations, the stabilization energy for (6c) is 28.8 kJ/mol. In the AM1 optimized complex, both hydrogens atoms in water molecule bind to O30. This is not the case in the PM3 Hamiltonian optimized complex.

Table 3 shows the Mulliken charges of budesonidewater complexes in O26 and isolated budesonide provided by different methods. Comparison between the PM3 and HF charges obtained for an isolated budesonide and a budesonide-water complex shows the enhancement of the negative charges of the oxygen atoms and the increase of the positive charge of the hydrogens which are making the hydrogen bond. In AM1 calculations, the charges of hydrogens in the water were equal. This was noticed in all complexes. This analysis shows that in the first solvation layer solvent molecules have influence on the electronic structure of budesonide. The solvent molecules are well organized around the hydroxyl groups according to the RDF analysis when one budesonide molecule is simulated in different solvents (Figures 7-8). The first two solvation shells are clearly seen in the RDF analysis. The Ohydroxyl-Osol distances in the first solvation shell are slightly shorter than 300 pm, which is the typical hydrogen bond distance. In the text and figures, we are using notations Ow, Ometh, and Oeth for oxygen atoms in the water, methanol, and ethanol solvents, respectively. The only differences between the solvation of hydroxyl groups are found in intensities. The different sizes of the simulated box causes changes for intensities. The solvation shells are remarkably better organized around the hydroxyl groups than around ketogroups, for example at O26 (Figure 9). A similar trend is seen at O30 (result not shown). The A- and C-rings are quite similar in budesonide and beclomethasone, so the solvation around polar groups is expected to be the same and omitted for further analysis. In the MD simulations, the RDF analysis shows that O28 and O29 atoms are not good acceptors for HB. This

Aggregation and Solvation of Steroid Molecules

Crystal Growth & Design, Vol. 2, No. 2, 2002 125

between budesonide and the solvent. The number of hydrogen bonds in the E-rings at each oxygen atom varies between 0 and 1 and in other acceptors between 1 and 2 in methanol and ethanol solvents. In the water environment, the number is usually increased by one, and also acceptor sites in E-rings are at all times forming hydrogen bonds. When 10 budesonide molecules were simulated, the number of HBs between budesonide molecules were 2-4, 2-6, and 2-4 in water, methanol, and ethanol solvents, respectively. The HBs were constructed mainly between hydroxyl groups and O36. The radius of gyration describes the compactness of the steroid aggregates

Figure 7. RDF of solvents around a hydroxyl group (O31) in a single budesonide simulation (iv).

Figure 8. RDF of solvents around a hydroxyl group (O27) in a single budesonide simulation (iv).

Figure 9. RDF of solvents around a keto group (O26) in a single budesonide simulation (iv).

is probable due to movements of the side chain of the budesonide molecule. Also, the partial charges in the oxygen atoms are less favorable for HB. In the crystal structure, the O28 atom is mentioned to be an HB acceptor,2 but it has to be a weak one, probably the weakest hydrogen bond in the crystal structure. The weak character of the hydrogen bonds in the E-ring is revealed also in the distribution of hydrogen bond distances where the maximum is at 220 pm. The maximum in the other parts of the molecule are about 180-200 pm, indicating clear, strong hydrogen bonds

Rg )

(

)

∑iri2mi 1/2 ∑imi

(1)

where mi is the mass of atom i and ri is the position of atom i with respect to the center of mass of the group. In the analysis of Rg, the center of mass of the group was chosen to be the center of all the 10 steroid molecules. In the budesonide simulations, Rg was converged to about 0.95, 2.40, and 2.20 in water, methanol, and ethanol simulations, respectively. This is somewhat unexpected because the amount of HBs in the methanol simulation is larger than in water simulations. This indicates that in the water environment, budesonide easily forms more aggregates to maximize the hydrophobic interactions and hydrogen bonding network with water. In the beclomethasone dipropionate, Rg was 1.0, 2.4, and 2.5 in water, methanol, and acetone simulations, respectively. The movements of side chains in steroid structures are interesting theoretically because these might affect the aggregation of molecules by hindering the formation of hydrogen bonds between solute molecules, for instance, forming complexes predicted by ab inito calculations. In the solvation of the budesonide, the hydroxyl group in the side chain could form a van der Waals interaction with the adjacent keto group as well as a hydrogen bond with the polar solvent molecules according to ab initio calculations. In the MD simulations, the dihedral angle (C20-C21-O31-H) was mainly at 180 degrees in all solvents forming hydrogen bonds to solvent. In water and methanol simulations, the distributions show that also 60 degrees was possible, forming perhaps some kind of complex with solvent, as predicted by ab initio calculations. The main point is that using only a single solvent molecule to describe the first solvation shell and its influence on structures may lead to wrong conclusions. In our MD simulations, we do not have any hydrogen bonding term in the force field, but the use of SPC water models have proved to achieve quite accurate results. Sillanpa¨a¨ et al. compared geometries of aluminum hydroxides using the supermolecule approach, continuum solvation model, and ab initio MD: They concluded that solvated species from the supermolecule approach might be totally different when compared to MD results.41 In beclomethasone, the distance between the hydrogen atom in the hydroxyl group and the ester group is ca. 0.6 and 0.5 nm in water and methanol solvents, respectively. In the crystallization of beclomethasone,

126

Crystal Growth & Design, Vol. 2, No. 2, 2002

a monohydrate form can be formed in which water is placed between a hydroxyl group and a keto group in the side chain. In the solution, the ester oxygen atom was closer (0.4 nm) to the hydroxyl group, indicating another form of solvation. In water environment, there is bigger flexibility in distances, up to 0.8 nm, and the dihedral angle (OdC-CH2-O) distribution shows that in the water environment the dihedral angle is 105 or 255 degrees, but in acetone only 105 degrees was detected. However, in our simulation only 1 ns was analyzed. In longer simulations, there is a possibility to see also other conformations; in any case, this shows the crucial effect of the solvent on conformations of molecules. It is worth noting that it could be interesting to test a polarizable force field in the calculations of budesonide. As is seen in the quantum chemical calculations, the electronic structure of the polar groups changes during the formation of a HB. The polarizable force fields could describe in some manner this phenomenon. Conclusion We have calculated geometries and stabilization energies for budesonide-water complexes with different methods. According to this study, the PM3 Hamiltonian gives fairly good geometries for budesonide--water complexes when comparing to complexes optimized at Hartree-Fock level with 6-31G** basis set. The AM1 optimized complexes differ dramatically from both PM3 Hamiltonian and Hartree-Fock optimized complexes. The AM1 Hamiltonian did not form hydrogen bonds but usual complexes. The AM1 stabilization energies are closer to HF stabilization energies, because in bifurcated structures there exist two proton donor-acceptor pairs that stabilize complexes. The PM3 Hamiltonian always underestimates stabilization energies as compared to HF stabilization energies. In the molecular dynamics simulations, budesonide and beclomethasone, clear differencies in aggregation were detected in different solvation. The conformation distribution of side chains was also different. This indicates that now it is possible to simulate the formation of aggregates and in the future perhaps actual nucleation using simple approximations. Acknowledgment. The authors thank the Center for Scientific Calculations (CSC) in Finland for providing the facilities for this work. References (1) Desiraju, R. G.; Steiner, T. The Weak Hydrogen Bond in Structural Chemistry and Biology; Oxford university Press Inc: New York, USA, 1999. (2) Albertsson, J.; Oskarsson, Å.; Svensson C. Acta Crystallogr. 1978, B34, 3027. (3) Duax, L. W.; Griffin, J. F.; Rohrer, D. C. J. Am. Chem. Soc. 1981, 103, 6705. (4) Thale´n, A. Acta Pharm Suec. 1987, 24, 97. (5) Thale´n, A. Acta Pharm Suec. 1982, 19, 327. (6) Thale´n, A.; Nylander, B. Acta Pharm. Suec. 1982, 19, 247. (7) Karelson, M. Advances in Quantum Chemistry, Academic Press Limited: New York, 1997; Chapter Quantum Chemical Treatment of Molecules in Condensed Disordered Media.

Pa¨iva¨rinta et al. (8) Wong, M. W.; Frisch, M. J., Wiberg, K. B. J. Am. Chem. Soc. 1991, 113, 4776. (9) Safi, B.; Choho, K.; De Proft, F.; Geerlings, P. Chem. Phys. Lett. 1999, 300, 85. (10) Miertus, S.; Scrocco, E.; Tomasi, J. J. Chem. Phys. 1981, 117, 55. (11) Klamt, A. J. Phys. Chem. 1995, 99, 2224. (12) Lahti, A.; Hotokka, M.; Neuvonen, K.; Karlstro¨m, G. J. Mol. Struct. (THEOCHEM) 1998, 452, 185. (13) Craw, J. S.; Guest, J. M.; Cooper, M. D.; Burton, N. A.; Hillier, I. H. J. Phys. Chem. 1996, 100, 6304. (14) Onsager, L. J. Am. Chem. Soc. 1936, 58, 1486. (15) Kirkwood, J. G. J. Chem. Phys. 1934, 2, 351. (16) Couthino, K.; Canuto, S. Advances in Quantum Chemistry, Academic Press Limited: New York, 1997; Chapter Solvent Effects from a Sequential Monte Carlo-Quantum Mechanical Approach. (17) Yasuda, T.; Ikawa, S. Chem. Phys. 1998, 238, 173. (18) Cummins, P. L.; Gready, J. E. J. Comput. Chem. 1996, 17, 1598. (19) Palmer, M. H.; Walker, I. C.; Guest, M. F. Chem. Phys. 1998, 238, 179. (20) Swaminathan, S.; Whitehead, R. J.; Guth, E.; Beveridge, D. L. J. Am. Chem. Soc. 1976, 99, 7817. (21) Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. Comput. Phys. Comm. 1995, 91, 43. (22) Taylor, R.; Kennard, O. J. Am. Chem. Soc. 1982, 104, 5063. (23) Mills, J. E. J.; Dean, P. M. CAMD 1996, 10, 607. (24) Ditchfield, R.; Hehre, W. J.; Pople, J. A. J. Chem. Phys. 1971, 54, 724. (25) Hehre, W. J.; Ditchfield, R.; Pople J. A. J. Chem. Phys. 1972, 56, 2257. (26) Hariharan, P. C.; Pople J. A. Theor. Chim. Acta (Berl.), 1973, 28, 213. (27) Frich, M. J.; Trucks, G. W.; Schlegel, H. B.; Gill, P. M. W.; Johnson, B. G.; Robb, M. A.; Cheeseman, J. R.; Keith, T.; Petersson, J. A.; Montgomery, K.; Raghavachari, K.; AlLaham, M. A.; Zakrzewiski, V. G.; Ortiz, J. V.; Foresman, J. B.; Peng, C. Y.; Ayala, P. Y.; Chen, W. J. S.; Defrees, D. J.; Baker, J.; Stewart, J. P.; Head-Gordon, M.; Gonzales, C.; Pople, J. A. GAUSSIAN94, Revision D.2, Gaussian, Pittsburgh, PA, 1995. (28) Wang, J.; Boyd, R. Chem. Phys. Lett. 1996, 259, 647. (29) Wang, J.; Boyd, R. J. Phys. Chem. 1996, 100, 16141. (30) Jorgensen, W. L.; Svensson, C. J. J. Am. Chem. Soc. 1985, 107, 569. (31) van Gunstren, W. F.; Billeter, S. R.; Eising, A. A.; Hu¨nenberger, P. H.; Kru¨ger, P. Mark, A. E., Scott, W. R. P., Tironi, I. G. Biomolecular Simulations: The GROMOS96 manual and user guide. Zu¨rich, Switzerland: Hochschulverlag AG an der ETH Zu¨rich, 1996. (32) van der Spoel, D.; van Buuren, A. R.; Apol, E.; Meulenhoff, P. J.; Tieleman, D. P.; Sijbers, A. L. T. M.; Hess, B.; Feenstra, K. A.; Lindahl, E.; van Drunen, R.; Berendsen, H. J. C. Gromacs User Manual version 2.0 Nijenborgh 4, 9747 AG Groningen, 1999. (33) Jorjensen, W. L.; Briggs, J. M.; Contreras, M. L. J. Phys. Chem. 1990, 94, 1682. (34) Freitas, L. C. G.; Cordeiro, J. M. M.; Laurenti, F. L. J. Mol. Liq. 1999, 79, 1. (35) Ferrario, M.; Haughney, M.; McDonald, I. K.; Klein, M. L. J. Chem. Phys. 1990, 93, 5156. (36) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (37) Petersen, H. G. J. Chem. Phys. 1995, 103, 3668. (38) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (39) Berendsen, H. C.; Grigera, J. R.; Straatsma T. P. J. Phys. Chem. 1987, 91, 6269. (40) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D. J. Chem. Phys. 1983, 79, 926. (41) Sillanpa¨a¨, A.; Pa¨iva¨rinta, J. Hotokka, M.; Rosenholm, J.; Laasonen, K. J. Phys. Chem. A 2001, 105, 10111.

CG010032D