Structural and Dynamic Properties of Hydrophobic Eutectic Solvents Based on Menthol and Fatty Acids: a Molecular Dynamics Simulation Study

The structural and dynamical properties of the binary mixture of Menthol (MEN) and Fatty acids (FAs) were investigated using molecular dynamics simulations. We focused on the relationship between the structural and dynamical properties of the eutectic mixtures of MEN and FAs with different molar percentages of FAs. Structural properties of the eutectic mixtures were characterized by calculating the combined distribution functions(CDFs), the radial distribution functions (RDFs), and the angular distribution functions (ADFs), and the Hydrogen bonding network between species and Spatial distribution functions (SDF). Further interaction between menthol and Caprylic acid molecules was confirmed by the results of these analyzes. Also, the transport properties of the mixtures were investigated by using the mean square displacement (MSD) of the centers of mass of the species, self-diffusion coefficients and vector reorientation dynamics (VRD) of bonds. The simulation results indicated that intermolecular interactions have a significant effect on the dynamic properties of species.

[ 1 ] This natural product has been used for medicinal purposes, cosmetics, as a flavor and, as an intermediary in the production of widely used industrial solvents (deep eutectic solvents) [2]. The deep eutectic solvents (DESs) are made by strong interactions between the Hydrogen bond donor (HBA) and Hydrogen bond acceptor (HBD) [3]. Eutectic mixtures have a melting point much lower than their individual components. DESs were first introduced by Abbot et al. (2003) for a mixture of choline chloride and urea [4]. It has been found that DESs and ILs have similar physical characteristics such as viscosity, conductivity, and surface tension [5]. DESs are a good alternative to ILs due to their biodegradability, lower toxicity, and reasonable price [6]. The various applications of eutectic solvents in chemical and pharmaceuticals make the use of DESs very appealing. So far, the most widely proposed solvents in the open literature are hydrophilic eutectic solvents thus are unstable in water. In 2015, Hydrophobic DES was receiving a great amount of attention due to the stable hydrophobic phase in contact with water and inexpensive and natural components that follows the principles of green chemistry [5]. Reports showed that the hydrophobic DESs based on ammonium salts and Decanoic acid were able to recovery of volatile fatty acids from aqueous solutions [7]. Marrucho et al. have reported several Hydrophobic deep eutectic solvents composed of natural, nontoxic compounds, i.e., terpenes and fatty acids [8]. Eutectic mixtures of menthol and saturated fatty acids derived from coconut oil and palm kernel oil [9]. Some DESs have found very interesting applications in the extraction of biomolecules and indium from aqueous solutions [8]. Because eutectic solvents have been widely used in the isolation of pollutants from aqueous solutions, predicting thermos-physical and structural properties is particularly important. MD simulation may provide molecular insight into the eutectic mixtures for medicinal chemists. In the present study, we were interested in focusing on the effect of the percentage composition of species and temperature on the dynamic and structural properties of the eutectic solvents based on menthol and fatty acid. Binary mixtures of Caprylic Acid: Menthol, Decanoic acid: Menthol, and Lauric Acid: Menthol were prepared. Investigation of structural and dynamic properties of eutectic solvents can provide a valuable perspective in the design of these solvents.

COMPUTATIONAL METHODS
The NAMD 2.13 software package was used to run all MD simulations [10]. Three fatty acids, Caprylic acid (CAP), Decanoic acid (DEC), and Lauric acid (LUA), and menthol, were selected. Packmol software package was used for the preparation of the eutectic mixtures with the different molar percentages of FAs [11]. The information about the simulated boxes and their different combinations (HBD and HBA) with abbreviations are listed in Table 1. The initial simulation box included with the mole percent of FAs at 20.0 % and 80.0 % of MEN was provided using distributed randomly. In the next step, the number of FAS molecules was increased for the preparation of the mixtures with the mole percent of FAs at 70.0 %. First, the initial configurations were minimized to remove bad contacts followed. The binary mixtures with the most intermolecular interaction between HBD and HBA were selected and their temperature was increased slowly from 0 to 353 K and then decreased back to 323 K and 300 K with a temperature gradient of 1 K/FS. The Nose-Hoover thermostat [12] and the Parrinello-Rahman pressure coupling [13] were used to maintain a constant temperature and pressure, respectively, during the simulation. Finally, the binary mixtures were equilibrated for 50 ns in the NPT ensemble with a time-step of 1 fs. The particles mesh Ewald (PME) method and the Lennard-Jones (LJ) 6-12 potential with a cut-off distance of 12 A˚ were applied for long-range and short-range interactions, respectively.

Radial Distribution Functions (RDFs)
Radial distribution functions (RDFs) were computed to characterize the various interactions  (green ellipse), and CT12 (red ellipse).

