Unravelling hydrogen bond network in methanol-propanol mixtures via molecular dynamics simulation and experimental techniques

ABSTRACT Hydrogen-bonded networks in 1-propanol-methanol binary mixture are studied using molecular dynamics simulation along with experimental techniques (FTIR, dielectric spectroscopy and refractive index measurements) for the entire concentration range. The structure of hydrogen-bonded networks is studied through radial distribution function, hydrogen bond statistics and graph theoretical analysis from molecular dynamics. It is observed that the probability of hydrogen bonding is maximum in methanol molecules followed by hydrogen bonding between methanol and 1-propanol molecules and minimum among 1-propanol molecules, as the concentration changes. The graph theory results show the formation of linear chain hydrogen bond networks, with no caged structures, which is validated from the experimental results. Further, the deconvolution of the FTIR OH peak suggests the presence of linear multimers in all concentrations of the binary system. Average degree of the hydrogen-bonded networks from graph theory indicates a complex networks among pure alcohols and dimers between methanol and 1-propanol in their binary mixtures at all concentrations. The negative values of excess permittivity indicate that the components of the mixture interact in a manner such that the net dipole polarisation decreases.


Introduction
Binary liquid mixtures have a variety of applications including catalyst separation [1], chemical analysis [2] and microfluidics [3,4].The interactions between the molecular species can be any of the intermolecular interactions such as hydrogen bonds and dipole-dipole interactions [5].Physico-chemical properties of such systems are highly correlated with the intermolecular interactions in them [5,6].Thus, understanding the nature and strength of hydrogen bonding in binary alcohol mixtures are essential both for predicting and for controlling their properties.For example, the ability to manipulate the hydrogen bonding interactions in these mixtures can be used to control the solubility of certain compounds or to improve the stability of pharmaceutical formulations [7].Additionally, binary alcohol mixtures can serve as model systems for studying the behaviour of more complex systems, such as proteins or DNA, which also rely heavily on hydrogen bonding interactions for their structure and function [8].
Alcohol-alcohol binary liquid mixtures represent a class of binary liquids in which both components are capable of hydrogen bonding [9].The complex interplay of the hydrogen bonding network (HBN) can lead to a variability in binary mixtures' physico-chemical properties.The role of hydrogen bonds is traditionally studied by both experimental (Fourier transform infrared spectroscopy [10][11][12] (FTIR), dielectric spectroscopy [13,14], nuclear magnetic resonance(NMR) [15][16][17] and ultrasonics [18,19]) and computational methods (Quantum mechanical methods such as Density Functional Theory [20], molecular dynamics (MD) simulations and Monte Carlo (MC) simulations [21,22]).Among modelling methods, MD is preferable due to its ability to capture the time evolution of the system and its relative cost-effectiveness [23].MD studies of ethanol-water [24] binary mixture showed the gradual breaking of hydrogen bonds between ethanol molecules with increasing water concentration and the HBN increases between water and ethanol.Similarly, the methanol-water and methanol-acetone system [25] showed the variation of aggregation of methanol as a concentration of water and acetone changes.
Previously, hydrogen-bonded structures such as larger multimers, alcohol trimer, anion-alcohol, ion pair-alcohol and the ion cluster-alcohol complexes were identified in the binary system of 1-butyl-3-methylimidazolium tetrafluoroborate and methanol/ethanol using FTIR [26].Similarly, refractive index studies of isoamyl alcohol with some alkoxyethanols binary systems found the interactions between the solvent molecules to increase with increasing chain length of alkoxyethanols [27].Dielectric spectroscopic studies of n-propyl alcohol-ethylenediamine indicate the presence of strong interactions between n-propyl alcohol and ethylenediamine molecules [28].
In our previous studies of the binary mixtures of methanol-ethanol [29] and ethanol-1-propanol [30], anomalies were found in the methanol-ethanol system at lower concentrations of methanol, while no anomalies were found in the ethanol-1-propanol system.Additionally, MD studies of acetone-alcohol and aniline-alcohol binary mixtures of alcohols with varying chain length (C R H 2 R + 1 -OH, R = 1 to 8) were investigated [31,32].The results show that as the alkyl chain length increases, the changes in HBN for R > 3 are less significant compared to R = 1, 2. As a continuation of our work in hydrogen-bonded liquids, in this paper, the binary mixture methanol (MeOH) and 1-propanol (PrOH) are studied over entire concentration using experimental (FTIR, dielectric, refractive index) and MD simulation.The characteristics of HBs are studied using radial distribution function g(r), hydrogen bond statistics and graph theoretical analysis (GTA).

