Molecular dynamics simulation of a binary mixture near the lower critical point

2,6-lutidine molecules mix with water at high and low temperatures but in a wide intermediate temperature range a 2,6-lutidine/water mixture exhibits a miscibility gap. We constructed and validated an atomistic model for 2,6-lutidine and performed molecular dynamics simulations of 2,6-lutidine/water mixture at different temperatures. We determined the part of demixing curve with the lower critical point. The lower critical point extracted from our data is located close to the experimental one. The estimates for critical exponents obtained from our simulations are in a good agreement with the values corresponding to the $3D$ Ising universality class.


I. INTRODUCTION
It is well recognised that solutions of water with organic molecules may show closedlooped phase diagrams with a miscibility gap. The occurrence of such phase diagrams with more than one critical point is very different from what is observed in mixtures of simple fluids. In simple fluids, the two fluids form a homogeneous mixed phase at higher temperatures, whereas they phase separate at lower temperatures. The appearance of an upper critical point (UCP) terminating the two-phase coexistence results from a competition between the entropy of mixing and the energy. In the case of mixtures of complex species, such as water and some organic molecules, the mechanism behind the occurrence of the miscibility gap is rather involved. As argued in Refs. [1,2], the existence of a lower critical point (LCP) can be due to formation of directional bonds between water and the organic molecules, e.g., hydrogen bonding. Below the LCP, the strong hydrogen bonding promotes mixing at the expense of rotational degrees of freedom of molecules, which are effectively "frozen out" by the bonding. At higher temperatures, thermal fluctuations "unfreezes" these rotational degrees of freedom so that the hydrogen bonding gets destroyed and the mixed phase separates. Theoretical studies of simple lattice models and within the Landau-Ginzburg approach have demonstrated that the miscibility gap in liquid mixtures may indeed emerge due to angular dependent attractive interactions on top of the spherically symmetric ones [2][3][4][5][6]. Experiments show that the shape of closed-loop phase diagrams of aqueous solutions of organic molecules near the LCP is very flat, i.e., concentrations of the species in the two coexisting phases near LCP vary strongly with temperature. Such behavior indicates that in both coexisting liquid phases some local structures are formed which are determined by hydrogen bonding between water and the organic molecules. The features of the miscibility gap, in particular its sensitivity to changes of the intermolecular interactions, has been studied theoretically [7], by computer simulation (for tetrahydrofuranwater mixtures) [9,10], and experimentally [7,8] by several techniques, e.g. by deuteration of water, addition of electrolytic impurities or hydrotrops. As follows from these studies, the details of the interactions affect the critical temperature as well as the detailed shape of the phase diagram.
The aqueous solution of 2,6-lutidine (2,6-dimethylpyridine) is a binary liquid mixture, which has gained considerable attention in context of its wetting behaviour at silica walls [11], porous glass [12], and colloids [13,14]. Of particular interest has been the effect of temperature changes on the reversible aggregation of colloidal particles dispersed in a 2,6-lutidine/water mixture [15][16][17]. More recently, 2,6-lutidine/water mixtures were used to determine critical Casimir interactions in colloidal systems [18][19][20][21]. The reason for the popularity of 2,6-lutidine/water mixture in these studies is that it possesses a closed loop phase diagram with a relatively wide temperature miscibility gap, i.e., the difference between the upper and lower critical point is large (≈197 • C). Further, the LCP is conveniently located near room temperature. The phase diagram as well as static and dynamic critical properties of 2,6-lutidine/water were studied intensively experimentally [22][23][24][25][26]. The LCP has been reported to occur at T c ≈ 307.1 K and at the lutidine mole fraction x lut ≈ 6.1% [22,24,25,27].
The location of the LCP is, however, strongly affected by impurities and, moreover, it is sensitive to the method of determination. This is why the values of the lower critical temperature and the critical mass, volume or mole fraction vary in the literature [26,[28][29][30][31].
Such uncertainty is especially troublesome for application of a 2,6-lutidine/water liquid mixture as solvent in colloidal systems tuned by critical Casimir interactions, where the precise knowledge of the deviation in temperature and concentration from the critical values is required. Computer simulations of the 2,6-lutidine/water mixture are thus highly desirable.
Moreover, such simulations can help to better understand the molecular mechanisms behind the lower critical point, and are a necessary prerequisite for studies of more complicated phenomena such as, for example, formation of of mesostructures in a mixture of water/organic solvent by adding an antagonistic salt, which is composed of hydrophilic cation and hydrophobic anion. Such mesostructures were observed recently in SANS experiments in a mixture of water/3-methylpyridine/NaBPh 4 near the lower critical point [32] and off-critical point [33,34]. A similar observation is reported for 2,6-lutidine/water mixture [35]. These phenomena are only partially explained theoretically [36,37]. Another interesting topic, which can be an extension of the present work is critical adsorption of 2,6-lutidine/water mixture containing salt (inorganic one) at a charged and selective wall. This phenomenon is of crucial relevance to recent experiments on critical Casimir interactions [18,20,38]. The findings of these experiments were not yet clarified in a satisfactory manner and there is some controversy about their origin. Some analytical studies of this problem are available in the literature [39][40][41]. However, these results were obtained within a mesoscopic model and by using various approximations and thus they need to be verified by means of microscopic simulations.
In the present paper we propose an atomistic description of the 2,6-lutidine molecule and apply it to study the bulk 2,6-lutidine liquid as well as the 2,6-lutidine/water mixture near the LCP by molecular dynamics simulations. The goal is to check whether our model of the 2,6-lutidine molecule is able to capture the main features of the bulk fluid as well as those of the aqueous solution. All simulations were performed by the Gromacs/4.6.7 package [42]. The Gromos54a7 force field [43] was applied for Lennard-Jones (LJ) pair potential parameters, bond lengths and bonded parameters for angles, and dihedrals. The Particle Mesh Ewald (PME) approach [44] was applied for electrostatic interactions, while a cut-off length r c = 1.2 nm was applied to the LJ interactions. Simulations in this work were either performed in the NpT ensemble (constant pressure) or NVT ensemble (constant volume); the type of ensemble is mentioned in the text in each case. In the NpT ensemble, the temperature and the pressure were controlled by a V-rescale thermostat [45] and a Parrinello-Rahman barostat [46] (isotropically to p = 1 atm), respectively. For the NVT ensemble, the temperature was controlled by V-rescale and no pressure coupling was applied. All bond lengths were constrained with the LINCS algorithm [47]. For water the TIP4P/2005 model was used [48].
The simulation outcomes were analyzed through our own programs, Gromacs and VMD plugins [49,50].

B. Parametrisation of the 2,6-lutidine molecule
We represent the 2,6-lutidine molecule, C 7 H 9 N, by 11 atoms wit the two CH 3 groups treated as single united atoms as shown in Fig. 1, while the hydrogen atoms that are attached to the ring carbons are explicitly included. The GROMOS force field does not provide partial charges on 2,6-lutidine molecule, therefore, we obtained an initial estimation of these charges from quantum chemistry simulations. Then, we varied the values of these charges until two goals were achieved: (i) an agreement with the experimental results for the density ρ and the heat of vaporization ∆H vap for liquid lutidine at temperatures around room temperature, and (ii) the existence of the LCP for the mixture of 2,6-lutidine and water in the experimental range (details of 2,6-lutidine/water simulation will be given in Subsec. II D). More precisely, we gradually scaled the partial charges of 2,6-lutidine molecule by a factor in order to change its Coulomb interaction with water until we obtained the LCP temperature close to experimental value. From experiments, we know that we should have a mixture at T = 280 K and a phase-separated system at T = 320K. Therefore, for each rescaled charge distribution, simulations were carried out at these two temperatures. The appearance of mixed and phaseseparated phases at these two temperatures ensures that the LCP is located somewhere in between. During the rescaling, we also took advantage of the fact that the molecule is symmetric and that the total charge is zero. Moreover, we kept the charges on the hydrogen atoms and CH 3 -groups consistent with the Gromos force field and just vary the remaining 3 parameters.

C. Validation of the lutidine model
Simulations were performed for 80 ns using a 2 fs time step in the NpT ensemble for the bulk system and in the NVT ensemble for the surface tension calculation. Results collected after an equilibration of 10 ns, are presented together with the corresponding experimental data taken from Refs. [51][52][53][54][55] in Table. I. The simulations provide data for the density ρ, the static dielectric constant ǫ and the heat of vaporization ∆H vap in a fair agreement with the experimental ones. The heat of vaporization was calculated as [56] where · denotes time average. U gas and U liquid are potential energies of lutidine in the gas and liquid phases. The gas phase is considered as ideal so U gas contains only interactions within the molecules. The surface tension γ was obtained from simulations of the liquid with two flat surfaces. This system was created by increasing the periodic box size of the equilibrated bulk system by one order of magnitude in one direction. The usage of a 1.2 nm cut-off for the LJ interactions has resulted in a too small surface tension. Therefore, we performed simulations of this system using PME treatment of LJ interactions (LJ-PME).
This gave the surface tension of 32.5 ± 0.2 mN/m in a better agreement with experiment.
With the fractional charges used in the present model, the dipole moment of lutidine is 2.5 D.
The experimental gas phase value is 1.7 D (no experimental value is available for liquid).
Similar differences are encountered for most water models; typically they predict about 30% larger dipole moment than possessed by water in a gas phase, although the difference here is slightly bigger. These differences are caused by the mutual polarization of molecules in the liquid phase. The good agreement between the experimental dielectric constant and the one calculated from fluctuations in the total dipole moment of the entire simulation box is reassuring and indicates that the electrostatic properties of the 2,6 lutidine fluids are properly modelled.
Finally, the heat capacities at constant pressure and constant volume, C p and C V , were calculated from the enthalpy and energy fluctuations in the simulations as The calculated value for C p given in the last row of the Table. I shows a substantial difference compared to the experimental value. The reason for this is that the classical treatment of the lutidine molecules allows too much energy to be taken up by degrees of freedom that in reality are quantum mechanical. A quantum-correction could, however, be added to the classical heat capacity assuming that the relevant degrees of freedom could be approximated as coupled harmonic oscillators. To do this we follow Refs.
[57]- [58] and determine the normal modes of the system from the velocity auto-correlation functions.
The normal mode distribution S(ν) was obtained from the Fourier transform of the velocity correlation functions and the quantum correction to the heat capacity C QM.corr V was obtained as the difference between the heat capacity of a quantum oscillator and a classical oscillator (k B ) integrated over the normal mode distribution [57] C QM.corr with u ≡ hν/k B T being the energy in thermal units.
The quantum corrections to the heat capacities at three temperatures were obtained accord-  drogen atoms in the CH 3 groups. These contributions are negligible. The corrected heat capacity C corr V,p was obtained as the sum of the classical value and the quantum correction In the last two columns of Table. II, the obtained quantum corrected heat capacities C corr p and the experimental ones can be compared.
For most fluids (as well as other condensed matter systems) (C p − C V )/C V is much smaller than 1. For liquid water for instance, this quotient is about 0.01. It is therefore a bit surprising that this quotient is as large as about 0.5 for the present system. As a consistency check, we use the exact thermodynamic relation to calculate C p − C V for lutidine from the molar volume, coefficient of thermal expansion The thermal expansion coefficient and the isothermal compressibility were obtained from the fluctuations in the NpT simulations. The isothermal compressibility was obtained from the volume fluctuations, while the thermal expansion coefficient was obtained from the cross correlations between volume and enthalpy fluctuations. The appropriate equations are discussed in Ref. [59] and   The data from the simulations shown in Table. III are consistent with the data in Table. II.
For comparison, the experimental data for liquid water (taken from [60]) are also shown in the table. The table shows that the main reason for the big difference between the two heat capacities for lutidine is the large thermal expansion coefficient (which is squared in the Eq. 5). It is worth mentioning that the similar big difference between C p and C V occurs in liquid benzene [61] and is reported experimentally in Ref. [62]. Experimental data for benzene are also given in the table.
D. Simulations of the 2,6-lutidine/water mixture In this section we present the simulation results for the 2,6-lutidine/water mixture. In order to simulate a 2,6-lutidine/water mixture, we took as initial configuration a box of the size (L, L, 7L), with L ≈ 3.8 nm, containing N lutidine = 2050 equilibrated bulk lutidine molecules, placed it at the center of the periodic box (L, L, 7L) with L ≈ 5.8 nm and filled (solvated) with N water = 31325 equilibrated water molecules described by TIP4P/2005 model. This configuration corresponds to a lutidine mole fraction x lut = 6.14%, which is close to the experimental value at the LCP. This initial configuration, which is neither a mixed phase nor a two-phases mixture, is shown in the top panel of Fig. 2; this configuration we used for all studied temperatures here. The non-cubic shape of the box with one side much longer than the other two sides has been chosen for two reasons. Firstly, we want "planar interfaces" separating lutidine-rich and lutidine-poor phases for temperatures as close to the critical temperature T c as possible. As discussed in details in Refs. [63][64][65], the larger the size ratio of the rectangular box, the closer one can approach T c keeping the slab structure of the lutidine-rich phase and, hence, the planar interface. If the size of the box in the z direction is not big enough, the slab structure is not stable close to T c . Rather, the lutidinerich phase forms a cylinder or a sphere. Within the simulation box that we have chosen, all volume fractions at all considered temperatures produce a planar interface. Secondly, we want the finite-size effects to be small in order to be able to determine the near-critical properties of a system as accurately as possible [66][67][68][69]. This is achieved by choosing the size of the box to be larger than the bulk correlation length for all studied temperatures.
Moreover, due to the enlargement of the periodic box in the z-direction the two interfaces in the slab structure do not influence each other.
We simulated the 2,6-lutidine/water mixture in the NpT ensemble for various temperatures; the simulation setting is given in Sec. II A. The time evolution of the local mass density of 2,6-lutidine for several temperatures are presented in Fig. 3. The plots show clearly a phase separated system at T = 320 and 380 K and a homogeneous mixture at T = 280K. The actual equilibration time is determined with the time at which not only the  density profiles remain the same within the statistical errors, but also the energies and all hydrogen bonds become stable. Fig. 4 presents the density profiles averaged over different time intervals after the equilibration, for two temperatures. Depending on temperature, stable equilibrium structures were reached after 0.7 − 3.7 µs. The simulations took several months, with average run of 50 (ns/day) for each temperature.
The results of the simulation are reported in the next section.

III. SIMULATION RESULTS
A. Phase behaviours of the 2,6-lutidine/water mixture near the lower critical point show that the two fluids mix at T = 280 K, while they phase separate at the higher temperature T = 380 K. In order to assure that the phase separation is not effected by the initial configuration, the simulations at these two temperatures were redone with different initial configurations. Although we used a mixed phase as the initial configuration for the higher temperature T = 380 K and a phase-separated initial configuration for the lower temperature T = 280 K, the same final configurations as given in Fig. 2 Table. IV. respectively) were reproduced.

