- Email: [email protected]
Statistical geometry characterization of local structure of TMAO, TBA and urea aqueous solutions
Journal of Molecular Liquids 286 (2019) 110870
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
Statistical geometry characterization of local structure of TMAO, TBA and urea aqueous solutions☆ E.D. Kadtsyn, A.V. Anikeenko, N.N. Medvedev ⁎ Voevodsky Institute of Chemical Kinetics and Combustion, SB RAS, 630090 Novosibirsk, Russia Novosibirsk State University, 630090 Novosibirsk, Russia
a r t i c l e
i n f o
Article history: Received 14 November 2018 Received in revised form 10 April 2019 Accepted 27 April 2019 Available online 6 May 2019
a b s t r a c t Recently, we have revealed that trimethylamine-N-oxide (TMAO) molecules in aqueous solutions are distributed generally like random spheres in space (A.V. Anikeenko et al. J. Mol. Liq. 245 (2017) 35–41). We have shown that the variance of the distribution of Voronoi region volumes for the TMAO molecules and the hard random spheres are identical at the corresponding concentrations up to the mole fraction x ~ 0.1 (or packing fraction η ~ 0.2). In the current work, we study clusters in TMAO solutions and calculate the weighted mean cluster size and some topological cluster characteristics. For all these parameters, a great matching of TMAO and random spheres is also observed. The same analysis for tert-butyl alcohol (TBA) and urea aqueous solutions is carried out. These molecules have a tendency to cluster even at very low concentrations. However, basic geometrical features of the clusters are the same both for solutions and for the system of random spheres, indicating that geometrical laws governing the allocation of non-overlapping spheres in space should be taken into account when describing the structure of solutions even at rather low concentrations. © 2019 Elsevier B.V. All rights reserved.
1. Introduction The impenetrability of molecules (atoms) is usually considered as a geometrical factor that defines the size and shape of a molecule, whereas physicochemical properties of a substance are regarded to be related to the specifics of intermolecular interactions, which are divided into dispersive, electrostatic, etc. [1]. Impenetrability, however, plays a significant role in the condensed phase, determining the mutual arrangement of atoms and molecules and imposing restrictions on their movement. In this case, the details of intermolecular interactions become less important. Thus, noble gases and many metals, such as gold, copper, lead, iron, form dense crystalline structures with high symmetry, in particular, bcc or fcc, regardless of the characteristics of their outer electronic shells. The same structure is obtained in the system of hard spheres, which have no other interactions except impenetrability. A similar situation is observed for liquids. All atoms mentioned, as well as compact molecules (for example, methane), exhibit similar radial distribution functions in a liquid state. The presence of the excluded volume inevitably leads to the fact that the location of such particles is
☆ The article belongs to the special issue on The problems of solvation, complex formation and crystallization challenging the development of smart materials. ⁎ Corresponding author at: Voevodsky Institute of Chemical Kinetics and Combustion, SB RAS, 630090 Novosibirsk, Russia E-mail address: [email protected] (N.N. Medvedev).
https://doi.org/10.1016/j.molliq.2019.04.147 0167-7322/© 2019 Elsevier B.V. All rights reserved.
governed by the laws of geometry. In particular, it was proved that the densest packing of identical hard spheres is the crystal structure of fcc or hcp [2], and in the case of non-crystalline systems (dense monoatomic liquids, glasses and dense packings of spheres), the socalled “polytetrahedral packing principle” works, where tetrahedral configurations of four atoms form linear and branched clusters [3,4]. It should be emphasized that the above structural patterns are not prerogative of hard spheres. They also hold for reasonably soft particles and compact molecules with only relative sphericity. An interesting example is the substances which have the plastic crystal phase. These are the compact molecules like neopentane, carbon tetrachloride, 2,3dimethylbutane, cyclohexane [5,6]. At low temperatures, each of them has its own specific crystal structure. However, when a liquid is cooled, a plastic crystal first appears: translational diffusion ceases, and a crystalline order corresponding to a bcc or fcc structure arises. The analysis of such structures shows that the molecules are located in the lattice sites, but are oriented differently, and rotate slowly relative to each other [6,7]. In the liquid phase, such molecules behave like atoms in simple liquids. In [8], it was shown that the radial distribution functions calculated for the centers of the molecules of cyclohexane and 2,3dimethylbutane coincide with those for the Lennard-Jones system. It was also shown that the similarity is also manifested in the characteristics of the “three-dimensional” structure, studied with the help of Delaunay simplices [9]. Thus, the high symmetry of plastic crystals, where the molecules actually “hang” freely in the lattice sites, and the
2
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
universal structure of the fluids of all spherically symmetric atoms and compact molecules, is a manifestation of the geometric laws governing the placement of impenetrable spherical particles in space. In contrast to condensed systems, little attention is paid to the study of the spatial distribution of particles in systems with low density. In pure form, they are gases or fluids, whose structure is usually not discussed. However, there are similar systems that are of great interest. These are solutions whose structure and properties are also studied at low concentrations of dissolved molecules. However, the chemistry of solutions focuses on the specific interactions and the impenetrability of the dissolved molecules at relatively low concentrations is usually not taken into account. As a rule, to describe the properties of solutions, the models are used in which the excluded volume of molecules is not taken into account, for example, the Debye-Hückel model [1]. It leads, in particular, to the fact that the occurrence of associates of dissolved molecules is explained solely by their specific interactions with each other or with a solvent. Such explanations neglect that even at low concentrations, the excluded volume may influence this process. In our previous work [10], we calculated the normalized variances of the Voronoi region volumes for randomly distributed equal hard spheres at different packing fractions (degrees of space filling). Zero packing fraction corresponds to a system of random points. For that zero packing fraction, there is a wide distribution of the Voronoi volumes, because there are large fluctuations in local environments. The dispersion of volumes decreases with an increase in the packing fraction of the system (with the radius of the spheres also increasing in a given model box). The system becomes more homogeneous: instead of pairs and groups of close random points, there are now clusters of spheres with the centers separated by the distance determined by the radius of that spheres [10]. When comparing the results for random spheres with the data obtained for aqueous solutions of TMAO, we found that over a wide range of solution concentrations (up to 0.2, in terms of the packing fraction η for the solute molecules) dispersion of volumes of Voronoi regions for TMAO molecules behave the same way as for random spheres at corresponding concentrations. This fact can be explained by the TMAO molecules being distributed in the solution independently, thus not exhibiting the specific interaction with each other. The question of why the TMAO molecules behave independently in solution remains open. These molecules have a large dipole moment, form hydrogen bonds, and affect water by their hydrophobic part, thus being able to interact both directly with each other and indirectly through water [11,12]. Apparently, the interactions of different nature somehow compensate each other. This question is worth being studied separately, especially, as we show that the similarity of the spatial distribution of TMAO in solution and random hard spheres is manifested not only in the dispersion of the volumes of the Voronoi regions but also in the structure of the emerging clusters. In this paper, in addition to TMAO solutions, we also study the solutions of tert-butanol and urea. These molecules tend to associate in aqueous solutions, and they have a high first peak of pair correlation function even at very low concentrations. At the same time, the mechanisms of association in TBA and urea are different. In the case of TBA, it is generally accepted that it is a hydrophobic interaction [13–15]. It is believed that water “expels” TBA molecules and they create relatively compact associates. In the case of urea, the situation is less clear. The urea molecule does not have a hydrophobic group, but hydrogen bonding is possible. The urea molecules are thought to be included in the structure of the hydrogen bonding network [16–19].
1 bar. A stochastic velocity rescaling thermostat [21] and ParinelloRahman barostat [22] were used. Integration step was 2 fs. The production runs were 100 ns length. Every 10th picosecond frame was used for analysis (104 independent configurations for each model). Parameters of force fields were taken from [23] for TMAO, from [24] for TBA, and for urea from [25]. TIP4P water was used [26]. The structure of these molecules is shown in Fig. 1. We prepared twenty models of TMAO and TBA aqueous solutions on the molar percentage interval from 0.5% to 10% with the step of 0.5%. Each model contained 100 solute molecules in the simulation box. For the lowest concentration, the box had an edge length of 8.5 nm and contained 19,900 water molecules, see Table 1. The more concentrated solutions were obtained by removing the required number of water molecules and decreasing the box size [10]. For urea, we also got twenty models but on the interval from 1% to 20% with the step of 1%. For the concentrations from 1 to 10% we used 100 urea molecules. However, for the interval 11–20% we used 200 molecules to keep the box size large enough at the higher concentrations. The larger mole fraction of urea was required because the urea molecule is smaller. So further we used the packing fraction as a measure of concentration. The relation between molar fraction, molar concentration, packing fraction and the corresponding number of water molecules in the models are shown in Table 1. The comparison of the calculated and experimental densities of the solutions is shown in Fig. 2. One can see quite a good agreement for all models for used intervals of concentration.
2. Methods
¼ NπD 6V box , where D. is the mean diameter of the spheres. It is easy to see
2.2. Systems of random spheres All models of random spheres were created by random sequential addition of hard spheres into empty simulation box with periodic boundary conditions. For cluster analysis, we have used only systems of monodisperse spheres. The diameter of these spheres was chosen to be 0.53 nm, and the packing fraction was determined by the variation of the size of the box. (Remember that packing fraction is the ratio of the 3
volume occupied by spheres to the volume of the system: η ¼ NπD 6V box , where Vbox is the volume of the model box, N is the number of spheres, and D is their diameter). In this work, we have used N = 100. For each packing fraction, 104 independent configurations for the similarity with solutions models in the statistical analysis were obtained. The systems of monodisperse spheres were obtained for packing fraction from η = 0.01 to 0.21 with the step of 0.01. We have also obtained the models of polydisperse spheres with Gaussian diameter distribution. These models were used for determining the effective radius of solute molecules, see below. We have varied both mean diameter D and its variance. Fig. 3 shows the pair correlation functions g(r) for the system of monodisperse and polydisperse spheres with the same diameter. One can see that vertical left side which is intrinsic for monodisperse spheres became smoother for polydisperse spheres. Note that the middle of this front corresponds well to the mean diameter of spheres at least for low concentrations. Such observation confirms reasonability of our suggestion used in [10], where values of effective diameters of TMAO and TBA were defined as the position of the middle of the left front of the first peak of g(r). Note that if σ is small enough (a few percent of diameter) then the packing fraction for the system of such polydisperse spheres systems can be calculated by using the formula for the monodisperse spheres η 3
2.1. Solutions. Molecular dynamics simulations The molecular dynamics simulation protocol was the same as in [10], using Gromacs 5.0 package [20] at temperature 300 K and pressure
that contributions from the spheres of the sizes being greater or less than the mean value compensates each other with the precision of the second order of smallness, which can be neglected. To estimate the volume of solute molecules in solution we represent them as spheres, see below. On the other hand, as is well known in
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
3
Fig. 1. Structure and van der Waals representation of TMAO, TBA, and urea molecules (from left to right).
chemistry and biology, the volume of a molecule can be calculated using van der Waals surface of atoms, or the molecular (Connolly) surface. Our choice is because we use a system of hard spheres as a reference. 2.3. The effective diameter of molecules Real molecules have a hard core, but they are not rigid themselves. To determine the effective diameter of the molecules, we use the models of random polydisperse spheres. The values of mean diameter D and its variance are chosen in such a way that the left front of the first peak of g(r) for the spheres system could match the left front of g (r) of the solute molecules studied. Note that in such comparison, the same number of molecules and packing fraction of the systems is used. Fig. 4 shows the g(r) for TMAO solution and random polydisperse spheres at different packing fractions. In both cases, the mean diameter of the spheres is 0.53 nm. The dispersion at both concentrations is also equal (σ = 0.035 nm) because the position and the steepness of the peak do not depend on the concentration. Note that in the previous work [10] we obtained the same value of molecule diameter. However, in the current work, we use the value of D, which was obtained from the fitting of the whole left front. We do not pay attention now to the area beyond the first peak. The curves of g(r) have different behavior there. For the random spheres at
the used packing fraction, far oscillations still do not appear, while for the solutions clear oscillations are observed. They are a result of “matrix influence”. In other words, they are determined by the structure of the surrounding water [30]. When determining the effective diameters of TBA and urea molecules, the additional question arises related to the fact that they have a high first peak of g(r). However, even in this case, the left front reflecting the variations of closest molecular contacts can be well described by using our approach with polydisperse spheres. The height of the first peak is caused by the larger number of the nearest solute molecules. The increase in the number of contacts is caused by the structure of solution, interactions with solvent and is not connected with direct intermolecular contacts and with the “hardness” of the molecules. In other words, the left front of the first peak can be considered without regard to its height. Therefore, we simply reduced the height of the first peak of the solution to the value of the system of random spheres by applying the corresponding normalizing factor. Fig. 5 presents the fitting of the left front of TBA and urea. We have obtained the values D = 0.52 nm for TBA and D = 0.39 nm for urea. For a higher concentration, the same good description is obtained, giving the same values of the effective diameter of the molecules. For these solutions, as well as for TMAO, the position and steepness of the left front remain almost unchanged in the concentration range used.
Table 1 3
Correspondence between molar fraction (x), molar concentration (C) and packing fraction (η), for TMAO, TBA and urea aqueous solutions. Packing fraction was calculated as η ¼ NπD 6V box , where Vbox is the volume of the model box, N – the number of spheres, and D is the effective diameter of the molecule: DTMAO = 0.53 nm, DTBA = 0.52 nm, Durea = 0.39 nm, see text. TMAO x, % C, mol/l pf,% Water molecules x, % C, mol/l pf, % Water molecules
0.5 0.273 1.31 19,900 5.5 2.613 12.25 1718
1.0 0.537 2.58 9900 6.0 2.815 13.16 1566
1.5 0.794 3.8 6566 6.5 3.015 14.04 1438
2.0 1.044 4.99 4900 7.0 3.205 14.89 1328
2.5 1.287 6.13 3900 7.5 3.393 15.71 1233
3.0 1.523 7.24 3233 8.0 3.575 16.49 1150
3.5 1.752 8.31 2757 8.5 3.756 17.28 1076
4.0 1.976 9.34 2400 9.0 3.927 18.03 1011
4.5 2.194 10.34 2122 9.5 4.083 18.77 952
5.0 2.407 11.31 1900 10.0 4.248 19.46 900
0.5 0.272 1.21 19,900 5.5 2.537 11.24 1718
1.0 0.534 2.37 9900 6.0 2.726 12.08 1566
1.5 0.788 3.49 6566 6.5 2.908 12.89 1438
2.0 1.033 4.58 4900 7.0 3.084 13.67 1328
2.5 1.270 5.63 3900 7.5 3.254 14.42 1233
3.0 1.500 6.65 3233 8.0 3.417 15.14 1150
3.5 1.721 7.63 2757 8.5 3.580 15.87 1076
4.0 1.935 8.58 2400 9.0 3.736 16.56 1011
4.5 2.143 9.5 2122 9.5 3.889 17.24 952
5.0 2.343 10.38 1900 10.0 4.031 17.87 900
7 3.52 6.58 1328 17 7.55 14.12 976
8 3.97 7.42 1150 18 7.9 14.77 911
9 4.41 8.24 1011 19 8.25 15.42 852
10 4.84 9.04 900 20 8.58 16.04 800
TBA x, % C, mol/l pf, % Water molecules x, % C, mol/l pf, % Water molecules Urea x, % C, mol/l pf, % Water molecules x, % C, mol/l pf, % Water molecules
1 0.55 1.02 9900 11 5.25 9.82 1618
2 1.08 2.01 4900 12 5.66 10.58 1466
3 1.59 2.98 3233 13 6.06 11.33 1338
4 2.09 3.91 2400 14 6.44 12.05 1228
5 2.58 4.83 1900 15 6.82 12.75 1133
6 3.06 5.72 1566 16 7.19 13.44 1050
4
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
system and the appearing clusters are shown in Fig. 7. The particles intersecting by their permeable shells are combined into the cluster.
3. Results
Fig. 2. The densities of aqueous solutions of urea, TMAO and TBA at T = 300 K and P = 1 bar vs. concentration. Solid lines refer to experimental data [27–29], large symbols show the simulation values.
Note that the value obtained for the TBA diameter slightly differs from the one for TMAO. In the [10] we used the same value D = 0.53 nm for both solutes for simplicity. However, such a distinction is not significant for our analysis. A small shoulder on the left slope of the TBA curve is due to the hydrogen bond formation between some TBA molecules. Fig. 5а shows that the fraction of such molecules is small, about 2–3%. This fact was also proven by our special analysis. Therefore, we do not pay any particular attention to these molecules.
2.4. Selection of the critical distances for cluster determination The choice of critical distances between the molecules to include them into one cluster is always arbitrary. In this work, we associate this distance with the molecular diameter D and use three values: 1.25 D, 1.30 D and 1.35 D. On the Fig. 6 all functions g(r) for our systems are represented in the units of D. One can see that critical distances used in our study go through the main part of the first peak, see vertical dotted lines. Recall that from the mathematical standpoint our systems are the systems of “hard spheres with permeable shells”, which used for studying the continuous percolation years ago [31–33]. An illustration of such
Fig. 8 demonstrates variances of Voronoi volumes for solute molecules in our solutions and in the hard spheres system at different packing fractions. Recall that Voronoi regions were calculated for the system of centers of the solute molecules, see details in [10]. In this case, we describe a general motif of the spatial distribution of the molecules. In current work, we represent the data for the new models and show the results for urea models, containing 100 and 200 solute molecules, see Table 1. The most surprising result is that curves for TMAO solutions and random spheres match almost exactly. TBA and urea curves show different behavior. We have related the increase of variance for TBA with the arising of inhomogeneity in its solutions. For urea, however, the variance decreases monotonically. It means that urea solutions become more homogeneous but to a lesser extent than in the case of TMAO or random spheres. In the following, we present the results of a cluster analysis which allow studying our systems from another side, namely from the standpoint of local structure. First, we calculated the simple mean cluster size, which shows the average number of particles in the cluster. It is defined as 〈i〉 = ∑i=1Nini/∑i=1Nni, where ni is the number of clusters containing i particles and the sum is taken over all existing clusters, which can include also single molecules (i = 1) or cluster formed by all solute molecules (i = N). It can also be calculated using standard gmx clustsize utility from the Gromacs package. This value is shown in Fig. 9 for our systems in dependence on packing fraction for different critical distance values. It should be emphasized that TMAO molecules in the solutions behave again like random spheres: black and red curves match for all critical distances. The clusters of TBA and urea differ from the TMAO ones and the random spheres. The larger size of these clusters can be explained by the higher number of the closest solute molecules in solutions. Indeed, both TBA and urea have high first peaks of g(r). This simple parameter, however, has a low sensibility for the cluster analysis. In the percolation theory, another mean value is commonly used: the second moment of the cluster size distribution normalized by the first moment: S = ∑i=1Ni2ni/∑i=1Nini [34–36]. In fact, it is a weighted (by the number of particles) the mean size of the cluster. Indeed, this formula can be rewritten as S = ∑i=1Nfii, where fi is the weight of the cluster with size i. It is defined as f i ¼ inNi , where ini is the number of particles in the all clusters of a given size i, and in the
Fig. 3. Function g(r) for systems of monodisperse random spheres of unit diameter (black) and for polydisperse random spheres with the same mean diameter and different dispersions: 0.05 (blue) and 0.03 (green). Models packing fraction η = 1% (a), and η = 15% (b).
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
5
Fig. 4. Comparison of g(r) for TMAO molecules in solution (red) and random polydisperse spheres system (black) with D = 0.53 nm, σ = 0.025 nm. TMAO concentration x = 0.05, spheres packing fraction η = 0.01 (a). TMAO x = 0.075, η = 0.15 (b), see Table 1.
denominator there is a full number of particles in the system, which was previously written as the sum of the particles over all clusters: ∑i=1Nini = N. In such definition of the weighted mean size, the more massive size clusters (which are always not high in number) are not suppressed by small clusters (which are always high in number). Note that such a parameter is used in chemistry in the study of complexation in the solutions, where the weight factor fi is called form fraction of i. [37]. In the percolation theory, one deals with infinitely large systems. In this case, the weighted mean size is calculated only for finite clusters. It is known, its value diverges according to a power law near the percolation threshold. However, for relatively small models, the percolation threshold is not well defined (large fluctuations appear at its determination) [36]. For this reason, we do not calculate the percolation cluster in our models and use the cluster size only for the comparison of the different systems. Fig. 10 shows the weighted mean cluster size normed by the full number of particles in the model. (Our models contain different number of urea molecules N. For comparing all models together, we use the normed cluster size S/N). The curves aim to the limit value of the unit with the increasing concentration. It means that almost all the solute molecules are included in one cluster at the high concentration in our solutions, namely at the packing fraction η ≈ 0.2. Therefore, our systems are beyond the percolation threshold at the given concentrations. In the work [33], percolation in the system of hard spheres with permeable shells was thoroughly studied with Monte-Carlo technique. The models with different shell thickness and packing fraction were investigated. For the shells corresponding to our critical distances (1.25, 1.30,
1.35), it was found that the percolation occurs at the packing fractions η are equal to 0.175, 0.15, 0.13 respectively. In Fig. 10, the vertical dotted lines mark these values. One can see that in all cases it corresponds approximately to the value of the normed mean size S/N ≈ 0.3. It seems reasonable to use this result for defining the critical concentration for our solutions. For TBA solutions, such a value of mean cluster size corresponds to the packing fraction close to 0.10 (at all critical distances). From Table 1, one can see that this value corresponds to the mole fraction х = 0.048. This value is in agreement with existing data. In the work [38], it was observed the percolation threshold of the TBA molecules at x ≈ 0.05. The authors obtained “islands” of TBA even at the lower mole fraction, while the large spanning cluster emerged above that mole fraction. A huge model of TBA solution was studied in [39]. It was noted that TBA aggregation could be detected at x ≈ 0.03 and a rapid increase in the size of the clusters is observed in the range x ≈ 0.03–0.06. By applying our indirect method, we can give even a more accurate value of the percolation threshold in TBA solutions. For the urea, the behavior of the clusters differs from the behavior for TBA. One can also see a difference in the various critical distances. It means a different cluster nature for urea solutions. Nevertheless, one can say that the percolation threshold for the urea solution is approximately at η = 0.12, corresponding to its molar fraction х = 0.14, that is much higher than for TBA. Unfortunately, we do not know any works that would study the percolation in urea solutions. Look at to the slight break at η ≈ 0.09 on the urea curves. It occurs because of a change in the size of the model. At a higher concentration of urea, we use a twice bigger box and worked with 200 solute
Fig. 5. Comparison of g(r) for TBA at x = 0.005 (a) and urea at x = 0.01 (b) with the random polydisperse spheres system at packing fraction of η = 0.01. The functions for TBA and urea are scaled to get the unit height of the first peak, see text. Fitting parameters are D = 0.52 nm, σ = 0.02 nm for TBA and D = 0.39 nm, σ = 0.017 nm for urea.
6
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
Fig. 6. Functions g(r) for TMAO, TBA, and urea in the units of the effective diameter of the molecule. Black curve corresponds to the random polydisperse spheres system. Red curves correspond to TMAO, blue to TBA, and green to urea. Vertical dashed line shows an effective diameter equal to the unit; black dotted lines show suggested critical distances for cluster selection.
molecules instead of 100 for low concentration (see above). The displacement of curves means that cluster composition in our models slightly differs. Indeed, for our rather small models, the fraction of bigger clusters and the fraction of single molecules in different sized boxes can be different. Fig. 9 does not show such a break at the change of the model size. That can be because of the low structural sensibility of that parameter - the simple mean size of the clusters. Fig. 11 demonstrates the mean degree of the nodes (the average number of the neighbors for a given molecule in the cluster) in dependence on the cluster size. We use here only rather small clusters that involve no more than 10 molecules. At relatively low concentration, there are quite a lot of such clusters in the solutions. As the concentration increases, their number significantly decreases due to the appearance of large clusters. For this reason, we work here only with small packing fractions: 0.04, 0.05 and 0.06 for hard spheres and with values for solutions close to them, see Table 1. Thus, for each solution and for random spheres system, three curves are shown corresponding to the given packing fractions, Fig. 11. In addition, a theoretical curve for tree clusters (cluster without cycles) is shown, which was calculated by the formula hdi ¼ 2−2 i .
Fig. 8. The variance of normed Voronoi region volume distributions for TBA (blue) urea (green) and TMAO (red) solute molecules in solution compared with random spheres (black). A break at η = 0.09 on the urea curve is due to size effect, see text.
Note that the degree of the nodes almost does not depend on the chosen critical distance, compare the figures (a)–(c). However, there is a clear difference between the systems. The highest degree of the nodes is for the TBA solutions (blue curves), the least one is for TMAO solutions and for the random spheres system. We emphasize that the curves for TMAO and for random spheres coincide very well also for these parameters. Note that for clusters formed by two particles, the degree of the node is 1. With the number of particles in a cluster increasing, the degree also increases. For trees, it approaches the theoretical limit equal to 2. However, for TBA solutions, it approaches the value 2.5, indicating the presence of cycles. The greater the mean degree of the nodes, the more cycles in the cluster of a given size. For the simplest case of a cluster of a size 3 there are only linear clusters or triangles. The fraction of triangles among all clusters of size 3 is shown in Fig. 12 as a function of concentration. The smallest amount of triangles is in the TMAO solutions and random spheres. In the TBA solutions, their fraction is more than twice greater. In the urea solutions, the fraction of triangles is smaller than in TBA solutions, suggesting again that the nature of clusters is different for these two systems. Note that this parameter also shows the structure
Fig. 7. Schematic representation of hard spheres with the permeable shell. The black areas define packing fraction η for the system (a). The clusters corresponding to the system (b).
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
7
Fig. 9. The mean size of clusters in solutions and the random sphere system depending on the packing fraction for different critical distances (see text): 1.25D (a), 1.3D (b) and 1.35D (c).
Fig. 10. The weighted mean size of clusters for different critical distances, see caption in Fig. 9. The vertical dotted lines mark the percolation thresholds for the corresponding random hard spheres with permeable shells obtained by Monte-Carlo technique [33], see text.
8
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
Fig. 11. Mean degree of the nodes in the cluster for clusters of different size for critical distances 1.25D (a), 1.3D (b) and 1.35D (c). Three groups of curves belong TMAO, TBA and urea solutions. Three curves in each group corresponding to the concentrations η ≈ 0.03, 0.04 and 0.06. Black curves represent theoretical values for trees, i.e. clusters without cycles.
closeness of TMAO solutions to the random sphere system, see black and red curves in Fig. 11. 4. Conclusions In this work, we have studied clusters of TMAO, TBA and urea molecules in aqueous solutions at different concentrations. In our previous work [10] we revealed that TMAO molecules are distributed in the solutions like the random non-overlapped spheres up to the packing fraction of η = 0.2. Now we have shown that TMAO clusters are also the
Fig. 12. The fraction of triangles among all clusters of size 3 depending on concentration. The results for critical distance 1.3D are shown.
same as clusters in the random spheres system. The question of the reason for such behavior of TMAO molecules remains open. They can interact with each other both directly and through the water due to having the large dipole moment, forming strong hydrogen bonds and having large hydrophobic moiety. The question of how all these interactions can be mutually compensated is worthy of a separate study. We want to emphasize that the similarity with random sphere system can be observed in a rather wide range of concentrations and manifests itself as in the spatial distribution of molecules as well as the features of intermolecular contacts. The distribution of TBA and urea in aqueous solution is different from one of random spheres. It is apparently the consequence of the specific intermolecular interactions that are not compensated like in TMAO solutions. The global distribution of these molecules is more inhomogeneous, and small clusters are larger and have more cycles than in the random spheres at the same concentrations. We estimated the percolation threshold in TBA and urea solutions is lower than in the random system. The reason is that these molecules begin to associate even at very low concentrations. However, the nature of the interactions of TBA and urea molecules is different, leading to a more homogeneous distribution of urea molecules and less number of cycles in the urea clusters. An interesting difference in the structure of these two solutions was noticed in the work [30]. It was shown that the height of the first peak of g(r) for urea decreases with concentration, while increases for TBA. It was suggested that urea clusters contain fewer cycles (are more ramified) than in the TBA case. Our analysis is in agreement with that explanation. However, despite the differences between TBA and urea solutions and their difference from TMAO and random spheres, there are common features of the cluster behavior with concentration for all these systems. The main cluster parameters follow to general geometrical laws for the systems with excluded volume, which is manifesting themselves even at low concentrations. Thus, geometric properties of the
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
arrangement of impermeable molecules should be always taken into account when studying the structure of solutions. Despite the difference in the properties of clusters for TBA and urea from the system of hard spheres, it should be noted a similar behavior of the curves in Figs. 9–12 for all these systems. The similarity is obviously a consequence of the existence of an excluded volume, i.e. mutual impermeability of molecules. This means that when explaining the structure of solutions, it is necessary to take into account the geometrical laws of the spatial arrangement of particles even at rather low concentrations. Acknowledgments This work was supported by Russian Foundation for Basic Research (project No. 18-03-00045) and partly by Russian Ministry of Science and Education under 5-100 Excellence Programme. References [1] P. Atkins, J. de Paula, Physical Chemistry: Thermodynamics, Structure, and Change, W. H. Freeman and Company, New York, 2002. [2] T.C. Hales, A proof of the Kepler conjecture, Ann. Math. (2005) 1065–1185. [3] J.D. Bernal, The Bakerian lecture, 1962. The structure of liquids, Proc. R. Soc. Lond. A Math. Phys. Sci. 280 (1382) (1964) 299–322, https://doi.org/10.1098/rspa.1964. 0147. [4] A.V. Anikeenko, N.N. Medvedev, Polytetrahedral nature of the dense disordered packings of hard spheres, Phys. Rev. Lett. 98 (23) (2007) 235504, https://doi.org/ 10.1103/PhysRevLett.98.235504. [5] C.N.R. Rao, Molecular motion in plastic crystals, Proc. Indian Acad. Sci. Chem. Sci. 94 (1) (1985) 181–199, https://doi.org/10.1007/BF02841265. [6] H. Farman, J.C. Dore, M.C. Bellisent-Funel, D.G. Montague, Structural studies of cyclohexane, C6D12 by neutron diffraction. II the plastic crystal phase, Mol. Phys. 73 (4) (1991) 855–871, https://doi.org/10.1080/00268979100101611. [7] J.H. Strange, Pressure and temperature dependence of molecular motion in organic plastic crystals, J. Chem. Phys. 68 (7) (1978) 3078–3088, https://doi.org/10.1063/1. 436174. [8] A.V. Anikeenko, A.V. Kim, N.N. Medvedev, Molecular dynamics simulation of the structure of c6 alkanes, J. Struct. Chem. 51 (6) (2010) 1090–1096, https://doi.org/ 10.1007/s10947-010-0167-z. [9] A.V. Anikeenko, A.V. Kim, N.N. Medvedev, Delaunay simplexes in liquid cyclohexane, Proceedings of Sixth International Symposium on Voronoi Diagrams in Science and Engineering (ISVD 2009), 23–26 July, Technical University of Denmark, Denmark 2009, pp. 271–277, https://doi.org/10.1109/ISVD.2009.10 , published by IEEE Computer Society. [10] A.V. Anikeenko, E.D. Kadtsyn, N.N. Medvedev, Statistical geometry characterization of global structure of TMAO and TBA aqueous solutions, J. Mol. Liq. 245 (2017) 35–41, https://doi.org/10.1016/j.molliq.2017.06.001. [11] L. Larini, J.E. Shea, Double resolution model for studying TMAO/water effective interactions, J. Phys. Chem. B 117 (2013) 13268–13277, https://doi.org/10.1021/ jp403635g. [12] M.V. Fedotova, S.E. Kruchinina, G.N. Chuev, Hydration structure of osmolyte TMAO: concentration/pressure-induced response, New J. Chem. 41 (2017) 1219, https:// doi.org/10.1039/c6nj03296f. [13] H.S. Frank, M.W. Evans, Free volume and entropy in condensed systems III. Entropy in binary liquid mixtures; partial molal entropy in dilute solutions; structure and thermodynamics in aqueous electrolytes, J. Chem. Phys. 13 (11) (1945) 507–532, https://doi.org/10.1063/1.1723985. [14] G. Onori, A. Santucci, Dynamical and structural properties of water/alcohol mixtures, J. Mol. Liq. 69 (1996) 161–181, https://doi.org/10.1016/S0167-7322(96)90012-4. [15] A. Di Michele, M. Freda, G. Onori, M. Paolantoni, A. Santucci, P. Sassi, Modulation of hydrophobic effect by cosolutes, J. Phys. Chem. B 110 (42) (2006) 21077–21085, https://doi.org/10.1021/jp068055w. [16] Y.L.A. Rezus, H.J. Bakker, Destabilization of the hydrogen-bond structure of water by the osmolyte trimethylamine N-oxide, J. Phys. Chem. B 113 (13) (2009) 4038–4044, https://doi.org/10.1021/jp805458p.
9
[17] J.K. Carr, L.E. Buchanan, J.R. Schmidt, M.T. Zanni, J.L. Skinner, Structure and dynamics of urea/water mixtures investigated by vibrational spectroscopy and molecular dynamics simulation, J. Phys. Chem. B 117 (42) (2013) 13291–13300, https://doi.org/ 10.1021/jp4037217. [18] Y.L.A. Rezus, H.J. Bakker, Effect of urea on the structural dynamics of water, PNAS 103 (49) (2006) 18417–18420, https://doi.org/10.1073/pnas.0606538103. [19] D. Bandyopadhyay, S. Mohan, S.K. Ghosh, N. Choudhury, Molecular dynamics simulation of aqueous urea solution: is urea a structure breaker? J. Phys. Chem. B 118 (40) (2014) 11757–11768, https://doi.org/10.1021/jp505147u (2014). [20] S. Pall, M.J. Abraham, C. Kutzner, B. Hess, E. Lindahl, Tackling exascale software challenges in molecular dynamics simulations with GROMACS, International Conference on Exascale Applications and Software, 1, Springer, Cham 2014, April, pp. 3–27, https://doi.org/10.1007/978-3-319-15976-8_1. [21] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (2007), 014101. https://doi.org/10.1063/1.2408420. [22] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182, https://doi.org/10.1063/1. 328693. [23] C. Hölzl, P. Kibies, S. Imoto, R. Frach, S. Suladze, R. Winter, ... S.M. Kast, Design principles for high–pressure force fields: aqueous TMAO solutions from ambient to kilobar pressures, J. Chem. Phys. 144 (14) (2016) 144104, https://doi.org/10.1063/1. 4944991. [24] C. Caleman, P.J. van Maaren, M. Hong, J.S. Hub, L.T. Costa, D. van der Spoel, Force field benchmark of organic liquids: density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant, J. Chem. Theory Comput. 8 (1) (2011) 61–74, https://doi.org/10. 1021/ct200731v. [25] S. Weerasinghe, P.E. Smith, A Kirkwood− buff derived force field for mixtures of urea and water, J. Phys. Chem. B 107 (16) (2003) 3891–3898, https://doi.org/10. 1021/jp022049s. [26] J.L.F. Abascal, C. Vega, A general purpose model for the condensed phases of water: TIP4P/2005, J. Chem. Phys. 123 (2005), 234505. https://doi.org/10.1063/1.2121687. [27] D.M. Makarov, G.I. Egorov, A.M. Kolker, Density and volumetric properties of aqueous solutions of trimethylamine N-oxide in the temperature range from (278.15 to 323.15) K and at pressures up to 100 MPa, J. Chem. Eng. Data 60 (2015) 1291–1299, https://doi.org/10.1021/je500977g. [28] G.I. Egorov, D.M. Makarov, Densities and volume properties of (water + tertbutanol) over the temperature range of (274.15 to 348.15) K at pressure of 0.1 MPa, J. Chem. Therm. 43 (2011) 430–441, https://doi.org/10.1016/j.jct.2010.10.018. [29] A.W. Hakin, C.L. Beswick, M.M. Duke, Thermochemical and volumetric properties of aqueous urea systems. Heat capacities and volumes of transfer from water to urea— water mixtures for some 1: 1 electrolytes at 298.15 K, J. Chem. Soc. Faraday Trans. 92 (1996) 207–213, https://doi.org/10.1039/FT9969200207. [30] E.D. Kadtsyn, A.V. Anikeenko, N.N. Medvedev, Structure of aqueous solutions of trimethylaminoxide, urea, and their mixture, J. Struct. Chem. 59 (2) (2018) 347–354, https://doi.org/10.1134/S0022476618020130. [31] A.L.R. Bug, S.A. Safran, G.S. Grest, Do interactions raise or lower a percolation threshold? Phys. Rev. Lett. 55 (18) (1985) 1896–1899, https://doi.org/10.1103/ PhysRevLett.55.1896. [32] I. Balberg, N. Binenbaum, Invariant properties of the percolation thresholds in the soft-core —hard-core transition, Phis. Rev. A 35 (12) (1987) 5174–5177, https:// doi.org/10.1103/PhysRevA.35.5174. [33] M. Rottereau, J.C. Gimel, T. Nicolai, D. Durand, 3d Monte Carlo simulation of sitebond continuum percolation of spheres, Eur. Phys. J. E 11 (2003) 61–64, https:// doi.org/10.1140/epje/i2003-10006-x. [34] A. Coniglio, U. De Angelis, A. Forlani, G. Lauro, Distribution of physical clusters, J. Phys. A Math. Gen. 10 (2) (1977) 219–228, https://doi.org/10.1088/0305-4470/ 10/2/011. [35] S.B. Lee, S. Torquato, Pair connectedness and mean cluster size for continuumpercolation models: computer-simulation results, J. Chem. Phys. 89 (10) (1988) 6427–6433, https://doi.org/10.1063/1.455411. [36] D. Stauffer, A. Aharony, Introduction to Percolation Theory, Taylor & Francis, London, 1992. [37] G.D. Christian, Analytical Chemistry, 5th edition Wiley, New York, 1994. [38] S. Banerjee, J. Furtado, B. Bagchi, Fluctuating micro-heterogeneity in water - tertbutyl alcohol mixtures and lambda-type divergence of the mean cluster size with phase transition-like multiple anomalies, J. Chem. Phys. 140 (2014), 194502. https://doi.org/10.1063/1.4874637. [39] R. Gupta, G.N. Patey, Aggregation in dilute aqueous tert-butyl alcohol solutions: insights from large-scale simulations, J. Chem. Phys. 137 (2012), 034509. https://doi. org/10.1063/1.4731248.
Contents lists available at ScienceDirect
Journal of Molecular Liquids journal homepage: www.elsevier.com/locate/molliq
Statistical geometry characterization of local structure of TMAO, TBA and urea aqueous solutions☆ E.D. Kadtsyn, A.V. Anikeenko, N.N. Medvedev ⁎ Voevodsky Institute of Chemical Kinetics and Combustion, SB RAS, 630090 Novosibirsk, Russia Novosibirsk State University, 630090 Novosibirsk, Russia
a r t i c l e
i n f o
Article history: Received 14 November 2018 Received in revised form 10 April 2019 Accepted 27 April 2019 Available online 6 May 2019
a b s t r a c t Recently, we have revealed that trimethylamine-N-oxide (TMAO) molecules in aqueous solutions are distributed generally like random spheres in space (A.V. Anikeenko et al. J. Mol. Liq. 245 (2017) 35–41). We have shown that the variance of the distribution of Voronoi region volumes for the TMAO molecules and the hard random spheres are identical at the corresponding concentrations up to the mole fraction x ~ 0.1 (or packing fraction η ~ 0.2). In the current work, we study clusters in TMAO solutions and calculate the weighted mean cluster size and some topological cluster characteristics. For all these parameters, a great matching of TMAO and random spheres is also observed. The same analysis for tert-butyl alcohol (TBA) and urea aqueous solutions is carried out. These molecules have a tendency to cluster even at very low concentrations. However, basic geometrical features of the clusters are the same both for solutions and for the system of random spheres, indicating that geometrical laws governing the allocation of non-overlapping spheres in space should be taken into account when describing the structure of solutions even at rather low concentrations. © 2019 Elsevier B.V. All rights reserved.
1. Introduction The impenetrability of molecules (atoms) is usually considered as a geometrical factor that defines the size and shape of a molecule, whereas physicochemical properties of a substance are regarded to be related to the specifics of intermolecular interactions, which are divided into dispersive, electrostatic, etc. [1]. Impenetrability, however, plays a significant role in the condensed phase, determining the mutual arrangement of atoms and molecules and imposing restrictions on their movement. In this case, the details of intermolecular interactions become less important. Thus, noble gases and many metals, such as gold, copper, lead, iron, form dense crystalline structures with high symmetry, in particular, bcc or fcc, regardless of the characteristics of their outer electronic shells. The same structure is obtained in the system of hard spheres, which have no other interactions except impenetrability. A similar situation is observed for liquids. All atoms mentioned, as well as compact molecules (for example, methane), exhibit similar radial distribution functions in a liquid state. The presence of the excluded volume inevitably leads to the fact that the location of such particles is
☆ The article belongs to the special issue on The problems of solvation, complex formation and crystallization challenging the development of smart materials. ⁎ Corresponding author at: Voevodsky Institute of Chemical Kinetics and Combustion, SB RAS, 630090 Novosibirsk, Russia E-mail address: [email protected] (N.N. Medvedev).
https://doi.org/10.1016/j.molliq.2019.04.147 0167-7322/© 2019 Elsevier B.V. All rights reserved.
governed by the laws of geometry. In particular, it was proved that the densest packing of identical hard spheres is the crystal structure of fcc or hcp [2], and in the case of non-crystalline systems (dense monoatomic liquids, glasses and dense packings of spheres), the socalled “polytetrahedral packing principle” works, where tetrahedral configurations of four atoms form linear and branched clusters [3,4]. It should be emphasized that the above structural patterns are not prerogative of hard spheres. They also hold for reasonably soft particles and compact molecules with only relative sphericity. An interesting example is the substances which have the plastic crystal phase. These are the compact molecules like neopentane, carbon tetrachloride, 2,3dimethylbutane, cyclohexane [5,6]. At low temperatures, each of them has its own specific crystal structure. However, when a liquid is cooled, a plastic crystal first appears: translational diffusion ceases, and a crystalline order corresponding to a bcc or fcc structure arises. The analysis of such structures shows that the molecules are located in the lattice sites, but are oriented differently, and rotate slowly relative to each other [6,7]. In the liquid phase, such molecules behave like atoms in simple liquids. In [8], it was shown that the radial distribution functions calculated for the centers of the molecules of cyclohexane and 2,3dimethylbutane coincide with those for the Lennard-Jones system. It was also shown that the similarity is also manifested in the characteristics of the “three-dimensional” structure, studied with the help of Delaunay simplices [9]. Thus, the high symmetry of plastic crystals, where the molecules actually “hang” freely in the lattice sites, and the
2
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
universal structure of the fluids of all spherically symmetric atoms and compact molecules, is a manifestation of the geometric laws governing the placement of impenetrable spherical particles in space. In contrast to condensed systems, little attention is paid to the study of the spatial distribution of particles in systems with low density. In pure form, they are gases or fluids, whose structure is usually not discussed. However, there are similar systems that are of great interest. These are solutions whose structure and properties are also studied at low concentrations of dissolved molecules. However, the chemistry of solutions focuses on the specific interactions and the impenetrability of the dissolved molecules at relatively low concentrations is usually not taken into account. As a rule, to describe the properties of solutions, the models are used in which the excluded volume of molecules is not taken into account, for example, the Debye-Hückel model [1]. It leads, in particular, to the fact that the occurrence of associates of dissolved molecules is explained solely by their specific interactions with each other or with a solvent. Such explanations neglect that even at low concentrations, the excluded volume may influence this process. In our previous work [10], we calculated the normalized variances of the Voronoi region volumes for randomly distributed equal hard spheres at different packing fractions (degrees of space filling). Zero packing fraction corresponds to a system of random points. For that zero packing fraction, there is a wide distribution of the Voronoi volumes, because there are large fluctuations in local environments. The dispersion of volumes decreases with an increase in the packing fraction of the system (with the radius of the spheres also increasing in a given model box). The system becomes more homogeneous: instead of pairs and groups of close random points, there are now clusters of spheres with the centers separated by the distance determined by the radius of that spheres [10]. When comparing the results for random spheres with the data obtained for aqueous solutions of TMAO, we found that over a wide range of solution concentrations (up to 0.2, in terms of the packing fraction η for the solute molecules) dispersion of volumes of Voronoi regions for TMAO molecules behave the same way as for random spheres at corresponding concentrations. This fact can be explained by the TMAO molecules being distributed in the solution independently, thus not exhibiting the specific interaction with each other. The question of why the TMAO molecules behave independently in solution remains open. These molecules have a large dipole moment, form hydrogen bonds, and affect water by their hydrophobic part, thus being able to interact both directly with each other and indirectly through water [11,12]. Apparently, the interactions of different nature somehow compensate each other. This question is worth being studied separately, especially, as we show that the similarity of the spatial distribution of TMAO in solution and random hard spheres is manifested not only in the dispersion of the volumes of the Voronoi regions but also in the structure of the emerging clusters. In this paper, in addition to TMAO solutions, we also study the solutions of tert-butanol and urea. These molecules tend to associate in aqueous solutions, and they have a high first peak of pair correlation function even at very low concentrations. At the same time, the mechanisms of association in TBA and urea are different. In the case of TBA, it is generally accepted that it is a hydrophobic interaction [13–15]. It is believed that water “expels” TBA molecules and they create relatively compact associates. In the case of urea, the situation is less clear. The urea molecule does not have a hydrophobic group, but hydrogen bonding is possible. The urea molecules are thought to be included in the structure of the hydrogen bonding network [16–19].
1 bar. A stochastic velocity rescaling thermostat [21] and ParinelloRahman barostat [22] were used. Integration step was 2 fs. The production runs were 100 ns length. Every 10th picosecond frame was used for analysis (104 independent configurations for each model). Parameters of force fields were taken from [23] for TMAO, from [24] for TBA, and for urea from [25]. TIP4P water was used [26]. The structure of these molecules is shown in Fig. 1. We prepared twenty models of TMAO and TBA aqueous solutions on the molar percentage interval from 0.5% to 10% with the step of 0.5%. Each model contained 100 solute molecules in the simulation box. For the lowest concentration, the box had an edge length of 8.5 nm and contained 19,900 water molecules, see Table 1. The more concentrated solutions were obtained by removing the required number of water molecules and decreasing the box size [10]. For urea, we also got twenty models but on the interval from 1% to 20% with the step of 1%. For the concentrations from 1 to 10% we used 100 urea molecules. However, for the interval 11–20% we used 200 molecules to keep the box size large enough at the higher concentrations. The larger mole fraction of urea was required because the urea molecule is smaller. So further we used the packing fraction as a measure of concentration. The relation between molar fraction, molar concentration, packing fraction and the corresponding number of water molecules in the models are shown in Table 1. The comparison of the calculated and experimental densities of the solutions is shown in Fig. 2. One can see quite a good agreement for all models for used intervals of concentration.
2. Methods
¼ NπD 6V box , where D. is the mean diameter of the spheres. It is easy to see
2.2. Systems of random spheres All models of random spheres were created by random sequential addition of hard spheres into empty simulation box with periodic boundary conditions. For cluster analysis, we have used only systems of monodisperse spheres. The diameter of these spheres was chosen to be 0.53 nm, and the packing fraction was determined by the variation of the size of the box. (Remember that packing fraction is the ratio of the 3
volume occupied by spheres to the volume of the system: η ¼ NπD 6V box , where Vbox is the volume of the model box, N is the number of spheres, and D is their diameter). In this work, we have used N = 100. For each packing fraction, 104 independent configurations for the similarity with solutions models in the statistical analysis were obtained. The systems of monodisperse spheres were obtained for packing fraction from η = 0.01 to 0.21 with the step of 0.01. We have also obtained the models of polydisperse spheres with Gaussian diameter distribution. These models were used for determining the effective radius of solute molecules, see below. We have varied both mean diameter D and its variance. Fig. 3 shows the pair correlation functions g(r) for the system of monodisperse and polydisperse spheres with the same diameter. One can see that vertical left side which is intrinsic for monodisperse spheres became smoother for polydisperse spheres. Note that the middle of this front corresponds well to the mean diameter of spheres at least for low concentrations. Such observation confirms reasonability of our suggestion used in [10], where values of effective diameters of TMAO and TBA were defined as the position of the middle of the left front of the first peak of g(r). Note that if σ is small enough (a few percent of diameter) then the packing fraction for the system of such polydisperse spheres systems can be calculated by using the formula for the monodisperse spheres η 3
2.1. Solutions. Molecular dynamics simulations The molecular dynamics simulation protocol was the same as in [10], using Gromacs 5.0 package [20] at temperature 300 K and pressure
that contributions from the spheres of the sizes being greater or less than the mean value compensates each other with the precision of the second order of smallness, which can be neglected. To estimate the volume of solute molecules in solution we represent them as spheres, see below. On the other hand, as is well known in
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
3
Fig. 1. Structure and van der Waals representation of TMAO, TBA, and urea molecules (from left to right).
chemistry and biology, the volume of a molecule can be calculated using van der Waals surface of atoms, or the molecular (Connolly) surface. Our choice is because we use a system of hard spheres as a reference. 2.3. The effective diameter of molecules Real molecules have a hard core, but they are not rigid themselves. To determine the effective diameter of the molecules, we use the models of random polydisperse spheres. The values of mean diameter D and its variance are chosen in such a way that the left front of the first peak of g(r) for the spheres system could match the left front of g (r) of the solute molecules studied. Note that in such comparison, the same number of molecules and packing fraction of the systems is used. Fig. 4 shows the g(r) for TMAO solution and random polydisperse spheres at different packing fractions. In both cases, the mean diameter of the spheres is 0.53 nm. The dispersion at both concentrations is also equal (σ = 0.035 nm) because the position and the steepness of the peak do not depend on the concentration. Note that in the previous work [10] we obtained the same value of molecule diameter. However, in the current work, we use the value of D, which was obtained from the fitting of the whole left front. We do not pay attention now to the area beyond the first peak. The curves of g(r) have different behavior there. For the random spheres at
the used packing fraction, far oscillations still do not appear, while for the solutions clear oscillations are observed. They are a result of “matrix influence”. In other words, they are determined by the structure of the surrounding water [30]. When determining the effective diameters of TBA and urea molecules, the additional question arises related to the fact that they have a high first peak of g(r). However, even in this case, the left front reflecting the variations of closest molecular contacts can be well described by using our approach with polydisperse spheres. The height of the first peak is caused by the larger number of the nearest solute molecules. The increase in the number of contacts is caused by the structure of solution, interactions with solvent and is not connected with direct intermolecular contacts and with the “hardness” of the molecules. In other words, the left front of the first peak can be considered without regard to its height. Therefore, we simply reduced the height of the first peak of the solution to the value of the system of random spheres by applying the corresponding normalizing factor. Fig. 5 presents the fitting of the left front of TBA and urea. We have obtained the values D = 0.52 nm for TBA and D = 0.39 nm for urea. For a higher concentration, the same good description is obtained, giving the same values of the effective diameter of the molecules. For these solutions, as well as for TMAO, the position and steepness of the left front remain almost unchanged in the concentration range used.
Table 1 3
Correspondence between molar fraction (x), molar concentration (C) and packing fraction (η), for TMAO, TBA and urea aqueous solutions. Packing fraction was calculated as η ¼ NπD 6V box , where Vbox is the volume of the model box, N – the number of spheres, and D is the effective diameter of the molecule: DTMAO = 0.53 nm, DTBA = 0.52 nm, Durea = 0.39 nm, see text. TMAO x, % C, mol/l pf,% Water molecules x, % C, mol/l pf, % Water molecules
0.5 0.273 1.31 19,900 5.5 2.613 12.25 1718
1.0 0.537 2.58 9900 6.0 2.815 13.16 1566
1.5 0.794 3.8 6566 6.5 3.015 14.04 1438
2.0 1.044 4.99 4900 7.0 3.205 14.89 1328
2.5 1.287 6.13 3900 7.5 3.393 15.71 1233
3.0 1.523 7.24 3233 8.0 3.575 16.49 1150
3.5 1.752 8.31 2757 8.5 3.756 17.28 1076
4.0 1.976 9.34 2400 9.0 3.927 18.03 1011
4.5 2.194 10.34 2122 9.5 4.083 18.77 952
5.0 2.407 11.31 1900 10.0 4.248 19.46 900
0.5 0.272 1.21 19,900 5.5 2.537 11.24 1718
1.0 0.534 2.37 9900 6.0 2.726 12.08 1566
1.5 0.788 3.49 6566 6.5 2.908 12.89 1438
2.0 1.033 4.58 4900 7.0 3.084 13.67 1328
2.5 1.270 5.63 3900 7.5 3.254 14.42 1233
3.0 1.500 6.65 3233 8.0 3.417 15.14 1150
3.5 1.721 7.63 2757 8.5 3.580 15.87 1076
4.0 1.935 8.58 2400 9.0 3.736 16.56 1011
4.5 2.143 9.5 2122 9.5 3.889 17.24 952
5.0 2.343 10.38 1900 10.0 4.031 17.87 900
7 3.52 6.58 1328 17 7.55 14.12 976
8 3.97 7.42 1150 18 7.9 14.77 911
9 4.41 8.24 1011 19 8.25 15.42 852
10 4.84 9.04 900 20 8.58 16.04 800
TBA x, % C, mol/l pf, % Water molecules x, % C, mol/l pf, % Water molecules Urea x, % C, mol/l pf, % Water molecules x, % C, mol/l pf, % Water molecules
1 0.55 1.02 9900 11 5.25 9.82 1618
2 1.08 2.01 4900 12 5.66 10.58 1466
3 1.59 2.98 3233 13 6.06 11.33 1338
4 2.09 3.91 2400 14 6.44 12.05 1228
5 2.58 4.83 1900 15 6.82 12.75 1133
6 3.06 5.72 1566 16 7.19 13.44 1050
4
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
system and the appearing clusters are shown in Fig. 7. The particles intersecting by their permeable shells are combined into the cluster.
3. Results
Fig. 2. The densities of aqueous solutions of urea, TMAO and TBA at T = 300 K and P = 1 bar vs. concentration. Solid lines refer to experimental data [27–29], large symbols show the simulation values.
Note that the value obtained for the TBA diameter slightly differs from the one for TMAO. In the [10] we used the same value D = 0.53 nm for both solutes for simplicity. However, such a distinction is not significant for our analysis. A small shoulder on the left slope of the TBA curve is due to the hydrogen bond formation between some TBA molecules. Fig. 5а shows that the fraction of such molecules is small, about 2–3%. This fact was also proven by our special analysis. Therefore, we do not pay any particular attention to these molecules.
2.4. Selection of the critical distances for cluster determination The choice of critical distances between the molecules to include them into one cluster is always arbitrary. In this work, we associate this distance with the molecular diameter D and use three values: 1.25 D, 1.30 D and 1.35 D. On the Fig. 6 all functions g(r) for our systems are represented in the units of D. One can see that critical distances used in our study go through the main part of the first peak, see vertical dotted lines. Recall that from the mathematical standpoint our systems are the systems of “hard spheres with permeable shells”, which used for studying the continuous percolation years ago [31–33]. An illustration of such
Fig. 8 demonstrates variances of Voronoi volumes for solute molecules in our solutions and in the hard spheres system at different packing fractions. Recall that Voronoi regions were calculated for the system of centers of the solute molecules, see details in [10]. In this case, we describe a general motif of the spatial distribution of the molecules. In current work, we represent the data for the new models and show the results for urea models, containing 100 and 200 solute molecules, see Table 1. The most surprising result is that curves for TMAO solutions and random spheres match almost exactly. TBA and urea curves show different behavior. We have related the increase of variance for TBA with the arising of inhomogeneity in its solutions. For urea, however, the variance decreases monotonically. It means that urea solutions become more homogeneous but to a lesser extent than in the case of TMAO or random spheres. In the following, we present the results of a cluster analysis which allow studying our systems from another side, namely from the standpoint of local structure. First, we calculated the simple mean cluster size, which shows the average number of particles in the cluster. It is defined as 〈i〉 = ∑i=1Nini/∑i=1Nni, where ni is the number of clusters containing i particles and the sum is taken over all existing clusters, which can include also single molecules (i = 1) or cluster formed by all solute molecules (i = N). It can also be calculated using standard gmx clustsize utility from the Gromacs package. This value is shown in Fig. 9 for our systems in dependence on packing fraction for different critical distance values. It should be emphasized that TMAO molecules in the solutions behave again like random spheres: black and red curves match for all critical distances. The clusters of TBA and urea differ from the TMAO ones and the random spheres. The larger size of these clusters can be explained by the higher number of the closest solute molecules in solutions. Indeed, both TBA and urea have high first peaks of g(r). This simple parameter, however, has a low sensibility for the cluster analysis. In the percolation theory, another mean value is commonly used: the second moment of the cluster size distribution normalized by the first moment: S = ∑i=1Ni2ni/∑i=1Nini [34–36]. In fact, it is a weighted (by the number of particles) the mean size of the cluster. Indeed, this formula can be rewritten as S = ∑i=1Nfii, where fi is the weight of the cluster with size i. It is defined as f i ¼ inNi , where ini is the number of particles in the all clusters of a given size i, and in the
Fig. 3. Function g(r) for systems of monodisperse random spheres of unit diameter (black) and for polydisperse random spheres with the same mean diameter and different dispersions: 0.05 (blue) and 0.03 (green). Models packing fraction η = 1% (a), and η = 15% (b).
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
5
Fig. 4. Comparison of g(r) for TMAO molecules in solution (red) and random polydisperse spheres system (black) with D = 0.53 nm, σ = 0.025 nm. TMAO concentration x = 0.05, spheres packing fraction η = 0.01 (a). TMAO x = 0.075, η = 0.15 (b), see Table 1.
denominator there is a full number of particles in the system, which was previously written as the sum of the particles over all clusters: ∑i=1Nini = N. In such definition of the weighted mean size, the more massive size clusters (which are always not high in number) are not suppressed by small clusters (which are always high in number). Note that such a parameter is used in chemistry in the study of complexation in the solutions, where the weight factor fi is called form fraction of i. [37]. In the percolation theory, one deals with infinitely large systems. In this case, the weighted mean size is calculated only for finite clusters. It is known, its value diverges according to a power law near the percolation threshold. However, for relatively small models, the percolation threshold is not well defined (large fluctuations appear at its determination) [36]. For this reason, we do not calculate the percolation cluster in our models and use the cluster size only for the comparison of the different systems. Fig. 10 shows the weighted mean cluster size normed by the full number of particles in the model. (Our models contain different number of urea molecules N. For comparing all models together, we use the normed cluster size S/N). The curves aim to the limit value of the unit with the increasing concentration. It means that almost all the solute molecules are included in one cluster at the high concentration in our solutions, namely at the packing fraction η ≈ 0.2. Therefore, our systems are beyond the percolation threshold at the given concentrations. In the work [33], percolation in the system of hard spheres with permeable shells was thoroughly studied with Monte-Carlo technique. The models with different shell thickness and packing fraction were investigated. For the shells corresponding to our critical distances (1.25, 1.30,
1.35), it was found that the percolation occurs at the packing fractions η are equal to 0.175, 0.15, 0.13 respectively. In Fig. 10, the vertical dotted lines mark these values. One can see that in all cases it corresponds approximately to the value of the normed mean size S/N ≈ 0.3. It seems reasonable to use this result for defining the critical concentration for our solutions. For TBA solutions, such a value of mean cluster size corresponds to the packing fraction close to 0.10 (at all critical distances). From Table 1, one can see that this value corresponds to the mole fraction х = 0.048. This value is in agreement with existing data. In the work [38], it was observed the percolation threshold of the TBA molecules at x ≈ 0.05. The authors obtained “islands” of TBA even at the lower mole fraction, while the large spanning cluster emerged above that mole fraction. A huge model of TBA solution was studied in [39]. It was noted that TBA aggregation could be detected at x ≈ 0.03 and a rapid increase in the size of the clusters is observed in the range x ≈ 0.03–0.06. By applying our indirect method, we can give even a more accurate value of the percolation threshold in TBA solutions. For the urea, the behavior of the clusters differs from the behavior for TBA. One can also see a difference in the various critical distances. It means a different cluster nature for urea solutions. Nevertheless, one can say that the percolation threshold for the urea solution is approximately at η = 0.12, corresponding to its molar fraction х = 0.14, that is much higher than for TBA. Unfortunately, we do not know any works that would study the percolation in urea solutions. Look at to the slight break at η ≈ 0.09 on the urea curves. It occurs because of a change in the size of the model. At a higher concentration of urea, we use a twice bigger box and worked with 200 solute
Fig. 5. Comparison of g(r) for TBA at x = 0.005 (a) and urea at x = 0.01 (b) with the random polydisperse spheres system at packing fraction of η = 0.01. The functions for TBA and urea are scaled to get the unit height of the first peak, see text. Fitting parameters are D = 0.52 nm, σ = 0.02 nm for TBA and D = 0.39 nm, σ = 0.017 nm for urea.
6
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
Fig. 6. Functions g(r) for TMAO, TBA, and urea in the units of the effective diameter of the molecule. Black curve corresponds to the random polydisperse spheres system. Red curves correspond to TMAO, blue to TBA, and green to urea. Vertical dashed line shows an effective diameter equal to the unit; black dotted lines show suggested critical distances for cluster selection.
molecules instead of 100 for low concentration (see above). The displacement of curves means that cluster composition in our models slightly differs. Indeed, for our rather small models, the fraction of bigger clusters and the fraction of single molecules in different sized boxes can be different. Fig. 9 does not show such a break at the change of the model size. That can be because of the low structural sensibility of that parameter - the simple mean size of the clusters. Fig. 11 demonstrates the mean degree of the nodes (the average number of the neighbors for a given molecule in the cluster) in dependence on the cluster size. We use here only rather small clusters that involve no more than 10 molecules. At relatively low concentration, there are quite a lot of such clusters in the solutions. As the concentration increases, their number significantly decreases due to the appearance of large clusters. For this reason, we work here only with small packing fractions: 0.04, 0.05 and 0.06 for hard spheres and with values for solutions close to them, see Table 1. Thus, for each solution and for random spheres system, three curves are shown corresponding to the given packing fractions, Fig. 11. In addition, a theoretical curve for tree clusters (cluster without cycles) is shown, which was calculated by the formula hdi ¼ 2−2 i .
Fig. 8. The variance of normed Voronoi region volume distributions for TBA (blue) urea (green) and TMAO (red) solute molecules in solution compared with random spheres (black). A break at η = 0.09 on the urea curve is due to size effect, see text.
Note that the degree of the nodes almost does not depend on the chosen critical distance, compare the figures (a)–(c). However, there is a clear difference between the systems. The highest degree of the nodes is for the TBA solutions (blue curves), the least one is for TMAO solutions and for the random spheres system. We emphasize that the curves for TMAO and for random spheres coincide very well also for these parameters. Note that for clusters formed by two particles, the degree of the node is 1. With the number of particles in a cluster increasing, the degree also increases. For trees, it approaches the theoretical limit equal to 2. However, for TBA solutions, it approaches the value 2.5, indicating the presence of cycles. The greater the mean degree of the nodes, the more cycles in the cluster of a given size. For the simplest case of a cluster of a size 3 there are only linear clusters or triangles. The fraction of triangles among all clusters of size 3 is shown in Fig. 12 as a function of concentration. The smallest amount of triangles is in the TMAO solutions and random spheres. In the TBA solutions, their fraction is more than twice greater. In the urea solutions, the fraction of triangles is smaller than in TBA solutions, suggesting again that the nature of clusters is different for these two systems. Note that this parameter also shows the structure
Fig. 7. Schematic representation of hard spheres with the permeable shell. The black areas define packing fraction η for the system (a). The clusters corresponding to the system (b).
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
7
Fig. 9. The mean size of clusters in solutions and the random sphere system depending on the packing fraction for different critical distances (see text): 1.25D (a), 1.3D (b) and 1.35D (c).
Fig. 10. The weighted mean size of clusters for different critical distances, see caption in Fig. 9. The vertical dotted lines mark the percolation thresholds for the corresponding random hard spheres with permeable shells obtained by Monte-Carlo technique [33], see text.
8
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
Fig. 11. Mean degree of the nodes in the cluster for clusters of different size for critical distances 1.25D (a), 1.3D (b) and 1.35D (c). Three groups of curves belong TMAO, TBA and urea solutions. Three curves in each group corresponding to the concentrations η ≈ 0.03, 0.04 and 0.06. Black curves represent theoretical values for trees, i.e. clusters without cycles.
closeness of TMAO solutions to the random sphere system, see black and red curves in Fig. 11. 4. Conclusions In this work, we have studied clusters of TMAO, TBA and urea molecules in aqueous solutions at different concentrations. In our previous work [10] we revealed that TMAO molecules are distributed in the solutions like the random non-overlapped spheres up to the packing fraction of η = 0.2. Now we have shown that TMAO clusters are also the
Fig. 12. The fraction of triangles among all clusters of size 3 depending on concentration. The results for critical distance 1.3D are shown.
same as clusters in the random spheres system. The question of the reason for such behavior of TMAO molecules remains open. They can interact with each other both directly and through the water due to having the large dipole moment, forming strong hydrogen bonds and having large hydrophobic moiety. The question of how all these interactions can be mutually compensated is worthy of a separate study. We want to emphasize that the similarity with random sphere system can be observed in a rather wide range of concentrations and manifests itself as in the spatial distribution of molecules as well as the features of intermolecular contacts. The distribution of TBA and urea in aqueous solution is different from one of random spheres. It is apparently the consequence of the specific intermolecular interactions that are not compensated like in TMAO solutions. The global distribution of these molecules is more inhomogeneous, and small clusters are larger and have more cycles than in the random spheres at the same concentrations. We estimated the percolation threshold in TBA and urea solutions is lower than in the random system. The reason is that these molecules begin to associate even at very low concentrations. However, the nature of the interactions of TBA and urea molecules is different, leading to a more homogeneous distribution of urea molecules and less number of cycles in the urea clusters. An interesting difference in the structure of these two solutions was noticed in the work [30]. It was shown that the height of the first peak of g(r) for urea decreases with concentration, while increases for TBA. It was suggested that urea clusters contain fewer cycles (are more ramified) than in the TBA case. Our analysis is in agreement with that explanation. However, despite the differences between TBA and urea solutions and their difference from TMAO and random spheres, there are common features of the cluster behavior with concentration for all these systems. The main cluster parameters follow to general geometrical laws for the systems with excluded volume, which is manifesting themselves even at low concentrations. Thus, geometric properties of the
E.D. Kadtsyn et al. / Journal of Molecular Liquids 286 (2019) 110870
arrangement of impermeable molecules should be always taken into account when studying the structure of solutions. Despite the difference in the properties of clusters for TBA and urea from the system of hard spheres, it should be noted a similar behavior of the curves in Figs. 9–12 for all these systems. The similarity is obviously a consequence of the existence of an excluded volume, i.e. mutual impermeability of molecules. This means that when explaining the structure of solutions, it is necessary to take into account the geometrical laws of the spatial arrangement of particles even at rather low concentrations. Acknowledgments This work was supported by Russian Foundation for Basic Research (project No. 18-03-00045) and partly by Russian Ministry of Science and Education under 5-100 Excellence Programme. References [1] P. Atkins, J. de Paula, Physical Chemistry: Thermodynamics, Structure, and Change, W. H. Freeman and Company, New York, 2002. [2] T.C. Hales, A proof of the Kepler conjecture, Ann. Math. (2005) 1065–1185. [3] J.D. Bernal, The Bakerian lecture, 1962. The structure of liquids, Proc. R. Soc. Lond. A Math. Phys. Sci. 280 (1382) (1964) 299–322, https://doi.org/10.1098/rspa.1964. 0147. [4] A.V. Anikeenko, N.N. Medvedev, Polytetrahedral nature of the dense disordered packings of hard spheres, Phys. Rev. Lett. 98 (23) (2007) 235504, https://doi.org/ 10.1103/PhysRevLett.98.235504. [5] C.N.R. Rao, Molecular motion in plastic crystals, Proc. Indian Acad. Sci. Chem. Sci. 94 (1) (1985) 181–199, https://doi.org/10.1007/BF02841265. [6] H. Farman, J.C. Dore, M.C. Bellisent-Funel, D.G. Montague, Structural studies of cyclohexane, C6D12 by neutron diffraction. II the plastic crystal phase, Mol. Phys. 73 (4) (1991) 855–871, https://doi.org/10.1080/00268979100101611. [7] J.H. Strange, Pressure and temperature dependence of molecular motion in organic plastic crystals, J. Chem. Phys. 68 (7) (1978) 3078–3088, https://doi.org/10.1063/1. 436174. [8] A.V. Anikeenko, A.V. Kim, N.N. Medvedev, Molecular dynamics simulation of the structure of c6 alkanes, J. Struct. Chem. 51 (6) (2010) 1090–1096, https://doi.org/ 10.1007/s10947-010-0167-z. [9] A.V. Anikeenko, A.V. Kim, N.N. Medvedev, Delaunay simplexes in liquid cyclohexane, Proceedings of Sixth International Symposium on Voronoi Diagrams in Science and Engineering (ISVD 2009), 23–26 July, Technical University of Denmark, Denmark 2009, pp. 271–277, https://doi.org/10.1109/ISVD.2009.10 , published by IEEE Computer Society. [10] A.V. Anikeenko, E.D. Kadtsyn, N.N. Medvedev, Statistical geometry characterization of global structure of TMAO and TBA aqueous solutions, J. Mol. Liq. 245 (2017) 35–41, https://doi.org/10.1016/j.molliq.2017.06.001. [11] L. Larini, J.E. Shea, Double resolution model for studying TMAO/water effective interactions, J. Phys. Chem. B 117 (2013) 13268–13277, https://doi.org/10.1021/ jp403635g. [12] M.V. Fedotova, S.E. Kruchinina, G.N. Chuev, Hydration structure of osmolyte TMAO: concentration/pressure-induced response, New J. Chem. 41 (2017) 1219, https:// doi.org/10.1039/c6nj03296f. [13] H.S. Frank, M.W. Evans, Free volume and entropy in condensed systems III. Entropy in binary liquid mixtures; partial molal entropy in dilute solutions; structure and thermodynamics in aqueous electrolytes, J. Chem. Phys. 13 (11) (1945) 507–532, https://doi.org/10.1063/1.1723985. [14] G. Onori, A. Santucci, Dynamical and structural properties of water/alcohol mixtures, J. Mol. Liq. 69 (1996) 161–181, https://doi.org/10.1016/S0167-7322(96)90012-4. [15] A. Di Michele, M. Freda, G. Onori, M. Paolantoni, A. Santucci, P. Sassi, Modulation of hydrophobic effect by cosolutes, J. Phys. Chem. B 110 (42) (2006) 21077–21085, https://doi.org/10.1021/jp068055w. [16] Y.L.A. Rezus, H.J. Bakker, Destabilization of the hydrogen-bond structure of water by the osmolyte trimethylamine N-oxide, J. Phys. Chem. B 113 (13) (2009) 4038–4044, https://doi.org/10.1021/jp805458p.
9
[17] J.K. Carr, L.E. Buchanan, J.R. Schmidt, M.T. Zanni, J.L. Skinner, Structure and dynamics of urea/water mixtures investigated by vibrational spectroscopy and molecular dynamics simulation, J. Phys. Chem. B 117 (42) (2013) 13291–13300, https://doi.org/ 10.1021/jp4037217. [18] Y.L.A. Rezus, H.J. Bakker, Effect of urea on the structural dynamics of water, PNAS 103 (49) (2006) 18417–18420, https://doi.org/10.1073/pnas.0606538103. [19] D. Bandyopadhyay, S. Mohan, S.K. Ghosh, N. Choudhury, Molecular dynamics simulation of aqueous urea solution: is urea a structure breaker? J. Phys. Chem. B 118 (40) (2014) 11757–11768, https://doi.org/10.1021/jp505147u (2014). [20] S. Pall, M.J. Abraham, C. Kutzner, B. Hess, E. Lindahl, Tackling exascale software challenges in molecular dynamics simulations with GROMACS, International Conference on Exascale Applications and Software, 1, Springer, Cham 2014, April, pp. 3–27, https://doi.org/10.1007/978-3-319-15976-8_1. [21] G. Bussi, D. Donadio, M. Parrinello, Canonical sampling through velocity rescaling, J. Chem. Phys. 126 (2007), 014101. https://doi.org/10.1063/1.2408420. [22] M. Parrinello, A. Rahman, Polymorphic transitions in single crystals: a new molecular dynamics method, J. Appl. Phys. 52 (1981) 7182, https://doi.org/10.1063/1. 328693. [23] C. Hölzl, P. Kibies, S. Imoto, R. Frach, S. Suladze, R. Winter, ... S.M. Kast, Design principles for high–pressure force fields: aqueous TMAO solutions from ambient to kilobar pressures, J. Chem. Phys. 144 (14) (2016) 144104, https://doi.org/10.1063/1. 4944991. [24] C. Caleman, P.J. van Maaren, M. Hong, J.S. Hub, L.T. Costa, D. van der Spoel, Force field benchmark of organic liquids: density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant, J. Chem. Theory Comput. 8 (1) (2011) 61–74, https://doi.org/10. 1021/ct200731v. [25] S. Weerasinghe, P.E. Smith, A Kirkwood− buff derived force field for mixtures of urea and water, J. Phys. Chem. B 107 (16) (2003) 3891–3898, https://doi.org/10. 1021/jp022049s. [26] J.L.F. Abascal, C. Vega, A general purpose model for the condensed phases of water: TIP4P/2005, J. Chem. Phys. 123 (2005), 234505. https://doi.org/10.1063/1.2121687. [27] D.M. Makarov, G.I. Egorov, A.M. Kolker, Density and volumetric properties of aqueous solutions of trimethylamine N-oxide in the temperature range from (278.15 to 323.15) K and at pressures up to 100 MPa, J. Chem. Eng. Data 60 (2015) 1291–1299, https://doi.org/10.1021/je500977g. [28] G.I. Egorov, D.M. Makarov, Densities and volume properties of (water + tertbutanol) over the temperature range of (274.15 to 348.15) K at pressure of 0.1 MPa, J. Chem. Therm. 43 (2011) 430–441, https://doi.org/10.1016/j.jct.2010.10.018. [29] A.W. Hakin, C.L. Beswick, M.M. Duke, Thermochemical and volumetric properties of aqueous urea systems. Heat capacities and volumes of transfer from water to urea— water mixtures for some 1: 1 electrolytes at 298.15 K, J. Chem. Soc. Faraday Trans. 92 (1996) 207–213, https://doi.org/10.1039/FT9969200207. [30] E.D. Kadtsyn, A.V. Anikeenko, N.N. Medvedev, Structure of aqueous solutions of trimethylaminoxide, urea, and their mixture, J. Struct. Chem. 59 (2) (2018) 347–354, https://doi.org/10.1134/S0022476618020130. [31] A.L.R. Bug, S.A. Safran, G.S. Grest, Do interactions raise or lower a percolation threshold? Phys. Rev. Lett. 55 (18) (1985) 1896–1899, https://doi.org/10.1103/ PhysRevLett.55.1896. [32] I. Balberg, N. Binenbaum, Invariant properties of the percolation thresholds in the soft-core —hard-core transition, Phis. Rev. A 35 (12) (1987) 5174–5177, https:// doi.org/10.1103/PhysRevA.35.5174. [33] M. Rottereau, J.C. Gimel, T. Nicolai, D. Durand, 3d Monte Carlo simulation of sitebond continuum percolation of spheres, Eur. Phys. J. E 11 (2003) 61–64, https:// doi.org/10.1140/epje/i2003-10006-x. [34] A. Coniglio, U. De Angelis, A. Forlani, G. Lauro, Distribution of physical clusters, J. Phys. A Math. Gen. 10 (2) (1977) 219–228, https://doi.org/10.1088/0305-4470/ 10/2/011. [35] S.B. Lee, S. Torquato, Pair connectedness and mean cluster size for continuumpercolation models: computer-simulation results, J. Chem. Phys. 89 (10) (1988) 6427–6433, https://doi.org/10.1063/1.455411. [36] D. Stauffer, A. Aharony, Introduction to Percolation Theory, Taylor & Francis, London, 1992. [37] G.D. Christian, Analytical Chemistry, 5th edition Wiley, New York, 1994. [38] S. Banerjee, J. Furtado, B. Bagchi, Fluctuating micro-heterogeneity in water - tertbutyl alcohol mixtures and lambda-type divergence of the mean cluster size with phase transition-like multiple anomalies, J. Chem. Phys. 140 (2014), 194502. https://doi.org/10.1063/1.4874637. [39] R. Gupta, G.N. Patey, Aggregation in dilute aqueous tert-butyl alcohol solutions: insights from large-scale simulations, J. Chem. Phys. 137 (2012), 034509. https://doi. org/10.1063/1.4731248.