Molecular dynamics simulations of 18-crown-6 aqueous solutions

Journal of Molecular Liquids 224 (2016) 825–831

Contents lists available at ScienceDirect

Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq

Molecular dynamics simulations of 18-crown-6 aqueous solutions Yuriy G. Bushuev ⁎, Tatiana R. Usacheva, Valentin A. Sharnin Department of General Chemical Technology, Ivanovo State University of Chemistry and Technology, Ivanovo, Russia

a r t i c l e

i n f o

Article history: Received 11 July 2016 Received in revised form 5 October 2016 Accepted 12 October 2016 Available online 17 October 2016 Keywords: Molecular dynamics 18-Crown-6 Water Aqueous solutions Diffusion

a b s t r a c t Physicochemical, thermodynamic and structural properties of 18-crown-6 aqueous solutions were investigated by molecular dynamics simulations. Force fields for flexible solute and solvent molecules were used. The calculated values of the density, hydration enthalpy, self-diffusion coefficients and their concentration dependences were found to be close to experimental data. The structure of hydration shells of crown ether molecules and the structure of water in the solutions were determined. A residence lifetime distribution of water molecules from the first hydration shell of a crown ether molecule shows a presence of long-lived states with typical lifetimes of hundreds picoseconds. © 2016 Elsevier B.V. All rights reserved.

1. Introduction The investigation of complexation in the macrocyclic compound solutions has attracted large attention from scientists due to their applications in various fields [1]. 18-Crown-6 (18С6) is one of the most extensively studied macrocyclic polyethers [2–9]. Due to high flexibility of a molecule and a presence of an intramolecular cavity, which is able to contain ions and molecular fragments, an 18С6 molecule can form complexes with ions and various organic molecules. The cavity is formed by the oxygen atoms carrying excess negative electrical charge. In aqueous solution this molecular fragment manifests hydrophilic properties. The outer surface of the molecule is mostly formed by hydrophobic molecular moieties \\CH2\\CH2\\. An 18С6 molecule adopts 190 ideal diamond-lattice conformations [10] but two of them are the most important. In crystal hydrates the molecule adopts the D3d conformation, which is characterized by the largest opening of the intramolecular circular cavity [11]. In the anhydrous 18С6 crystal the molecules are found in the Ci conformation, where atoms are more densely packed, and the cavity is closed [12]. In solutions a population of the conformational states is determined by dynamic equilibrium and depends on energy, thermodynamic conditions (T, p), the nature of the solvents, and a complexation of 18C6 with ions or organic molecules [13,14]. Aqueous solutions of 18С6 are important for investigation because they are model systems for study of the complexation in solutions of biological macromolecules [1]. Computer simulations allow to characterize the studied systems at atomistic scale and structural, thermodynamic properties of solutions can be determined. Aqueous solutions of ⁎ Corresponding author at: Ivanovo State University of Chemistry and Technology, pr. Sheremetevskiy, 7, Ivanovo, Russia. E-mail address: [email protected] (Y.G. Bushuev).

http://dx.doi.org/10.1016/j.molliq.2016.10.062 0167-7322/© 2016 Elsevier B.V. All rights reserved.

crown ethers have been simulated previously [3–6,8,15–17]. However, due to the complex structure and flexibility of the 18С6 molecule, modeling of crown ether solutions does not an easy task. It is difficult to take into account all types of intramolecular movements, and often the molecule or specific fragments of the molecules (CH2 groups, for example), are considered as solids [5,13,18,19]. The simplifications may compromise the quality of the computer models. To obtain a representative set of molecular configurations for completely flexible solute and solvent molecules it is necessary to simulate large systems consisting of thousands of atoms. Due to slow intramolecular relaxation the dynamics must be run for a few nanoseconds. A separate problem is a selection of charges, qi, assigned to atoms of an 18С6 molecule. This problem has been discussed recently [4,8]. Y. Sun and P. Kollman reported [20] that an 18C6 molecule adopts 30 low energy conformations (ΔE b 10.5 kJ mol−1) in a gas phase. The authors used the AMBER force field supplemented with the parameters of Billeter and Howard [21]. According to the force field, atoms of 18C6 have atomic charges: qO = −0.406, qC = 0.244, and qH = −0.021 ǀeǀ [22]. This model, denoted 18C6(−0.4), has been widely used for simulations of 18C6 containing systems [4,8,13,17,21,22]. The charges qO = −0.340, qC = 0.086, and qH = 0.042, obtained from an ESP HF/6-31G (d,p) calculation, have been shown to satisfactorily predict the complexation of K+ in water and in methanol [4]. However, it was demonstrated that the stability of the Na+, K+ and Cs+ crown ether complexes in the gas phase and in water is underestimated when these charges are used [8]. To mimic the polarization of 18C6 by ions the charges from the 18C6(− 0.4) model were scaled by a factor of 1.5 [4,8]. This model is named 18C6(− 0.6). Aqueous solutions of the 18C6 play a role of reference systems, when complexation in multicomponent aqueous solutions is investigated. The basic properties of the 18C6(− 0.6) model are less known in comparison with the more studied 18C6(−0.4) model.