Coordination number
The coordination number is defined as the number of nearest neighbours of a central molecule in the mixture [14]. Generally, the CN can be gained by integrating the surface region surrounded at the first solvation shell (rfs). The CNs are calculated by Equation (1).
, where and represent the number density and the first minimum in the RDFs [15].
Coordination numbers of the molecules around each other in the binary mixtures are listed in Table S1.

Hydrogen bond analysis
In  Table S2).
The H-bond number between FAs and MEN was reduced in the binary mixtures with a higher percentage of FAs. Most likely, Decreasing the average number of hydrogen bonds is due to the distribution of the number of hydrogen bonds between the carboxyl groups of FAs (see Fig.   5a). The hydrogen-bonding network is directly affected by the temperature of the mixture.
Investigation of the hydrogen bonding network in the binary mixture at the temperature range of 300 to 353 also showed that the H-bond number was slightly decreased with raising the temperature (see Fig. 5b). The total time that a hydrogen-bond network has established between HBD and HBA is described using the percent occupancy [18]. This analysis has depicted the stability and persistence of hydrogen bonds between the possible HBD and HBA pairs [19].
The stabilization of the hydrogen bonding network between the different species is provided in

3.3.Spatial Distribution Functions (SDFs)
For further visualization of the local structure of the eutectic mixture, spatial distribution functions (SDFs) were analyzed [20]. SDFs can present the desirable sites for intermolecular hydrogen-bond interactions. Three-dimensional spatial distribution functions of different species around a reference molecule were plotted in Figure 7

Tracing [MEN] [FAs] H-bonds Using Combined Distribution Function (CDF)
The distance criteria of H-bond in the binary mixtures was determined from the first minimum of the peaks of the corresponding RDFs. However, the distance criteria may not be adequate to capture the presence of a strong H-bond in the binary mixtures [16]. Therefore, the combined radial/angular distribution functions were obtained for all possible pairs. The composed CDFs of the RDF between the hydroxyl hydrogen (HA) site of MEN and the OA atom of the FAs as the X-axis, the Angular Distribution Function (ADF) between two vectors (R1 and R2) as the Y-axis are shown in Figure 8a.
H-bonds (weak or strong) were found at the angles range between 130 and150° and distances of less than 2 Å. Figure 8a shows that there is more menthol around CAP than the other acids in the binary mixtures (see Fig.8b and Fig.8c)

Non-bonded interaction energy
The non-bonded interaction energy is an important analysis for understanding microscopic properties in MD simulations. The non-bonded energy includes two terms, short-range van der Waals (vdW) interactions and long-range electrostatic interactions, which are calculated Lennard-Jones (12 -6) function and Coulombic equation, respectively [22]. The vdW and Coul interactions are described by the following equations: The values of the non-bonded energies are reported in Table S3. It could be seen from Table   S3 temperatures from 300 to 353 K (see Fig. 10). According to Figure S1, increasing the length of the alkyl chain of FAs may prevent further interaction between the HBA and HBD.

Dynamical and transport properties
The mean-squared displacements (MSDs) of species were measured to describe the migration of species [23]. The mean square displacements (MSDs) of species as a function of time were plotted for binary mixtures of FAs and MEN in the same ratio in Figure 11a-11b. The low slope of the MSD curve of menthol molecules indicates a strong interaction between the two components in the binary mixture of menthol and Caprylic acid compared to the mixture of MDA. The self-diffusivity coefficient is the macroscopic dynamical property which is calculated from MD simulation [24]. The self-diffusion coefficients of species were obtained from the slopes of the lines fitted to the MSD curves using the Einstein relation: Where, 〈∆| ( )| 2 〉 is the mean-square displacement (MSD) of species i. To ensure the measured values of self-diffusion coefficients, the diffusive regime was identified and sampling was performed in this regime [25]. The parameter was used to determine the location of the diffusive regime in the binary mixtures [26]. This parameter, β is calculated according to: The reduced viscosity at higher temperatures has led to more rapid migration of species in the diffusive regime. The calculated self-diffusion coefficient values from the diffusive regime are collected in Table 2

Vector Reorientation Dynamics (VRD)
The vector reorientation dynamics is a common quantity that can be determined by both experimental techniques [27] and molecular dynamics simulation methods. VRD is computed from equation (4), , where ⃗ ( ) and ⃗ ( + ) are a vector selected on a molecule in step t and at the later time t + τ, respectively [28]. The reorientation dynamic of the bonds vector of the MEN molecules is displayed in Fig.12a.

Thermo-physical Properties Analysis
The shear viscosity η can be obtained from the xy elements of the molecular pressure tensor using the Green-Kubo expression [29]. η is given by , where V and T are stands for the molar volume and the temperature, respectively. KB is the Boltzmann constant [30].
Here, p i x and indicate the momentum component and the interaction potential, respectively. x i is a component of the radius vector [31]. The measured viscosity values are collected in Table 3.

CONCLUSION
In the recent work, the correlation between the structural and dynamical properties was