Computational details
MD simulations are carried out using GROMACS software [23,[33][34][35][36][37][38].MeOH and PrOH are modelled using OPLS/AA force field [39,40].The binary mixtures of MeOH and PrOH are modelled with 1000 molecules packed in both pure form and mixture in the number ratios of 9:1 to 1:9 increasing and decreasing the number of molecules by 100 in each concentrations.The initial coordinates of MeOH and PrOH and the binary mixture are generated using packmol software [41].The steepest descent algorithm is employed for energy minimisation [42].After energy minimisation, the system is annealed from 498 K to 298.15 K, with the temperature being reduced by 20K/step over 2ps.The equilibrated systems are then subjected to isothermal-isobaric (NPT) and canonical (NVT) ensembles for 5 ns each, using the Berendsen barostat [43] and the V-rescale thermostat [44], respectively.To obtain reliable results, the equilibrated systems are then used to run MD production for 25 ns in a microcanonical (NVE) ensemble, with coordinates and velocities being saved at an interval of 1ps.Finally, the VMD software [45] is used to visualise the system.The density versus concentration of MeOH is plotted in Figure 1.Simulated density values are close to the literature values of pure liquids, showing less than 3% error which validates our simulations are in close agreement between the simulated and experimental values (Figure 1).
Using the gmx rdf module and gmx hbond module of GROMACS, respectively, the radial distribution function g(r) and number of HBs are computed for all systems.The NetworkX [46] package in Python is used for GTA.The average degree of the graphs is determined using coordinates obtained from MD simulation at various time frames of the trajectory in intervals of 100 ps.

Materials
MeOH and PrOH of purity 99.5% and 99.9% from FINAR and Sigma Aldrich, respectively, are used as such without any further purification.Binary mixtures are made at 11 different volume concentrations from 0 to 100% in steps of 10% of MeOH in PrOH.

Apparatus and procedure
Refractive index of the binary liquids is measured with abbes refractometer, which is calibrated with distilled water and sodium light is used as the light source.Rhode and Schwarz Z series Vector Network Analyzer with a Dielectric Assessment Kit (DAK) is used to study the dielectric properties in the frequency range 100 MHz to 20 GHz.Calibration of the dielectric probe kit is done with open, short and load (deionised water).FTIR measurements are made using the ATR (Attenuated Total-internal Reflection) technique on a spectrometer of Perkin Elmer 400.The spectral measurements are done in the wavenumber range 4000 cm −1 to 400 cm −1 with a resolution of 1 cm −1 .Origin software is used for the deconvolution of IR spectra.All the experiments are performed at room temperature.

Radial distribution function, g(r), estimates the likelihood of finding hydrogen atoms donor (H Phydrogen of PrOH and H M -hydrogen of MeOH) in reference to acceptor oxygen atoms (O M and O P -oxygen atom in MeOH and PrOH, respectively). The g(r) of O M -H M , O P -H P, O M -H P & O P -H M is
shown in Figure 2. The first peak for all concentrations in Figure 2 demonstrates the dominance of OH -O interaction indicating the hydrogen bond dimers [47].The first solvation peak (solvation maxima) for O-H is found to occur between 0.15 nm and 0.27 nm, with a width of 0.12 nm.The solvation maxima of the O M -H P and O P -H M interaction decrease with an increase in the concentration of MeOH.The g(r) showed an increase in probability in the order of O M -H M > O P -H M > O M -H P > O P -H P with an increase in concentration of MeOH.This suggests that the interaction between the PrOH molecule has a lower probability of forming a hydrogen bond than MeOH-MeOH and MeOH-PrOH hydrogen bonds, and this probability decreases with an increase in the concentration of MeOH (Figures 2 and 3).
The second solvation peak in Figure 2 corresponds to the second hydrogen bond neighbour indicating the formation of chain cluster among alcohols [48].
The  (COD n ) is calculated, and the calculated values are shown in Table 1.The COD n increases as the concentration of MeOH increases, the COD n of MeOH-PrOH also increases as the MeOH concentration increases in a similar way indicating strong HBN is formed between MeOH and PrOH.The more probability of forming hydrogen bond network among MeOH followed by the MeOH-PrOH molecules and minimum PrOH molecules is due to the steric effect.The absence of the peak corresponding to the second neighbour in Figure 3c, compared to 3a and 3b, is due to the steric effect since the probability of PrOH to form hydrogen bonds beyond the first neighbour is less compared to MeOH.As the hydrophobic chain increases, it decreases the probability of forming hydrogen bonds beyond the first neighbour which has been also seen in our previous works of acetone-alcohol and aniline-alcohol binary mixtures [31,32].
P(s) denotes the probability of cluster distribution of oxygen atoms as a function of cluster size (s).P(s) vs. s for the entire concentration range in the MeOH-PrOH mixture is given in Figure S1 of the supporting information document.The cut-off distance used to define the cluster is 3.5 Å between the oxygen atoms of the MeOH, PrOH and between MeOH-PrOH molecules.The criteria to define the cut-off distance are taken from the minima of g(r) in Figure 3. gmx clustsize module in gromacs is used to calculate the hydrogen bond cluster distribution.The cluster distribution plots show various hydrogen bond clusters with smaller clusters having higher distribution.

Hydrogen bond statistics
The average number of strong hydrogen bonds for O M -H M , O M,P -H P,M , O p -H P and the total number of hydrogen bonds for the entire concentration range in the MeOH-PrOH mixture   Within the binary concentrations, the average number of hydrogen bonds formed between MeOH andPrOH mixtures is highest at 50% concentration, with decreasing probabilities at either side of the concentration range.This is because at 50% concentration, the number of molecules of both MeOH and PrOH are the same, providing an optimal probability for hydrogen bond formation between them.The total number of hydrogen bonds in the system increases as the concentration of MeOH increases.

Graph theory analysis
The trajectory of MD simulations is analysed using graph theory to investigate the HBN in each system.The snapshots of the simulation frames depict the HBN, providing a twodimensional map of the network's formation and structure.The nodes of the graph represented oxygen and hydrogen atoms, and edges were drawn between hydrogen-bonded atoms to depict the hydrogen bond.The cut-off distance between oxygen and hydrogen atoms was taken to be the first solvation maxima of the RDF, which is seen to be between 0.15 and 0.27 nm.For each simulated concentration, a corresponding network graph is plotted and the graph theory plots are added.Figure 5(d) shows a graph theory representation of all HBN in 50% MeOH concentration.
The HBN of all MeOH-PrOH and HBN among MeOH, MeOH-PrOH and PrOH for all concentrations of MeOH are given Table S1 of supporting information document.It is evident from both the HBN and from the rdf calculations that dimers and chain clusters dominate the MeOH-PrOH mixture.The HBN shows that most of the molecules have two hydrogen-bonded neighbours, which is in accordance with the COD n .
Calculating the average degree from this 2-D representation allows to quantify the degree of networking in these graphs.The average degree is calculated using the relation The average degree of HBN calculated for MeOH-MeOH, MeOH-PrOH and PrOH-PrOH is plotted in Figure 6.MeOH has the maximum average degree in the pure mixtures, the average

Dielectric properties
Dielectric investigations are sensitive to conformational changes that take place in binary systems of liquids because they evaluate changes in molecule dipoles.The dielectric permittivity (ε') and loss (ε'') of the binary mixture are given in Figure 7(a, b), respectively.The dielectric loss peaks are observed to shift to higher frequencies as the concentration of methanol increases, which is a sign of a shorter relaxation period as methanol concentration increases.
The experimental values of complex permittivity are fitted in the following Debye equation: where ε ∞ , ε s and τ are the permittivity at high frequency, the static dielectric constant, and the relaxation time constant, respectively.The values of these parameters obtained through fitting is given in Table 2.The square of the refractive index (from the values of refractive index in Figure 8) is taken as the permittivity at high frequency.The values of static permittivity of MeOH and PrOH are found to be in agreement with the literature values [49,50].The excess permittivity is calculated using the relationship below: where X is the mole fraction, the subscript m, A and B represent the mixture and the pure components of the mixture, respectively.The values of (ε s ) E give information about the type of interaction occurring among the molecules in the binary [51] which is tabulated in Table 3.
In the MeOH-PrOH binary mixture, the negative values of excess permittivity indicate that the components of the binary are interacting through hydrogen bonding in such a way that the net dipole moment decreases.
Kirkwood correlation factor provides information about the orientation of electric dipoles in polar liquids.The Kirkwood equation for the mixture is expressed by the following (Equation 4): where X A and X B are the mole fractions of liquids A and B in the mixture, respectively, ρ A and ρ B are the densities of liquid A and liquid B, and μ A and μ B are the gaseous-phase dipole moments of the liquids A and B, respectively.g eff is the effective Kirkwood correlation factor for the mixture, which is an index of solute-solvent interactions.An ideal non-interacting mixture will have a g eff value of unity, and any deviation from unity is the indication of heterogeneous interaction between the components in the mixture [52].
The values of dipole moment and density of MeOH and PrOH from literature [51] used in the calculation are tabulated in Table 4.The value of g eff is high for PrOH while compared to MeOH.This indicates that PrOH favours parallel orientation of dipoles among themselves as compared to MeOH.This is expected since the hydrogen bonding between MeOH molecules is stronger and their networks are larger, the average dipole moment of the hydrogen-bonded species will be more random, in comparison with the longer PrOH.The parameter g eff decreases as the concentration of MeOH in the binary mixture increases as is shown in Table 3.This trend leads to the inference that heterogeneous interaction between the compounds occurs forming hydrogen-bonded multimers with anti-parallel orientation of the electric dipoles [53].

Refractive index
The refractive index of the MeOH-ProH binary mixture is given in Figure 8.The refractive index readings of pure MeOH and pure PrOH from the literature were found to be 1.33057 and 1.3840, respectively [54], which are close to the experimental values.A linear decrease in the refractive index with the increase in the concentration of MeOH in the solution is seen, indicating the possible absence of caged (closed) hydrogen-bonded structures [55].

FTIR
The FTIR spectrum for the entire concentration of the binary mixture is given in Figure 9. Alcohols self-associate through hydrogen bonding to form clusters which results in a broad absorption in the 3500-3200 cm −1 range [56,57].Due to the severe overlap of OH bands, this broad OH band can be deconvoluted to find different hydrogen-bonded structures in the binary [29].Following this, the OH band of the FTIR spectrum obtained is further deconvoluted which is given in Figure 10 and the three deconvoluted peaks obtained coincide with the peak positions of tetramers, pentamers and hexamers from literature [30,58].Structures of these are also observed in the MD simulations.The figures of the found multimers and peak positions of the same are given in Figure 11.Along with these, several other structures such as dimers, trimers and linear multimers chains are observed via MD simulations (Figures 9,10,11).Hydrogen numbers calculated from MD simulations show that the maximum number of heteroclusters occur at equimolar concentration.Previous studies also indicate that the largest deviation from the excess absorbance spectra for the binary mixtures of methanol with aliphatic alcohols is observed at equimolar concentration and this is attributed to the involvement of the 50% of the molecules in heterocluster formation [20].Thus, the reason for the red-shift seen around equimolar concentration can be attributed to the high probability to form heteroclusters in this region.

Conclusion
The binary mixtures of methanol with 1-propanol are studied over the entire concentration using molecular dynamics simulation and experimental techniques.The radial distribution function shows with an increase in methanol concentration, the g(r) indicates an increase in probability of hydrogen bond formation in the following order: O M -H M > O P -H M > O M -H P > O P -H P .The coordination number, hydrogen bond statistics and average degree calculation also show that methanol forms a more complex hydrogen-bonded network than 1-propanol.The network graph plotted using the graph theoretical approach shows a similar linear chain of hydrogen bond network in maximum concentration of both methanol and 1-propanol having average degree maximum for methanol.The deconvoluted FTIR peaks confirm the presence of tetramers, pentamers and hexamers, which is also validated by MD simulations.The linear decrease of the refractive index values for the entire concentration ranges is indicative of the absence of any caged structures, which is also evident from hydrogen bond network plots.The negative values of excess static permittivity indicate that the components of the mixture interact through hydrogen bonding such that the net effective dipole polarisation decreases.
g(r) of O M -O M , O P -O P and O M -O P is shown in Figure 3.The g(r) of O-O interactions show that the solvation maxima begins at 0.24 nm and ends at 0.35 nm, with a width of 0.11 nm.The RDF analysis indicates that the probability of interaction follows the order of O M -O M > O M -O P > O P -O P similar to O-H interactions.From Figure 3, coordination number

Figure 1 .
Figure 1.(a) simulated density vs concentration of MeOH (b) molar volume vs concentration of MeOH.

Figure 2 .
Figure 2. RDF of (a) O M -H M , (b) O M -H P , (C) O P -H M , (d) O P -H P .

Figure 3 .
Figure 3. RDF of (a) O M -O M , (b) O M -O P , (C) O P -O P .

Figure 6 .
Figure 6.Average degree versus concentration of MeOH.

Figure 7 .
Figure 7. (a) variation of dielectric permittivity with frequency for all concentrations of the methanol-propanol binary mixture, (b) variation of dielectric loss with frequency for all concentrations of the methanol-propanol binary mixture.

Figure 11 .
Figure 11.Variation of peak position of deconvoluted peaks with concentration of MeOH (%).

Table 2 .
Dielectric parameters for the binary mixture of MeOH-PrOH.

Table 4 .
Molecular parameters for MeOH and PrOH.