826

Y.G. Bushuev et al. / Journal of Molecular Liquids 224 (2016) 825–831

The main purpose of this study is to determine the physicochemical and thermodynamic properties of aqueous solutions by molecular dynamics simulations of the two models, where both water and 18С6 molecules are fully flexible. To examine the quality of the models calculated and available experimental data were compared. The structure of 18С6 hydration shells, densities of the solutions, the enthalpy of hydration, self-diffusion coefficients of the molecules, and energies of intra- and inter-molecular interactions are discussed in the next sections. Particular attention is focused on a collective movement of ether-water molecular complexes.

the force fields are used for optimization, but these values significantly differ from the result of DFT calculation (−15.2 kJ mol−1). The force fields have been especially designed to simulate condensed phases of the substances. They effectively include molecular polarization and multibody interactions. In comparison with TIP3P the flexible SPC/Fw water model better reproduces a self-diffusion coefficient of liquid water [25]. It was an additional reason to select the SPC/Fw model. 3. Results and discussions 3.1. Structure of hydration shells

2. Methods Classical molecular dynamics (MD) simulations have been employed to investigate 18C6 aqueous solutions. Production simulations were performed in NpT ensemble at T = 300 K and p = 0.1 MPa. All calculations were performed with the DL_POLY Classic 1.9 program [23]. The Nosé-Hoover (Melchionna) thermostat and barostat with relaxation times of 0.1 ps (T) and 1 ps (p) was employed. An integration timestep was of 1 fs. The cubic boxes, taken with the periodic boundary conditions, contained 1200 water and 1, 3, 6, and 12 crown ether molecules. The cutoff for short-range interactions was set to 12 Å. Electrostatic interactions were treated using the Ewald method. After the relaxation period lasting ca. 1 ns, molecular configurations were saved for later analysis. The total simulation time of each system usually exceeds 3 ns. The AMBER force field, a well-proven in the simulations of large organic molecule, was applied to calculations of intramolecular interactions of an 18C6 molecule [24]. The potential energy of the molecule is calculated by the equation: 2

U ¼ ∑ kl ðl−l0 Þ þ ∑ kθ ðθ−θ0 Þ2 þ ∑ ∑ V n ð1 þ cosðnφ−γÞÞ bond angle dihedrals n "
 6
 12 # Rij Rij qi q j ; ð1Þ −2ε ij þ εij þ∑ Rij Rij Rij ib j

An 18С6 molecule is very flexible and in aqueous solution it adopts various conformations. The structure of hydration shells of the 18С6 molecule depends on the conformational state of the molecule, which has hydrophobic and hydrophilic fragments. It has been analysed, using different criteria for definition of the shells [15,18]. In the present work we have focused on hydration of strongly interacted oxygens of 18С6 molecules. Two examples of a molecule hydration are presented in Fig. 1. These molecular configurations were selected from the calculated trajectory. The central-symmetric configuration is shown in Fig. 1a, where the 18С6 molecule adopts the D3d conformation. The molecule forms six H-bonds with water molecules located above and below the plane of the macrocycle. This configuration is typical for crystal hydrates [11]. The other configuration is shown in Fig. 1b, where some oxygen atoms of the 18С6 molecule are in open positions and, therefore, more accessible to hydration. The oxygens have a large electric charge and able to interact strongly with water molecules, playing a role of hydrogen acceptors, when intermolecular H-bonds are formed. Standard functions, which describe the structure of liquids, are the radial distribution function (RDF), gij(r), and the running coordination number, Nj(r): Zr g ij ðr Þ ¼ ρ j ðr Þ=ρav

g ij ðr Þr 2 dr

N j ðrÞ ¼ 4πρav

ð2Þ

0