(bottom and center
To quantify the phase separation, the mass fraction of 2,6-lutidine w lut (z) has been calculated as a function of z-coordinate from simulations at different temperatures, see Fig. 5.
'Classical' theories for the interface separating two coexisting phases such Cahn and Hilliard [70] or Landau-Ginzburg theory, predict a hyperbolic-tangent shaped interfacial density profiles. This has later been verified in simulations of interfaces [66,[71][72][73][74][75][76][77][78]. Since there are two interfaces in the present system, we fit the density profile to the function with w r lut and w p lut being the mass fractions of 2,6-lutidine in the lutidine-rich and lutidinepoor phases. λ is a measure of the width of the interface (softness of the transition between the two regions) and is proportional to the correlation length ξ, which is defined from the decay of the density-density correlation function. c is half width of the lutidine-rich region and z 0 is center of the lutidine-rich phase. Fits Eq. 7 to the profiles are shown in Fig. 5 as dashed lines. The parameters obtained from the fits are given in Table. IV while the coexistence curve obtained from the data in the table is shown in Fig. 6. In order to compare with experiments [22,24,25] one may need to convert between mass fractions and where w, x indicates mass and mole fractions respectively, while i, j refers to water and 2,6lutidine molecules with molar masses M i,j . The mass and mole fractions of 2,6-lutidine are shown in the left and right vertical axes in Fig. 5 for several temperatures, and by circles and squares in Fig. 6.
When T c is approached from above, λ increases and will eventually not be much smaller than c. Then, the system will be too small to accommodate a saturated lutidine-rich phase.
We note that this starts to occur at T ≈ 320 K. This makes it difficult to determine the lutidine mass fraction w r lut in the lutidine-rich phase. At higher temperatures we could just read off the value of w r lut from the flat part of the mass fraction profiles. The alternative way to determine w r lut is by fitting the parameters in Eq. 7 to simulation data. At high temperatures, this method gives the same values for w r lut as the ones extracted from the flat parts of the mass fraction profiles, whereas for 315 and 318 K it gives the values that are higher than the maxima in Fig. 5. The question is now whether this procedure could be trusted or not. First we note that the fit of the mass fraction profile to the functional form of Eq.7 looks very good. We tested the validity of this procedure further by running  Table. IV, and mole fractions (squares) for the poor and rich phases obtained using Eq. 8.
simulations of a smaller system in which the width of the lutidine-rich region is about one third of the present one, at two temperatures 330 and 380 K. For the smaller system the density profiles do not exhibit flat regions corresponding to the lutidine-rich phases at these temperatures, unlike the original simulations. By fitting the mass fraction profiles to Eq. 7, we do however obtain the same values for w r lut in the lutidine-rich regions as in the large system at the same temperatures, despite that the maximum of the mass fraction profiles are about 20% smaller compared to the larger system. Although we do know that it eventually may break down close to T c , but we are clearly not close enough for that. Now, we turn to the critical properties of the system. We fully realize that an accurate calculation of critical exponents and T c would require simulations closer to T c , larger systems and a finite-size scaling analysis [79]. This would for the present fairly complicated system call for an unjustifiable amount of computer resources. Therefore, based on our simulations, the calculated exponents are obtained within certain amount of statistical errors.
As mentioned earlier, in order to minimize the finite-size effects we limited our simulations to the range of temperatures further away from T c , for which the bulk correlation length of the mixture is distinctively smaller than the linear dimension of the simulation box (see below) and the order parameter is relatively large. Under such conditions we do not expect the mean field behaviour to occur [80,81]. On the other hand, in this temperature range the corrections to the critical scaling are relevant and therefore we will employ them. The critical  Table. V.
The green diamond symbols represent the experimental data [23]. point and the shape of the coexistence curve were determined by using well established procedures [82][83][84][85][86][87][88]. By employing the Wegner expansion [89], the order parameter (OP) which in this case is the mass fraction difference of lutidine in the rich and poor phases w r lut − w p lut , can be written as where τ = T −Tc Tc . The rectilinear diameter can be similarly expressed as w r lut + w p In Eqs. 9-10 w c is the mass fraction of lutidine at the critical point, A i , i = 0, 1 and B i , i = 0, 1 are non-universal constants, while β, ∆ and Γ are universal critical exponents. For the 3D Ising universality class relevant for the present study, the exponents are approximately β = 0.326, ∆ = 0.5 and Γ = 0.1 [90][91][92][93]. We fitted the obtained w p lut and w r lut from the simulations (Table. IV Table. V are shown in Fig. 7 together with the experimental coexistence curve [23]. The obtained curve exhibits a slight shift to the right as compared to the experimental one. In Fig. 8   as a function of the reduced temperature together with the fit to Eq. 9 with (left panel) and without (right panel) correction-to-scaling term. One can see that the OP is fitted quite well to the power law with the 3D Ising exponent β (after ignoring two data points furthest away from T c ); including the correction-to scaling makes this fitting work also for temperatures further away from T c . As a check for consistency of our estimates, we performed fitting treating the exponent β as a free parameter. This fit provides a value of β that is only slightly different from the 3D Ising exponent (see Table. V).
We also estimated the temperature interval in which the UCP of the mixture is located.
This was done by running simulations (with the similar setting as in the original ones) for a smaller system (L x = 5nm, L y = 5nm, L z = 12nm) at several higher temperatures. The simulations data show a phase-separated mixture at 450 K, while a mixed liquid phase at 510K. This indicates that UCP is located between 450-510 K, in agreement with experiments [22,25,27].
B. The surface tension and the correlation length of 2,6-lutidine/water mixture We computed the surface tension in the phase-separated system (above the LCT) as a function of temperature from the simulations as an integral of the difference between the normal and tangential components of the pressure (stress) tensor across the interface [94]; the result as a function the reduced temperature is plotted in Fig. 9. The surface tension and the correlation length have similar scaling behaviour as the OP, but with different universal exponents [85,[94][95][96]  without correction-to-scaling, B 1 = 0, and β fixed to the Ising value (0.326), while we ignored two temperatures furthest away from T c (inset is in log-log scale).
where C i , i = 0, 1 and λ i , i = 0, 1 are non-universal constants. We fitted the simulation data of γ and λ to these expressions using the value of ∆ fixed to its 3D Ising universality class value. In the Tables. VI-VII, the first rows are fits with the exponents µ and ν fixed to the 3D Ising universality class values, while the second rows show the fits with µ and ν as free fitting parameters. The last rows show fits obtained by imposed T c = 310.5 ± 1.5 as estimated from the coexistence curve fit (   Table. VI). The inset is in log-log scale.
(right) The surface tension versus τ determined by a fit to Eq. 11 without correction-to-scaling, C 1 = 0, and ignoring two temperatures furthest away from T c .
FIG . 10: (left) The thickness of the interface, λ, extracted from simulations (symbols) together with a fit to Eq. 12 (parameters from the first row in Table. VII). The inset is in log-log scale. (right) The thickness of the interface versus τ determined by a fit to Eq. 12 without correction-to-scaling, λ 1 = 0, and ignoring two temperatures furthest away from T c .
C. Interactions between water and 2,6-lutidine in the mixture    water and lutidine becomes stronger upon decreasing temperature. Fig. 11 (right) shows the number of hydrogen bonds between water and 2,6-lutidine molecules per lutidine molecule.
From the figure it is seen that upon decreasing temperature the number of hydrogen bonds between water and 2,6-lutidine molecules increases, in line with the behaviour seen in Fig. 11 (left).

IV. CONCLUSION
In this study, we have proposed an atomistic description of the 2,6-lutidine molecule which we have shown is able to successfully describe bulk 2,6-lutidine liquid. We then   have employed this model together with the TIP4P/2005 water model to study the phase behaviour of 2,6-lutidine/water mixture near the LCP. We conclude that by using these models for molecules it is possible to describe the critical properties of the mixture well.
From the density profiles computed in simulations we have obtained phase diagram of the mixture with the lower critical temperature 310.5 ± 1.5K, which is just a couple of degrees higher than the experimental value [22,24,25,27]. We have found that the UCP is located between 450-510 K, in agreement with experiments [22,25,27]. We have computed the order parameter, the surface tension and the correlation length as a function of temperature.
As expected, the order parameter and the surface tension vanish upon approaching the LCP from above, while the correlation length increases. Moreover, we have found that close to T c the temperature dependence of these quantities is well described by power laws. The calculated exponents deviate less than about 0.02 from those of the 3D Ising universality class [90][91][92][93] to which the studied system belongs. However, the estimated errors [0.02−0.18] are clearly larger than this. A more accurate calculation of the critical exponents and of Tc would require simulations closer to Tc , larger systems and a finite-size scaling analysis [79].
V. ACKNOWLEDGMENTS