J. Phys. Chem. B 2008, 112, 2085-2094
2085
Monte Carlo Structure Simulations for Aqueous 1,4-Dioxane Solutions Peter I. Nagy,*,† Gergely Vo1 lgyi,‡ and Krisztina Taka´ cs-Nova´ k‡ Department of Medicinal and Biological Chemistry and the Center for Drug Design and DeVelopment, The UniVersity of Toledo, Toledo, Ohio 43606-3390, and Department of Pharmaceutical Chemistry, Semmelweis UniVersity, Ho¨gyes E. u. 9, H-1092 Budapest, Hungary ReceiVed: July 17, 2007; In Final Form: NoVember 1, 2007
Monte Carlo simulations in the NpT ensembles have been performed for the structure exploration of aqueous 1,4-dioxane solutions. Three different systems with all-atom dioxane:TIP4P water molar compositions of 2:500 (code:D2), 8:465 (D8), and 17:425 (D17) modeled solutions of 0.22, 0.88, and 1.86 mol/dm3 concentrations, respectively, at T ) 298 K and p ) 1 atm. The calculated solution densities increase from 0.992 to 1.002 g/cm3 with increasing dioxane concentration and approach the experimentally determined densities within 1%. This close agreement was achieved by utilizing RESP charges fitted to the in-solution IEF-PCM/B3LYP/6-31G* electrostatic potential of dioxane taken in its chair conformation and recently developed C, H steric parameters for ethers for calculations with a 12-6-1 all-atom potential. Solution structure analyses pointed out that the dioxane molecules arrange in the solutions with favorable distances of 4-8 Å for the ring symmetry centers. Within this range not only pairs of rings but triangular triads and tetrads have also been observed with center-center distances 4.5 in a linear regression analysis, the calculated r2 is equal to 0.63. Although the 12-6-1 potential depends only on atom-atom distances, the overall interaction energy of two dioxane molecules at a given Xs‚‚‚Xs separation depends on the orientation of the rings as well. The importance of this latter effect reveals from a breakdown of the EXX term: the Coulomb component averaged over 0.1 M configurations is slightly positive in some cases. The dioxane-dioxane interaction, however, was always dominated by the negative 12-6 energy term accounting for the van der Waals interactions and keeping the overall EXX also negative. A necessary feature of a converged simulation is the stability of the macroscopic thermodynamic parameters. Figures S1 and
Nagy et al.
Figure 3. Changes of the ESX solute-solvent interaction energies for the D8 and D17 models in the last 45 and 75 M configurations, respectively. Values were sampled in 0.5 M steps averaged over 0.1 M configurations. For codes D17/1 and D17/2, see Figure 2.
S2 (Supporting Information) show that both the density and the system energy are stable after considering 105 M configurations. For the D17 system, the total energies averaged over the last 75 M configurations are -4690.3 ( 3.3 and -4684.4 ( 2.3 kcal/mol and do not differ significantly at the 95% (2σ) level from the two simulations. The corresponding densities (Table 1) are 1.002 ( 0.002 and 1.000 ( 0.001 g/cm3. Thus the two simulations provide statistically equivalent macroscopic system parameters for the D17 model. This conclusion is important with reference to the forthcoming analyses of the EXX and ESX (solute-solvent) energies. Figures 2 shows the solute-solute energies for the D8 and D17 models. The values were sampled in 0.5 M steps (averaged over the preceding 0.1 M configurations) in the 105-150 M and 105-180 M ranges. The D8 curve varies between -2 to -9 kcal/mol. The EXX terms are considerably more negative for the D17 system, which is reasonable by considering almost twice as many solutes in the latter model. The two EXX curves for the D17 model differ, however, conspicuously. Although the simulations started from different solute arrangements, the EXX curves run close to each other in the 120127.5 M range. This finding has been interpreted so that the different starting structures do not remarkably affect the statistical and structural results after considering 105-120 M configurations. After 127.5 M, however, the courses of the curves are different. Both curves show fluctuations by a few kcal/mol, but the EXX values generally differ by about 10 kcal/ mol. Since neither the densities nor the total energies of the systems deviate significantly (Figures S1, S2), the systems have been considered as two different microstates of the same macrostate (or those of two very close macrostates). The nearly equal total energies, Etot, can come into existence by compensation upon two other components according to the relationship Etot ) EXX + ESX + ESS, where ESS is the solvent-solvent interaction energy. The ESX curves are compared in Figure 3. The fluctuations of the curves in Figure 3 are acceptable, and the variations of the energy values do not indicate a systematic trend, which is just the character of an equilibrium state. The increase in ESX for the D17 system in the second simulation is just the indication of the partial EXX and ESX compensation. When EXX becomes more negative, then the solute molecules get closer to each other
Aqueous 1,4-Dioxane Solutions
Figure 4. O(ether)‚‚‚H(water) radial distribution functions for the D2, D8, and D17 models. The nontypical rdf was taken from the 135142.5 M range of the D17/2 simulation (for the code, see Figure 2).
and allow only a smaller number of the solutes to locate very close to some solutes. As a consequence, ESX becomes less negative. Solution Structures, General Analysis. The D2, D8, and D17 systems with calculated concentrations of 0.22, 0.86, and 1.86 mol/dm3, respectively, might be considered as characteristic models for aqueous dioxane solutions in the medium concentration range. The most interesting structural question regarding these solutions is whether the dioxane solutes form associations in the aqueous phase or are fairly evenly dissolved in the solvent. Radial distribution functions, rdfs, shown in Figure 4, were obtained by averaging over 7.5 M configurations. The selected ones are, however, typical and representative rdfs for D2, D8, and D17 models in the last 45-75 M configuration range. The maxima of the curves were found at 1.75-1.80 Å O(dioxane)‚ ‚‚H(water) separation for all three systems. These values are close to the optimized O(dioxane)‚‚‚H(water) distance in the dioxane‚‚‚water complex, where the O‚‚‚H separation and the O‚‚‚H-O angle were calculated at 1.72 Å and 179.5°, respectively, using the set 2 and TIP4P water parameters. The shift of the peak site by a few hundredths of an angstrom compared to the optimal O(dioxane)‚‚‚H(water) separation is a natural consequence of the thermal disordering and the effect of the interactions of the elements of the pairs with other molecules in the system. Integration of the four D2 rdfs (with only one representative in Figure 4) until their first minima at about 2.5 Å provide 0.9-1.5 water hydrogens around each dioxane oxygen in the first solvation shell. The average number is 1.16 over 45 M configurations. This O(dioxane)/H(water) coordination number provides an upper limit for the hydrogen bonds per dioxane oxygen. Only 12 O(dioxane)/H(water) rdfs (the maximum number of rdfs allowed by the BOSS program for a simulation session) out of the possible 16 and 34 were followed for each of the D8 and D17 systems, respectively, and one representative is provided for each system in Figure 4. For the calculated rdfs, the first peak values scatter generally in the 0.6-1.3 range. The O/H coordination numbers were calculated in the range of 0.61.6 for both systems, with most individual values close to the middle of the range throughout the last 45 and 75 M configurations for the D8 and D17 models, respectively. By comparison of the 135-180 M range of the two simulations for the D17
J. Phys. Chem. B, Vol. 112, No. 7, 2008 2089 model, where the EXX curves differ considerably, the O/H coordination numbers are still in the 0.57-1.55 range. The only exception was found in the 135.0-142.5 M range of the second simulation, where one of the O/H coordination numbers was as small as 0.37. The corresponding rdf has been indicated in Figure 4 as a nontypical curve. The changes in the heights of the rdfs suggest that the locations of water hydrogens around different ether oxygens are not equally strict. Nonetheless, no rdf has been found permanently with relatively low peak value. It is to be emphasized that the rdfs were averaged over 7.5 M configurations. If a less strict location of the water hydrogen appeared in a 7.5 M fraction of the simulation (lower peak), the corresponding rdf showed also higher peak values in later 7.5 M units. This finding thus suggests that the ether oxygens have not been consistently buried in a fairly hydrophobic pocket possibly formed by the C4 planes of close dioxane molecules in a triangular triad or a tetrad (see next section). Snapshots, for example in Figures 5 and 6, support this conclusion. Each typical rdf in Figure 4 has a narrow and sharp first peak and a broader and lower second peak at 3.0-3.2 Å. The first peak refers to the water hydrogen in hydrogen bond to the dioxane oxygen, whereas the second peak has been attributed to the other hydrogen of the same water molecule. The stable hydrogen bond keeps the other hydrogen at a preferred O‚‚‚H distance of about 3.0-3.2. The peak at about 4.6 Å for D17 and the plateaulike course of the rdfs in the 4.5-5.5 Å range for the D2 and D8 systems may indicate some structural differences for the systems with different dioxane concentrations. As mentioned, the first peak refers to a water‚‚‚dioxane intermolecular hydrogen bond. The 4.5-5.5 Å range of the rdfs has been attributed to the separation of the oxygen in question and the water hydrogen(s) in hydrogen bond with the other ether oxygen. Also this peak may have contributions from water molecules in the second solvent shell of the referred ether oxygen. The rdf values increase in the D2, D8, D17 series when the 4.5-5.5 Å range is considered. This means that the location of the water molecules around the ether solutes is stronger with increasing concentration. The resolved third peak for the D17 rdf supports this interpretation. With more solutes in the system, a specific water molecule experiences more, mostly favorable, interactions from the close solute molecules (see below). These interactions help more stably localize the water molecules in the first solvation shell of a solute, as suggested by the developed third peak in the D17 rdf. Number of the hydrogen bonds (nHB) could be, however, better estimated by analyzing the solute-solvent pair-energy distribution function (pedf) per dioxane molecule (Figure 7). Points of the pedf provide the average number of the water molecules in interaction energy with a solute in the range of E ( ∆E. ∆E was set to 0.25 kcal/mol in the present simulation, and each pedf in Figure 7 was averaged over the last 45 or 75 M configurations. Integral of the pedf until its first minimum provides the number of solvent molecules in strong interaction with the solute(s). In case of possible hydrogen bonds, this integral value has been interpreted as the nHB value for many systems.15bc,16a-c,18 The onset of the pedf for an infinitely dilute solution with a single solute cannot be at a more negative value than the optimized solute-solvent interaction energy. The threshold for the pedf is, in fact, generally at a slightly less negative value than the optimal one because the thermal disorder prevents the optimal solute-solvent interaction most of the time. In solution
2090 J. Phys. Chem. B, Vol. 112, No. 7, 2008
Nagy et al.
Figure 5. The arrangement of the D8 rings in configuration 148.5 M. Rings 2, 5, and 6 form a triangular trimer with R(Xs‚‚‚Xs) < 8 Å. The distances of the centers for rings 1, 3, and 7 are slightly larger than 8 Å. Rings 4 and 8 are more separated from the two triads.
models with at least two solute molecules, however, this restriction does not hold. The pedf calculates the interaction energy of each individual solvent molecule in interaction with all solutes. If two or more solutes are close enough to the selected solvent and most pair interactions are favorable (geometries with very unfavorable solute-solvent interactions cannot stably exist in an equilibrated system), then the overall one-solvent/all-solutes interaction energy may be more negative than the optimal interaction energy for the solute-solvent dimer. Indeed, the -11.0 kcal/mol threshold for the D2 pedf is considerably more negative than the optimized water-dioxane interaction energy of -5.65 kcal/mol. Integration of this pedf until its minimum at -3.0 kcal/mol provides about two water molecules strongly bound to a dioxane solute. In the present study, this integral value was accepted as providing the nHB value. The strongest solute-solvent interaction energies that actually emerged in the simulations are -12.0 and -14.0 kcal/mol for the D8 and D17 systems, respectively. The pedf integrals up to the first minima at E ) -3.0 kcal/mol provide nHB ) 2.11 and 2.06-2.18 dioxane‚‚‚water hydrogen bonds, respectively. The corresponding value is 2.03 hydrogen bonds for the D2 system. Overall, each dioxane solute forms about two hydrogen bonds, on average, to a water molecule in the 0.22-1.86 molar range. By consideration of the typical O(dioxane)‚‚‚H(water) coordination numbers, nHB ≈ 2 per solute correspond, on average, to one hydrogen bond per dioxane oxygen. Nontypical rdfs appearing through the D17/2 simulations provide a rationale for the small decrease of the number of the strongly bound solvent molecules to the more aggregated solutes in the D17/2 (nHB ) 2.06) compared to the D17/1 (nHB ) 2.18) simulation. However, nontypical rdfs were not found frequently in the second simulation, resulting in only a small decrease in the assigned nHB value. The shift of the thresholds for the D8 and D17 pedfs to more negative values is reasonable. The water:dioxane ratio of 250 for D2 decreases to 25 in the D17 model. In a more concentrated solution, a number of water molecules can stay close to more
than one dioxane solute. In a favorable hydrogen bond with the closest one, the dispersion interaction with some more remote dioxane molecules increases the solute-(one-solvent) interaction energy. This effect emerges in parallel with a slight increase of the hydrogen bonds per dioxane molecule and is in accord with the interpretation of the 4.5-5.5 Å region of the rdfs: water molecules in the first hydration shells of the dioxane solutes are more strongly localized for the D17 system than for the other two models. Solution Structures, Special Considerations. The distributions of the Xs‚‚‚Xs distances were calculated on the basis of 91 (D8) and 151 (D17) solution configurations sampled in 0.5 M steps in the last 45 and 75 M configurations, respectively. The distributions showed already fast convergence when only 46 and 61 elements were considered in subsets. There are 28 and 136 R(Xs‚‚‚Xs) distances for the pairs with 8 and 17 centers in the D8 and D17 systems, respectively. From the point-ofview of a possible dioxane aggregation, consideration of the shortest R values is important. The normalized average probability densities, p(RXX) ) (1/N)*(∆Ni(Ri)/∆Ri), (N ) 28, 136) for the center distances over 0.5 Å wide ∆Ri ranges are shown in Figure 8. The standard deviations for the ∆Ni(Ri) values are up to 10% with the exception of the hardly populated smallest and largest separations. The D8 probability density curve shows four resolved peaks. With a much larger number for the possible pairs, the D17 curves are smoother, showing only two resolved peaks and shoulders where the D8 curve has peaks. Each curve shows, however, a well-resolved first peak with local maxima at 5.75 and 6.25 Å for the D17 and D8 systems, respectively, and a first local minimum at 7.25 Å. Thus the favored value for the Xs‚‚‚Xs separation decreases by 0.5 Å for the closest dioxane centers when the concentration increases from 0.88 to 1.86 molar. Integration of the probability density function up to the first minimum indicates that there are about 3 pairs in this range, on average, for the D8 system, corresponding to 11% of the total numbers of the pairs.
Aqueous 1,4-Dioxane Solutions
J. Phys. Chem. B, Vol. 112, No. 7, 2008 2091
Figure 6. (a) The arrangement of the D17 rings in configuration 148.5 M of the first simulation. Rings (5, 7, 9), (5, 7, 15), (5, 9, 15), (7, 9, 15), and (8, 11, 16) form triangular trimers with 8.0 Å side limits. Rings (5, 7, 9, 15) form a tetrad. Most of the center separations for the pairs of rings 3, 10, 13, 14, and 17 are in the range of 7-10 Å. (b) The arrangement of the D17 rings in configuration 172 M of the second simulation. There are five triads (3, 11, 14), (4, 6, 7), (11, 14, 15), (11, 15, 16), (14, 15, 17) and a tetrad (1, 12, 15, 16) with R(Xs‚‚‚Xs) < 8 Å.
Statistics for individual solutes indicate one isolated dioxane ring, five rings with about one neighbor, and two rings with 1.5-1.8 neighbors with R < 8.0 Å for the D8 system, on average. If a ring has only one neighbor within the limit than the local structure depends on the neighbor number for the partner. If the partner ring has also only about one neighbor, then this pair of rings may be considered as forming a dimer with respect to the R limit. Upon the calculated distance matrix of the centers, appearance of so-called triangular triads (also referred to as N3 structures henceforth) and tetrads (N4) were analyzed. Three and four rings form a triangular triad and a tetrad, respectively, if the distance
of the centers for any pair of the involved rings is below an R limit. By setting R to 8.0 Å, three rings with two neighbors for each (ring numbers 2, 5, and 6, see Figure 5) formed an N3 structure in 20 snapshots (22%), and other triangular trimers (involving rings 1, 3, and 7 as well) were found in 3 further snapshots (3%). Since the triangular trimers were formed and disrupted back and forth, their elements could not move too far away from each other. Figure 5 shows that the rings 1, 3, and 7 are close to each other and that ring 7 is close to rings 2 and 5 as well. These six close rings (including ring 6) have formed three different N3 structures in the sampled configurations. Rings 4 and 8 were never found participating in a triangular
2092 J. Phys. Chem. B, Vol. 112, No. 7, 2008
Figure 7. Solute-solvent pair energy distribution functions averaged over the last 45 and 75 M configurations for the D2, D8, and the D17 models, respectively. For codes D17/1 and D17/2, see Figure 2. The dN/dE curves for the D17 systems have almost zero values in the E ) -12 to -14 kcal/mol range, and this tail has not been indicated.
Figure 8. Numerical calculation of the probability density functions for the R(Xs‚‚‚Xs) separations of the symmetry centers with 0.5 Å wide ranges for the R values. The functions were obtained by utilizing distances from 91 sample configurations out of the last 45 M structures (D8) and from 151 sample configurations out of the last 75 M structures (D17). For codes D17/1 and D17/2, see Figure 2.
triad. Thus altogether in 25% of the sampled configurations there was an N3 structure, but none of snapshots showed an N4 structure. By consideration of the average numbers of the closest neighbors, the remaining 75% of the configurations without a triangular trimer may be still interpreted as different combinations of close dimers and bent trimer(s) (with a distance larger than the R limit for two involved centers) or loose zigzag chains of the ring centers and a few isolated rings. The average solution volume per dioxane molecule is 1886 Å3 in the D8 system, which corresponds to a cube with an edge of about 12.4 Å. In a latticelike arrangement with the dioxane symmetry center located at the center of the closely packed cubes for a very large volume of the solution in three dimensions, the closest-neighbor distance would be 12.4 Å. The calculated pair-distance distribution indicates, however, a number of much shorter separations for the centers, indicating
Nagy et al.
Figure 9. Percentage distributions (conf %) of the numbers of triangular trimers (triads) per configuration, calculated by considering 151 solution structures sampled in 0.5 M steps throughout the last 75 M structures. R ) 6.0 and 8.0 Å refer to the side limits for the triads. For codes D17/1 and D17/2, see Figure 2.
a trend for aggregation for a large fraction of the solute molecules in the 0.88 molar aqueous dioxane solution. Similar statistical analyses have been performed for the D17 system. At this concentration, the per dioxane volume corresponds to a “solution cube” with edges of about 9.6 Å. Thus for the D17 system, the trend for aggregation would require significantly shorter closest-neighbor distances than 9.6 Å. Although the two p(RXX) curves for the D17 systems are remarkably different in the shortest Xs‚‚‚Xs region (Figure 8), the results characterize microstates underlying very similar macrostates, as mentioned above. The distribution function from the second simulation has a peak value at 5.75 Å, twice as high as from the first simulation. By consideration of both simulations, 15-25 dioxane pairs meet the limit of R < 8 Å, on average, out of the total of 136 combinations. This result indicates that at least a large fraction of the 17 solute molecules forms dimers with separations shorter than in a regular lattice. (The pairs necessarily include a number of rings at relatively short separation from more than one other rings. 17 solute molecules could form only 8 distinct pairs.) The average numbers of the neighbors not farther away than 8 Å are 1-4 for the most solutes. The analysis results for the N3 structures are summarized in Figure 9. The total number of the triads as well as the percentage distribution of the number of the N3 structures found in a specific solution configuration depends strongly on the R limit. Increasing the limit from R ) 6.0 Å to R ) 8.0 Å, the total number of the N3 structures increased from 8 to 285 and 173 to 1702 from the first and the second simulations, respectively, corresponding to 0.1 to 1.9 and 1.1 to 11.3 triads, on average, per configuration. (The considered number of the configurations was 151 in all cases.) However, the actual number of the triads in the configurations is much different. The percentage distributions of the number of the N3 structures (Figure 9) show that no triad exists in 4295% of the snapshots with R ) 6.0 Å, whereas this fraction has reduced to 0-23% with R ) 8.0 Å. Thus the latter limit allows at least one triad for at least 77% of the samples. The distribution functions changed in shape from a monotonically decreasing form to a curve with a maximum when the R limit was increased. Integrals of the corresponding curves nearly in the range of the middle-third of the calculated N3 values
Aqueous 1,4-Dioxane Solutions indicate a considerable shift in the number of configurations. By consideration of the second simulation, 28% of the configurations contain 2-4 triads (N3max ) 6) with R ) 6.0 Å, whereas 61% of the configurations include 8-14 triads (N3max ) 21) with R ) 8.0 Å. From the first simulation, 33% of the configurations include 2-3 triads (N3max ) 5) with R ) 8.0 Å. Reducing the limit to R ) 6.0 Å, only 5% of the configurations included one triad, and no triad was found in 95% of the cases. Numbers of the tetrads are reasonably smaller than triads; two tetrads were found at a maximum in any configuration. No tetrad was assigned with a Xs‚‚‚Xs limit of 6.0 Å, and tetrads appear only at the R limit of 8.0 Å from the first simulation. Even at this limit, 99% of the configurations remained without a tetrad, and a single tetrad was found only for two snapshots out of 151 ones. In snapshots from the second simulation, 41 and 112 tetrads were assigned at R ) 7.0 Å and R ) 8.0 Å, respectively. Sometimes two tetrads appeared in a configuration, and in total 26 and 58% of the sampled structures included tetrads meeting the R limits of 7.0 and 8.0 Å, respectively. Thus, formation of tetrads becomes feasible with an R limit larger than 6.0 Å. Tetrads then may appear in a considerable fraction of the configurations, but our calculations predict two dioxane tetrads per configuration at most. Parts a and b of Figure 6 may provide some illustration of the analysis results above. Figure 6a shows the solute configuration 148.5 M from the first simulation. By use of the indicated ring numbering, the solutes could be assigned to three groups for characterizing the solute structures in the sampled 151 configurations. Elements of both group 1 (rings 1, 2, 5, 7, 8, 9, and 15) and group 2 (rings 3, 10, 13, 14, and 17) formed N3 structures primarily with rings assigned to the same group. Rings 4, 6, 11, 12, and 16 (group 3) form triangular trimers with each other very rarely but are sometimes involved in N3 structures with group 1 rings. Whereas group 1 and 2 rings rarely form mixed triads, group members could comprise more than one triangular trimer at a time. The solute configuration in Figure 6a corresponds to N3max ) 5 triads (R ) 8.0 Å) composed of the (5, 7, 9), (5, 7, 15), (5, 9, 15), (7, 9, 15), and (8, 11, 16) rings. The elements of the first four triads formed by the group 1 rings comprise a tetrad with Xs‚‚‚Xs < 8.0 Å. The fifth triad is an example for a “mixed” N3 structure by group 1 and 3 elements. The maximum of five N3 structures were found in eleven configurations with 2-3 trimers involving rings from the same group. The composition of the N3 trimers within a group changed but returned from time to time. This means that, e.g., rings 5 and 7 (group 1) or ring 17 (group 2) formed N3 structures with Xs‚‚‚Xs distances