which includes contributions from the distortion of the lengths of valence bonds (l0), distortions of valence (θ0) and dihedral angles, and the contributions from electrostatic and van der Waals interactions. The 1–4 van der Waals and Coulombic interactions were decreased by factor 2.0 and 1.2, respectively. There are 42 chemical bonds, 78 valence angles and 132 dihedral angles in an 18С6 molecule. Therefore, to calculate the potential energy of the molecule it is necessary to use dozens of parameters kl, kθ, Vn, Rij⁎, εij, n, and γ. This is a difficult task, so usually a part of the terms was not taken into account. It is necessary to apply special programs simplifying the process of writing the FIELD file for the DL_POLY package. The DL_FIELD version 3.2 program has been employed to write the force field for the crown ether molecules and the force field for intermolecular interactions. To calculate the interactions of water molecules the SPC/Fw potential [25] was selected, which takes into account the flexibility of water molecules. The parameters of the 18C6 force field are presented in Tables S1, S2 (Supplementary information) and in a fragment of a DL_POLY OUTPUT file (S3). To verify the force fields, geometrical and energetic characteristics of D3d single out monodentate and bidentate water-ether motifs (Fig. S1) were optimized for the three charge distributions on an 18C6 molecule (Table S2). The TIP3P and SPC/Fw water models were used for testing. Results of calculations are presented in Table S3. The geometrical parameters of the water-ether complexes weakly depend on water model. The bidentate complex is more stable than the monodentate one for all investigated models and the 18C6(−0.6) model demonstrates the largest stability. For all these models intermolecular oxygen-oxygen distances are shorter and energy gaps are larger than the corresponding values calculated by the DFT method [10]. The energy gap between the bidentate and monodentate motifs is varied from −25.8 to −32 kJ mol−1 when

where ρj(r) is the local number density of j atoms in a spherical layer of thickness dr at the distance r from the i atom, ρav is the average density. Among variety of atom-atom RDFs, which describe hydration shells of 18С6 molecules and the solvent structure, the functions calculated for oxygen atoms of water molecules (Ow) with respect to oxygen atoms of 18С6 molecules (Os), gOsOw, and oxygens of water, gOwOw, are particularly important. Fig. 2 shows the gOsOw and NOsOw functions calculated for both the 18C6(−0.4) and 18C6(−0.6) models. To calculate the running coordination numbers of the hydrophilic moiety of 18C6 a special computer program was used. According to Eq. (2), N(r) describes a local

Fig. 1. Snapshots of molecular configurations illustrating hydration of 18С6 molecules: (a) the compact structure which is typical for crystal hydrates; (b) the conformer with open access to the oxygen atoms. The distances from 18С6 oxygens to water oxygens (red spheres) do not exceed 3.2 Å. All hydrogen atoms are not shown in (b) for clarity. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Y.G. Bushuev et al. / Journal of Molecular Liquids 224 (2016) 825–831

827

8

a

1.0

b

18C6 (-0.6)

7

m=0.278 mol kg-1 18C6 (-0.4) 18C6 (-0.6)

0.5

NOsOw

gOsOw

6 5 4 3

m, mol kg-1 0.046 0.139 0.278 0.555

2

18C6 (-0.4) 1 0.0 2

4

6

8

10

12

14

0 2.4

2.6

2.8

r, Å

3.0

3.2

3.4

3.6

r, Å

Fig. 2. (a) Radial distribution functions, g(r), and (b) corresponding running coordination numbers, N(r), for water oxygens (Ow) with respect to oxygens of 18С6 (Os). The vertical lines mark the position of specific points of the RDFs. The horizontal lines show the corresponding N(r) values.

0

g ðr0Þρð0Þ ¼ g ðrmÞρðmÞ; gðr0Þ ¼ g ðrmÞ ¼ g ðrmÞ

V s ðmÞ ; Vw

ð3Þ

where m is the ether concentration, Vs(m) and Vw = Vs(0) are the volumes of the simulation boxes for the solutions and water, respectively. All the boxes contain the same number of water molecules. Fig. 4b shows that the renormalized RDFs do slightly depend on concentration at r b 3.2 Å. Therefore, the analysis of the functions has shown only a small change of water structure in these aqueous solutions. Approximately a half of water molecules belongs to the first hydration shells of 18C6 molecules at m = 0.56 mol kg−1, but the solute only slightly distorts the structure of water. 3.2. Density and self-diffusion coefficients These two physicochemical properties of solutions are important to estimate the quality of the models. The experimental [7] and calculated densities of the aqueous solutions of 18C6 are presented in Fig. 5. The calculated density of SPC/Fw water is about 1% larger than the experimental value, but any water model predicts some properties with errors [25]. The presented functions show almost linear increase in density of the solution with increase in concentration of 18С6. The slopes of the lines, fitting experimental and calculated data, are only slightly different. The calculated densities were shifted down on the value of the

2.73

16

m=0, pure water m=0.56 mol kg-1

12

3.2 2

8

NOw

3

gOwOw

environment of each Os, but some water molecules belong to coordination shells of several Os atoms of the same 18C6 molecule. The hydration shells of atoms are overlapped and this effect is especially pronounced at large r. It is important to know a number of water molecules near any 18С6 molecule, and we took into account only the water molecules, which are in a sphere of radius r from any Os atom of the same 18C6 molecule. Two well-resolved peaks of gOsOw are observed at 2.75–3.0 Å and 4.7 Å. In average 2.5–3 water molecules are in close contact with ether oxygens. The first hydration shell is determined by the position of the first minimum of gOsOw, which is at 3.2 and 3.5 Å for the 18C6(− 0.6) and 18C6(− 0.4) model, respectively. In a dilute aqueous solution, the first hydration shell of the hydrophilic moiety of an ether molecule contains ca. 5.5–6.1 water molecules, and this number decreases with concentration due to solute-solute interactions resulting in overlap of the hydration shells. This effect is more pronounced for 18C6(−0.4) model. The position of the second minimum of gOsOw at ca. 6 Å corresponds to a boarder of the whole first hydration shell of an 18C6 molecule, which are hydrated approximately by 62 water molecules. Beyond discussed specific points, no additional features of any significance can be distinguished in Fig. 2a. Functions gOwOw, presented in Fig. 3, show only small changes of the pure water structure in the concentrated aqueous solution of 18C6(−0.6). At r b 3.2 Å a number of nearest neighbours of a water molecule in water and in solutions is slightly larger than in the ice Ih and Ic, where it is equal to four. A degree of tetrahedral ordering is characterized by position and height of the second maximum of gOwOw. In this region (r ~ 4.5 Å) RDFs of pure water and water in solutions are very close to each other. The difference is noticeable in the values NOw(r) outside the first hydration shell (r N 3.2 Å, Fig. 3). This is the effect of excluded volume. The first peak of gOwOw is higher for the concentrated solution than for pure water. This may be associated with an increase in the degree of water structuring, but the running coordination numbers, calculated for pure water and solution, are almost equal at r b 3.2 Å (Fig. 3). This is possible when the product g × ρav (Eq. (2)) does not depend on concentration. In average, a local environment of each water molecule is almost the same in pure water and in solutions. Reducing the average number density of water in the solutions is compensated by increasing of RDF. This conclusion is illustrated in Fig. 4 where the first peaks of the gOwOw functions for pure water and the solutions are presented. To demonstrate the compensatory effect, the RDFs were renormalized according to the equations:

1 4

4.05 1.5

0 2.4

2.8

3.2

3.6

4.0

4.4

0 4.8

5.2

r, Å Fig. 3. Radial distribution functions, g(r), and running coordination numbers N(r) for water oxygens in pure water and concentrated aqueous solutions of 18С6, gOwOw. N(r) curves are represented by red lines. The vertical lines mark the position of specific points of the RDFs. The horizontal red lines show the corresponding N(r) values. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

828

Y.G. Bushuev et al. / Journal of Molecular Liquids 224 (2016) 825–831

a

3

pure water m=0.28 mol kg-1 m=0.56 mol kg-1

3

b pure water m=0.28 mol kg-1 m=0.56 mol kg-1

2

gOwOw

g'OwOw

2

1

1

0 2.4

2.6

2.8

3.0

3.2

3.4

0 2.4

2.6

r, Å

2.8

3.0

3.2

3.4

r, Å

Fig. 4. (a) The first peaks of RDFs gOwOw calculated for water and 18C6(−0.6) solutions. (b) The RDFs renormalized according to Eq. (3).

difference between the calculated and experimental densities of pure water. The 18C6(− 0.4) model corresponds better to experimental data than 18C6(−0.6). It is noteworthy that the corrected densities, calculated for the 18C6(− 0.6) model, and the density of crystal hydrate with composition 1:12 (18С6/H2O) [11] lay on the fitting line. We may conclude that deviations from experimental data are not large and both models are reasonable. When an equilibrium state is reached, the self-diffusion coefficient D is computed by taking the slope of mean-squared displacement (MSD) at long time * + N 1 d 2 lim ∑ ðr ðt Þ−r i ð0ÞÞ D¼ 6N t→∞ dt i¼1 i

ð4Þ

where N is number of molecules, ri(t) is the position of the molecule i at time t. Fig. 6a shows the MSD dependence on time for systems containing six and twelve 18С6 molecules in the simulation boxes. Observed deviations of the curves from the straight lines indicate the presence of a high statistical noise due to a small number of 18С6 molecules in the simulation boxes (Eq. (4)). Nevertheless, the calculated diffusion coefficients for the centre of mass of an ether molecule are: D = 0.47 ± 0.08 (m = 0.28 mol kg−1) and 0.41 ± 0.07 × 10−9 m2 s−1 (m = 0.56 1.05 Calc. corrected 18C6(-0.6) Expt. Crystal Hydrate 1:12 Linear fit

1.2

d

1.04 1.1

d, g cm-3

1.03

1.0 0

1

2

3

4

5

m

1.02

mol kg− 1) for 18C6(− 0.6), and D = 0.43 ± 0.06 × 10− 9 m2 s−1 (m = 0.56 mol kg−1) for the 18C6(−0.4) model. Thus, within statistical uncertainties both models predict the same value of the diffusion coefficient. The experimental coefficient, measured for a dilute aqueous solution, is 0.56 × 10− 9 m2 s− 1 [26,27]. With increasing crown ether concentration a decrease in D is observed [27]. We estimated the coefficient from the plots presented in the article [27], D = 0.26 × 10−9 m2 s− 1 at m = 1 mol kg− 1. Thus, the calculated diffusion coefficients close to experimentally measured ones. The diffusion coefficients of water molecules are several times larger than the diffusion coefficient of 18С6 molecules. In pure water it is equal to 2.3 × 10−9 m2 s−1 [28], whereas the calculated value of 2.4 × 10−9 m2 s−1 is only slightly larger than experimental one. This characterizes the quality of the used water model. Water diffusion decreases with increase in concentration of 18С6 as it is presented in Fig. 6b. These large organic molecules strongly influence on diffusion of water molecules which are in the first hydration shells. In order to determine the nature of water molecules movement, distributions of their residence lifetimes near hydrophobic and hydrophilic moieties of 18C6 molecules were calculated. During the lifetime the 18С6 and water molecules are moving together. This type of movement is different from the recently studied kinetics of 18C6 conformational transition or lifetimes of H-bonds in water-ether adduct, which take place on a time scale of several picoseconds [10]. Each molecular configuration stored with a time step of 10 ps, was analysed to calculate the residence lifetimes. A list of water molecules which are at a distance rOsOw b 3.2 Å from any Os atom of the same water molecule was composed for each time tl. Then, intervals of continuous location (the residence lifetime, ti = tj − tk) of each water molecule in the first hydration shell of the hydrophilic moiety of each 18C6 molecule were determined. The other criterion was used for selection of water molecules belonging to the first hydration shell of the hydrophobic moiety: rOsOw N 3.2 Å and rCOw b 4.2 Å. The fractions F(t) are calculated according to Eq. (5):

1.01

1.00

0.99

F ðt i Þ ¼ Nðt i Þ=∑i Nðt i Þ;

Calc. 18C6(-0.6) Calc. corrected 18C6(-0.6) Calc. 18C6(-0.4) Calc. corrected 18C6(-0.4) Expt.

0.0

0.1

0.2

0.3

0.4

0.5

0.6

m, mol kg-1 Fig. 5. The experimental [7] and calculated and corrected densities of the aqueous solutions of 18С6. The density of crystal hydrate (1:12) and corrected densities of the 18C6(−0.6) model are presented in the insertion.

ð5Þ

where the N(ti) is a number of water molecules with the residence lifetime ti. The distributions of the fractions versus the lifetime are shown in Fig. 7. A large proportion of the water molecules has short lifetimes in the first hydration shell of 18C6 (ca. 10–30 ps). The fractions of long-lived states are small, and they decrease with increasing the lifetime. However, it was found that the distribution is very broad, and some water molecules have the lifetimes of hundreds picoseconds. The longest lifetime

Y.G. Bushuev et al. / Journal of Molecular Liquids 224 (2016) 825–831

829

70

MSD/(6N), Å2

50 40 30

b

18C6 (-0.6) m=0.56 Linear fit m=0.28 Linear fit 18C6 (-0.4) m=0.56 Linear fit

Calc. Expt. Linear fit

2.4

2.2

D x 109, m2s-1

a 60

2.0

20 1.8

10 1.6

0 0.2

0.4

0.6

0.8

1.0

1.2

1.4

0.0

1.6

0.1

0.2

0.3

0.4

0.5

0.6

m, mol kg-1

t, ns

Fig. 6. (a) Time evolution of mean square displacement of 18С6. (b) The diffusion coefficients of water in 18C6 aqueous solutions.

was N1.3 ns. The time of productive runs was a few nanoseconds, and it is not possible to determine the exact fraction of such states due to poor statistics. Water molecules from the hydration shell of hydrophobic moieties have residence times in order of magnitude smaller than the water

a

b

0.008

0.6

400

600

800

1000

0.01

1200

1400

t, ps

m, mol kg-1 0.046 0.139 0.278 0.555

0.2

Fraction

0.000 200

0.4

0.02

0.6

0.046 0.139 0.278 0.555

0.004

Fraction

molecule from the hydrophilic shell, ca. 0.17 ns. The same tendencies are observed for the 18C6(− 0.4) model (Fig. 7c). An 18C6 molecule slows down an exchange of water molecules between the first and second hydration shells and this effect is more pronounced for the hydrophilic moiety of the 18C6 molecule. This result correlates with the

0.00

0.4

60

80

100

120

140

160

t, ps m, mol kg-1 0.046 0.139 0.278 0.555

0.2

0.0

0.0 10

20

30

40

50

60

70

80

90

100

10

20

30

40

50

t, ps

t, ps 0.8

c Fraction

0.6

Fraction

0.002

0.001

0.000

200

0.4

300

400

500

600

t 18C6(-0.4) m=0.278 mol kg-1 Os-Ow C - Ow

0.2

0.0 10

20

30

40

50

60

70

80

90

100 110 120

t, ps Fig. 7. The distribution of water molecules versus their residence lifetime: (a) for the molecules from the first hydration shell of oxygens; (b) for the molecules from the first hydration shell of carbons, the 18C6(−0.6) model is used; (c) for molecules from the hydrophobic and hydrophilic moieties, the 18C6(−0.4) model is used.

830

Y.G. Bushuev et al. / Journal of Molecular Liquids 224 (2016) 825–831

-40

Δ H, kJ mol-1

-80

-120

-160

Calc. 18C6 (-0.6) Calc. 18C6 (-0.4) Expt.

-200

-240 0.0

0.1

0.2

0.3

0.4

0.5

0.6

m, mol kg-1 Fig. 8. Calculated (squares) and experimental (dashed line) [7,30] hydration enthalpies of 18С6. Bars show the standard error of mean only for the 18C6(−0.6) model. They are the same for the 18C6(−0.4) model.

experimentally observed decrease of water diffusion and with the increase of viscosity of solutions with increase in concentration of 18С6. [29]. 3.3. Energetics Enthalpy of hydration is defined by changes in enthalpy of molecule when it transfers from the gas phase into solution. The enthalpies of hydration (ΔH) are calculated by the equation:

excess enthalpy on the 18C6 concentration [7] were used to present the experimental enthalpies in Fig. 8. The enthalpies, calculated for 18C6(− 0.6) model, deviate from the experimental line, but the line are within the standard error intervals. The enthalpies of hydration are systematically shifted up in the case of the 18C6(− 0.4) model. The error intervals are the same for both models. The force field used for the calculations contains several terms (Eq. (1)). Dependences of the intra- and inter-component energies of interactions on concentration were calculated for the 18C6(− 0.6) model, and they are presented in Fig. 9. The energies of valence angle, dihedral angle, and bond distortions are positive (see Fig. 9a), since the zero energy levels correspond to the optimum values of angles and lengths. The large distortions are expectable for the cyclic crown ether molecule. The energies in the plots were fitted by straight lines. Statistical uncertainties are large at small concentrations of 18С6. All slope coefficients are small in the investigated concentration range, indicating no significant intramolecular conformational changes caused by the interaction of 18C6 molecules with one another or by changing their hydration state. Fig. 9b shows the contributions of inter-molecular interactions to the total energy of the system. Coulombic and van der Waals interactions of water molecules and inter-component interactions are specially highlighted. The slopes of the fitting lines again demonstrate weak dependences of the energies on the concentration of the solutions. Except the van der Waals interactions of water molecules, all other interactions are negative and they provide the system stability. Thus, these suggest that there are no compensatory effects, which could explain the weak dependence of hydration enthalpy. This kind of dependence is also the property of each contribution to the total energy.

4. Conclusions ΔH ðmÞ ¼ H sol ðmÞ−Hw −H g ;

ð6Þ

where Hsol, Hw, and Hg are the enthalpies of a solution, pure water, and 18С6 molecules in a gas phase, respectively. The enthalpy of ether in the gas phase was calculated separately. To this purpose several ether molecules were placed in a large simulation box and a calculation was carried out in NVT-ensemble. Results of calculations of hydration enthalpies are presented in Fig. 8. A standard error of mean is large especially at small 18С6 concentrations, even the enthalpy is measured during several nanoseconds. The experimental enthalpy of hydration at infinite dilution (− 154.7 kJ mol−1) [30] and the dependence of

The MD simulations were performed for flexible 18С6 and water molecules. Two models with different atomic charge distributions were investigated. Both models adequately reproduce the properties of actual aqueous solutions of 18-crown-6. The calculated self-diffusion coefficients of the molecules in solutions correspond to the experimental data. Densities, calculated for the 18C6(− 0.4) model, are closer to experimental values, but the 18C6(− 0.6) model better reproduces experimental enthalpies of hydration. Both the force fields, where all molecules are flexible, may be recommended for simulations of multicomponent systems containing water and crown ether molecules. The 50

a

100

-50

80

E, kJ mol-1

90

E, kJ mol-1

b

0

valence angle bond dihedral angle Linear fit

70

-100

Van der Waals water-water crown-water Coulombic water-water crown-water Linear fit

-150 -200 -250

60

-300 50

-350 0

2

4

6

8

Nc

10

12

0

2

4

6

8

10

12

Nc

Fig. 9. Energies of interactions in aqueous solutions of crown ether, 18С6(−0.6): (a) the energies of intramolecular interactions corresponding to distortions of valence, dihedral angles and bond lengths; (b) the energies of van der Waals and Coulombic intermolecular interactions. Bars indicate the standard deviations. Nc is the number of 18С6 molecules in the simulation box.

Y.G. Bushuev et al. / Journal of Molecular Liquids 224 (2016) 825–831

selection depends on a specific task and polarization properties of solutes and co-solvents. The analysis of the results demonstrates that the structure of water is only slightly distorted by the presence of dissolved crown ether molecules. Hydrophilic moieties of 18С6 molecules coordinate in the first hydration shell in average 5.5–6.1 water molecules depending of the used model. The diffusion coefficient of water decreases with increase in concentration of 18С6. Distributions of residence lifetime of water from the first hydration shell of an 18С6 molecule demonstrate a presence of long-lived states. This means that the water and ether molecules can move in solutions together for a long time. The hydrophilic moiety of an ether molecule can strongly bind water. Some water molecules from the first hydration shells of the hydrophobic moieties move together with 18C6 molecules but their lifetimes are several times less than in the case of water from the first hydration shells of the hydrophilic moieties. Disclosure statement No potential conflict of interest was reported by the authors. Funding This work was supported by the Ministry of Education and Science of the Russian Federation under project number 2293. Acknowledgements This work was done at Institute of Thermodynamics and Kinetics of Chemical Processes of Ivanovo State University of Chemistry and Technology. The reported study was supported by the Supercomputing Centre of Lomonosov Moscow State University [project number 1756]. DL_FIELD is a program package written by C. W. Yong and has been obtained from STFC's Daresbury Laboratory via the website http://www.ccp5.ac.uk/DL_FIELD. Appendix A. Supplementary data Supplementary data to this article can be found online at http://dx.doi. org/10.1016/j.molliq.2016.10.062. References [1] G.W. Gokel, W.M. Leevy, M.E. Weber, Crown ethers: sensors for ions and molecular scaffolds for materials and biological models, Chem. Rev. 104 (2004) 2723–2750, http://dx.doi.org/10.1021/cr020080k. [2] V.P. Solov’ev, N.N. Strakhova, O.A. Raevsky, V. Rüdiger, H.-J. Schneider, Solvent effects on crown ether complexations 1, J. Organomet. Chem. 61 (1996) 5221–5226, http://dx.doi.org/10.1021/jo952250h. [3] T. Kowall, A. Geiger, Molecular dynamics simulation study of 18-crown-6 in aqueous solution. 1. Structure and dynamics of the hydration shell, J. Phys. Chem. (1994) 6216–6224, http://dx.doi.org/10.1021/j100075a026. [4] G. Benay, G. Wipff, Ammonium recognition by 18-crown-6 in different solutions and at an aqueous interface: a simulation study, J. Phys. Chem. B 118 (2014) 13913–13929. [5] S. Krongsuk, T. Kerdcharoen, S. Hannongbua, The hydration structure of 18-crown6/K+ complex as studied by Monte Carlo simulation using ab initio fitted potential, J. Mol. Graph. Model. 25 (2006) 55–60, http://dx.doi.org/10.1016/j.jmgm.2005.11. 002. [6] H. Kim, Monte Carlo simulation study of solvent effect on selectivity of 18-crown-6 to between La3+ and Nd3+ Ion, Bull. Kor. Chem. Soc. 24 (2003) 751–756.

831

[7] W. Zielenkiewicz, O.V. Kulikov, I. Kulis-Cwikla, Excess enthalpies and apparent molar volumes of aqueous solutions of crown ethers and cryptand(222) at 25°C, J. Solut. Chem. 22 (1993) 963–973. [8] G. Benay, G. Wipff, Liquid–liquid extraction of alkali cations by 18-crown-6: complexation and interface crossing studied by MD and PMF simulations, New J. Chem. 40 (2016) 2102–2114, http://dx.doi.org/10.1039/C5NJ02609A. [9] T.R. Usacheva, V.A. Sharnin, I.V. Chernov, E. Matteoli, Calorimetric investigation of the complex formation reaction of 18-crown-6 ether with d,l-alanine in water– ethanol mixtures, J. Therm. Anal. Calorim. 112 (2013) 983–989, http://dx.doi.org/ 10.1007/s10973-012-2625-7. [10] M. Olschewski, S. Knop, J. Seehusen, J. Lindner, P. Vöhringer, Ultrafast internal dynamics of flexible hydrogen-bonded supramolecular complexes, J. Phys. Chem. A 115 (2011) 1210–1221, http://dx.doi.org/10.1021/jp110729d. [11] A. Albert, D. Mootz, Formation and crystal structures of the hydrates of 18-crown-6 [1], Z. Naturforsch., B: J. Chem. Sci. 52 (1997) 615–619. [12] J.D. Dunitz, P. Seiler, 1,4,7,10,13,16-Hexaoxacyclooctadecane, Acta Crystallogr. Sect. B. 8924 (1974) 2739–2741. [13] L. Troxler, G. Wipff, Conformation and dynamics of 18-crown-6, cryptand 222, and their cation complexes in acetonitrile studied by molecular dynamics simulations, J. Am. Chem. Soc. 116 (1994) 1468–1480, http://dx.doi.org/10.1021/ja00083a036. [14] K. Fukuhara, K. Ikeda, H. Matsuura, Normal coordinate analysis and force field of 18crown-6, J. Mol. Struct. 224 (1990) 203–224. [15] A. Chaumont, R. Schurhammer, P. Vayssière, G. Wipff, Simulations of the dynamics of 18-crown-6 and its complexes: from the gas phase to aqueous interfaces with SC-CO2 and a room-temperature ionic liquid, in: K. Gloe (Ed.) Macrocycl. Chem. Curr. Trends Futur. Perspect., Springer Netherlands, Dordrecht 2005, pp. 327–348, http://dx.doi.org/10.1007/1-4020-3687-6_21. [16] L.X. Dang, P.A. Kollman, Free energy of association of the 18-crown-6:K+ complex in water: a molecular dynamics simulation, J. Am. Chem. Soc. 112 (1990) 5716–5720, http://dx.doi.org/10.1021/ja00171a006. [17] L.X. Dang, Free energies for association of Cs+to 18-crown-6 in water. A molecular dynamics study including counter ions, Chem. Phys. Lett. 227 (1994) 211–214, http://dx.doi.org/10.1016/0009-2614(94)00810-8. [18] S. Krongsuk, T. Kerdcharoen, S. Hannongbua, How many water molecules in the hydration shell of 18-crown-6? Monte Carlo simulations based on ab initio-derived potential energy surface, J. Phys. Chem. B 107 (2003) 4175–4181, http://dx.doi. org/10.1021/jp026762t. [19] T. Kowall, A. Geiger, Molecular dynamics simulation study of 18-crown-6 in aqueous solution, 2. Free energy profile for the association 18C6⋯K+ in water, J. Phys. Chem. 99 (1995) 5240–5246. [20] Y. Sun, P.A. Kollman, Conformational sampling and ensemble generation by molecular dynamics simulations: 18-crown-6 as a test case, J. Comput. Chem. 13 (1992) 33–40, http://dx.doi.org/10.1002/jcc.540130105. [21] M. Billeter, A. Howard, A new technique to calculate low-energy conformations of cyclic molecules utilizing the ellipsoid algorithm and molecular dynamics: application to 18-crown-6, J. Am. Chem. Soc. 110 (1988) 8385–8391, http://dx.doi.org/10. 1021/ja00233a016. [22] A. Howard, U. Singh, Many-body potential for molecular interactions, J. Am. Chem. Soc. 110 (1988) 6984–6991, http://dx.doi.org/10.1021/ja00229a009. [23] W. Smith, T.R. Forester, DL_POLY_2.0: a general-purpose parallel molecular dynamics simulation package, J. Mol. Graph. 14 (1996) 136–141, http://dx.doi.org/10. 1016/S0263-7855(96)00043-4. [24] W.D. Cornell, P. Cieplak, C.I. Bayly, I.R. Gould, K.M. Merz, D.M. Ferguson, D.C. Spellmeyer, T. Fox, J.W. Caldwell, P.A. Kollman, A second generation force field for the simulation of proteins, nucleic acids, and organic molecules, J. Am. Chem. Soc. 117 (1995) 5179–5197, http://dx.doi.org/10.1021/ja00124a002. [25] Y. Wu, H.L. Tepper, G.A. Voth, Flexible simple point-charge water model with improved liquid-state properties, J. Chem. Phys. 124 (2006) 24503, http://dx.doi. org/10.1063/1.2136877. [26] O. Mayzel, Y. Cohen, Diffusion coefficients of macrocyclic complexes using the PGSE NMR technique: determination of association constants, J. Chem. Soc. Chem. Commun. 1901 (1994) http://dx.doi.org/10.1039/c39940001901. [27] K.J. Patil, S.R. Heil, M. Holz, M. Zeidler, Self-diffusion coefficient and apparent molar volume studies of crown ethers in aqueous (D2O) and CDCl3 solutions, Ber. Bunsenges. Phys. Chem. 101 (1997) 91–95, http://dx.doi.org/10.1002/bbpc. 19971010112. [28] W.S. Price, H. Ide, Y. Arata, Self-diffusion of supercooled water to 238 K using PGSE NMR diffusion measurements, J. Phys. Chem. A 103 (1999) 448–450, http://dx.doi. org/10.1021/jp9839044. [29] D. Dagade, R. Pawar, K. Patil, Viscosity behavior of 18-crown-6 in aqueous and carbon tetrachloride solutions at different temperatures and at ambient pressure, J. Chem. Eng. Data 49 (2004) 341–346, http://dx.doi.org/10.1021/je034188o. [30] L.-E. Briggner, I. Wadsö, Some thermodynamic properties of crown ethers in aqueous solution, J. Chem. Thermodyn. 22 (1990) 143–148, http://dx.doi.org/10. 1016/0021-9614(90)90077-4.