Microsolvation of phenol in water: structures, hydration free energy and enthalpy

ABSTRACT In this work, we have studied the microsolvation of phenol in water. We started by identifying initial configurations of phenol-water clusters using classical molecular dynamics. The configurations are optimised at the ωB97XD/aug-cc-pVDZ level of theory. To understand the interaction between phenol and the solvating water molecules, we performed a quantum theory of atoms in molecules (QTAIM) analysis. The results show that the structures of phenol-water clusters are similar to those of neutral water clusters. The QTAIM analysis shows that the structures of phenol-water clusters are stabilised by strong OH··· O hydrogen bondings, weak CH··· O hydrogen bondings, and OH bonding interactions. The located structures of phenol-water clusters have been used to calculate the absolute hydration free energy and enthalpy of phenol for temperatures between 20 and 400 K. The hydration energies are calculated using the cluster continuum solvation model. It has been found that the explicit solvation has negligible effects on the hydration free energy and enthalpy of phenol. Furthermore, the hydration free energy of phenol is found to be linearly varying with increasing temperature, while the hydration enthalpy is found to be temperature independent. The estimated hydration free energy of phenol is slightly underestimated as compared to a previously reported experimental estimate.


Introduction
Phenol is an organic aromatic compound with a hydroxyl OH group soluble in water. It is an important compound used in industry for numerous purposes. Phenol and its derivatives are used, for example, in the synthesis of herbicides (used in agriculture), and several pharmaceutical drugs. Thus, understanding its solvation in water, especially its close interactions with solvent water molecules, is of primary importance. That is the reason why the solvation of phenol in water or the phenolwater cluster has received significant attention in the cluster community. As outlined in the next paragraph, phenol-water clusters have been investigated from several perspectives.
The first perspective was the understanding of the hydrogen bond networks of PhOH(H 2 O) n . The study of hydrogen bond networks of the phenol-water clusters, PhOH(H 2 O) n , has been reported by several authors [1][2][3][4][5][6][7][8][9][10][11]. Early investigations of the phenol-water clusters have been reported by Watanabe and Iwata [1] at the HF/6-31G(d) level of theory. The authors have reported different configurations of the phenol-water clusters, PhOH(H 2 O) n , for n = 1−5. Clustering energies, clusters free energies, binding energies, as well as vibration frequencies analysis have been also reported for the studied clusters [1]. In addition, Benoit and Clary [2] reported a similar study using an improved computational level of theory, B3LYP/6-311++G(d,p). Guedes et al. [3] used three different functionals of DFT (density functional theory) to study the solvated phenol clusters, for n = 1−6. They recommended the B3LYP and B3PW91 hybrid functionals over the BLYP DFT functional, which underestimates the PhO-H bond dissociation energy. Parthasarathi et al. [8] reported the investigation of PhOH(H 2 O) 1−3 using ab-initio methods and DFT functionals. The authors have also investigated the hydrogen bonding using the quatum theory of atoms in molecules (QTAIM) analysis. Beside ab-initio and DFT functionals, molecular dynamics simulations have been also used to study the hydrogen bond networks of the solvated phenol in water [9,11].
The second perspective was the understanding of the vibrational features of the solvated phenol in water. Several authors have reported the electronic, and the infrared or vibrational spectroscopy of PhOH(H 2 O) n for different cluster sizes [12][13][14][15][16][17][18]. Early experimental spectroscopic study of the phenol-water clusters was reported by Fuke and Kaya [12]. The authors reported the p − p * electronic absorption of PhOH (H 2 O) 0−3 . Janzen et al. [14] reported the experimental vibrational spectroscopy of PhOH(H 2 O) 7,8 in the OH-stretching region. They found that the vibrational spectra are the combination of features from different isomers. Their ab-initio calculations showed that the structures of the phenol-water clusters are very similar to those of neutral water clusters [14]. Largesized solvated phenol clusters, PhOH(H 2 O) n , 10 ≤ n ≤ 50, have been reported by Mizuse, Hamashima and Fujii [15,16], experimentally. The authors measured the IR spectra of the clusters in the OH-stretching region. The results show that the band associated with the free OH stretching is blue-shifted with the increasing cluster size n [15]. Recently, Katada and Fujii [18] reported the experimental IR spectroscopy of the protonated phenol-water clusters for n = 1−5.
The third perspective was the evaluation of the solvation free energy and the solvation enthalpy of phenol in water. Few authors have reported the study of the solvation free energy of phenol in water using different approaches [19][20][21][22][23]. Hydration free energies (i.e. solvation free energies in water) of several molecules have been evaluated by Gallicchio et al. [19] based on surface generalised Born (SGB) continuum dielectric electrostatic model. The model used by the authors is termed as experiment, and the evaluated hydration free energy of phenol is -6.6 kcal/mol [19]. Reddy and Erion [20] have reported the relative solvation free energy of phenol→benzene using QM/MM (quantum mechanics and molecular mechanics). The relative solvation free energy of phenol→benzene is evaluated to be 4.9 kcal/mol. In addition, hydration energies of phenol as well as other molecules have been evaluated using different empirical solvation models by Sharma and Kaminski [21]. The hydration energy of phenol is evaluated to be −6.6, −6.7 and −5.6 kcal/mol using Fuzzy-border, Poisson-Boltzmann, and Generalised Born models, respectively [21]. Besides, enthalpies of solvation of phenol and substituted phenols in acetonitrile, tetrahydrofuran, and 1,4-dioxane have been reported by Nagrimanov et al. [22].
Observation of the literature reveals that the study of the structures and the hydrogen bond networks of PhOH(H 2 O) n has been limited to small-sized clusters. Moreover, even for the small-sized clusters, no evidence of global minimum potential energy surfaces (PESs) exploration has been reported. It is worth noting that the accuracy of any study on the solvated phenol clusters depends on the accurate identification of all possible configurations of the clusters. Therefore, it becomes necessary to revisit the study of PhOH(H 2 O) n clusters. Thus, we explored the PESs of the clusters of PhOH(H 2 O) n , for n = 1−12, using classical molecular dynamics. The obtained configurations are then optimised using a dispersion-corrected DFT functional, ωB97XD, associated with the aug-cc-pVDZ basis set. In addition, QTAIM analysis has been performed on the most stable configurations of the phenol-water clusters to identify the interactions between phenol and water molecules. Thermodynamics properties and absolute solvation energies are also provided for temperatures ranging from 20 to 400 K.

Sampling of configurations
To achieve the objectives of this work, we need to identify lowlying energy structures of the solvated phenol clusters, PhOH (H 2 O) n , n = 1−12. Thus, we started by sampling possible structures of different clusters. For each cluster size, initial configurations are generated using the ABCluster code of Zhang and Dolg [24,25]. ABCluster samples all possible configurations and classifies them from the most stable to the least stable configuration based on a classical energy. The classical energy used by ABCluster is constituted of Lennard-Jones potential as well as electrostatic potential. The parameters used are based on CHARMM's force field [26]. Details on how the configurations are generated can be found in our recent works [27][28][29][30][31]. In addition, the reader is advised to read the original papers of Zhang and Dolg [24,25] for more details on ABCluster.

Solvation free energy and enthalpy
In this work, solvation free energy and solvation enthalpy refer to the absolute solvation free energy and the absolute solvation enthalpy of phenol, respectively. The solvation free energy and the solvation enthalpy are calculated in this work using the cluster continuum solvation model (CCM). In the CCM approach, the phenol is explicitly solvated with few water molecules, while the remaining solvent is represented by a dielectric continuum medium. The phenol solvation can be represented by Equation (1).
Then, the absolute solvation free energy and the absolute solvation enthalpy of phenol can be calculated using Equations (2) and (3), respectively.
where the subscript s stands for the water solvent and the subscript g stands for the gas phase. Thus, the solvation free energy and the solvation enthalpy of phenol in water are obtained when the calculated values of DG s (Phenol) n and DH s (Phenol) n are not varying with increasing cluster size n (which corresponds to the convergence of the solvation free energies). The CCM approach can be summarised by the schematic representation of Figure 1. Therefore, it comes out from the above equations and scheme that the calculation of the solvation free energy and enthalpy of phenol depend on the determination of the structures of the solvated phenol clusters, Phenol(H 2 O) n , as well as the structures of neutral water clusters. The structures of the solvated phenol clusters are generated and optimised in this work, while the structures of neutral water clusters are retrieved from our previous works [32,33]. It should be noted that the cluster continuum solvation model has been successfully used in previous works to determine the absolute solvation free energies of the proton and other ions [34][35][36][37][38][39][40][41][42][43].
Although not explicitly shown, Equations (2) and (3) are temperature dependent. To calculate the free energy and the enthalpy of Phenol(H 2 O) n and (H 2 O) n , we use the free energies of different possible isomers, G k (T), weighted by their relative probabilities, W k (T). The relative probabilities are and Thus, the contribution of the isomers to the final free energy of a given cluster is dictated by its Boltzmann weight W k (T). As will be seen later, most of the isomers do not contribute to the cluster's population or they contribute with negligible probability. It is worth mentioning that the free energies G k (T) of different isomers as well as the weight W k (T) are calculated using the Tempo code of Fifen and coworkers [44]. The program Tempo computes the free energies G k (T) for different values of temperature using the canonical ensemble. To calculate G k (T), Tempo uses the electronic energies and the harmonic frequencies from the Gaussian output files. With the calculated free energies G k (T), Tempo computes the weight W k (T) using the Boltzmann distribution of Equation (5).

Computational details
After sampling all the configurations of the solvated phenol clusters using ABCluster, the generated configurations have been fully optimised at the ωB97XD/aug-cc-pVDZ level of theory. To ensure that we have located the most stable structures, frequencies have been calculated at the same level of theory. It has been checked that no imaginary frequency is found for the reported structures. Besides, the frequencies are also used to calculate the thermodynamic properties (free energies and enthalpies) of the solvated phenol clusters. The optimisations and the calculations of frequencies have been performed using the Gaussian 16 suite of program [45]. For accurate optimisations, we used the tight option and ultrafine grid for accurate integrals. All optimisations are performed in the implicit solvent phase using the polarisable continuum solvation model [46]. Single-point calculations have been performed at the SCS-MP2/aug-cc-pVTZ level of theory. SCS-MP2 [47] calculations are performed using Orca computational chemistry program [48]. The single-point calculations are performed in the implicit solvation using the SMD solvation model [49]. The quantum theory of atoms in molecules (QTAIM) analysis has been performed using the AIMAll program [50]. QTAIM analysis has been performed to understand the non-covalent bondings between the water molecules and phenol.

Results and discussions
As stated in the methodology section, to compute the solvation free energy and enthalpy of phenol we need the structures of the solvated phenol clusters, PhOH(H 2 O) n . Thus, we start this section by presenting the located structures of the solvated phenol clusters, PhOH(H 2 O) n , n = 1−12 (see Section 3.1). After presenting the structures, we performed a QTAIM analysis on the most stable structures to understand the nature of non-covalent bondings between the solvating water molecules and phenol (see Section 3.2). The cluster continuum solvation approach presented in the previous section is temperaturedependent. Among the located isomers, only those having non-negligible probabilities can contribute to the population of the cluster, and therefore only those isomers can contribute to the calculation of the solvation free energy and enthalpy. Consequently, we presented in Section 3.3 the probability of all the isomers, for different cluster sizes, as a function of temperature. Finally, we reported in Section 3.4 the calculated absolute solvation free energy and absolute solvation enthalpy of phenol in water. These solvation energies are reported at room temperature as well as at temperatures ranging from 20 to 400 K.

Microsolvation of phenol
The optimised structures of the solvated phenol monomer and dimer are reported in Figure 2.
After generating the initial configurations using ABCluster, optimisation has been performed at the ωB97XD/aug-cc-pVDZ level of theory. For the solvated phenol monomer, PhOH(H 2 O) 1 , only one configuration is found to be stable (see PW1 in Figure 2). In PW1, the water molecule establishes one strong OH··· O hydrogen bond and one weak CH··· O hydrogen bond with the phenol molecule. In addition, the water molecule acts as a hydrogen bond acceptor, while the OH group of phenol is the hydrogen bond donor (see Figures 2 and 8). For the case of the solvated phenol dimer, PhOH(H-2 O) 2 , eight initial configurations have been optimised; among which, five are found to be different from each other covering an energy cutoff of 2.7 kcal/mol (see Figure 2). The most stable isomer of the solvated phenol dimer is PW2_1. In PW2_1, the two water molecules establish one strong OH··· O hydrogen bond and two weak CH··· O hydrogen bonds with the phenol molecule. In addition, the two water molecules are linked by a strong OH··· O hydrogen bond (see Figures 2,and 8). The second and third most stable isomers, PW2_2 and PW2_3, lie 0.9 and 1.0 kcal/mol above the global Figure 2. (Colour online) Structures of the solvated phenol with one and two explicit water molecules, as optimised at the ωB97XD/aug-cc-pVDZ level of theory. Numbers in blue color represent the zero-point corrected relative electronic energy of isomers in kcal/mol. minimum energy structure. In these two isomers, the two water molecules and the OH group of phenol form a cyclic OH··· O hydrogen bondings network (see Figure 2).
It is should be noted that previous works have also reported the PW1 isomer as the most stable configuration of the phenol-water monomer [1,3,5,8]. Watanabe and Iwata [1] have also reported another stable structure at the HF/6-31G level of theory. That isomer was reported to lie about 3.0 kcal/mol above PW1, highlighting its less stability [1]. As far as the phenol-water dimer is concerned, all previous works have reported the cyclic isomers, PW2_2 and PW2_3, to be the most stable isomer. This result is in contradiction with our finding which shows that the isomer PW2_1 is the most stable. It is worth mentioning that the isomer PW2_1 is reported here for the first time. There are two possible reasons that can explain the difference in the located most stable configuration of the phenol-water dimer. First, the structures reported in this work are optimised in the solvent phase, while previous works reported gas phase structures. Second, no global optimisation has been performed previously, which could have led to missing isomers (such as PW2_1) on the PES of the phenol-water dimer. The isomers PW2_4 and PW4_5 are located in this work for the first time.
Fifteen initial geometries of the solvated phenol trimer, PhOH(H 2 O) 3 , have been optimised. After optimisation, 11 structures have been located to be a different one from another. The retained structures are reported in Figure 3 covering an energy cutoff of 4.2 kcal/mol. The results show that there are three iso-energetically most stable isomers of the solvated phenol trimer, PW3_1, PW3_2, and PW3_3 (see Figure 3). For the three isomers, the three water molecules and the phenol OH group form a cyclic OH··· O hydrogen bondings network. The difference between the isomers lies in the orientation of the free OH vis-à-vis the OH··· O plane. This follows the same configurations as the most stable structures of neutral water tetramer [32,51,52]. It has been noted that apart from the three most stable structures, the other isomers establish only one OH··· O hydrogen bond with the phenol OH group. We concluded that the stability of the solvated phenol trimers is affected by the number OH··· O hydrogen bonds established with phenol. The isomers in which the phenol OH group is both hydrogen bond donor and acceptor are found to be more stable than the isomer in which the phenol OH group establishes only one hydrogen bond, although sometimes with additional OH· · · p or CH··· O interactions.
All previous works reported the isomer PW3_1 to be the most stable configuration of the phenol-water trimer, in agreement with our result [1,3,5,8]. However, two more similar isomers, PW3_2 and PW3_3, have been reported additionally in this work. The three isomers, PW3_1, PW3_2 and PW3_3 have the same hydrogen bond networks as the most stable isomers of the neutral water tetramer [33]. furthermore, most of the configurations of Figure 3 are located in this work for the first time.
Regarding the solvated phenol tetramer, 16 stable configurations are located on its PES. The located isomers have their relative energies spanning from 0.0 to 2.4 kcal/mol (see Figure 4). The most stable structure, PW4_1, has a cyclic OH configuration. In addition to the five strong OH··· O hydrogen bondings, PW4_1 has two weak CH··· O hydrogen bondings and one OH· · · p interaction (see Figure 8). The high stability of PW4_1 is attributed to the fact that it establishes extra bonding interactions in addition to the cyclic OH··· O hydrogen bondings. An isomer with cyclic OH configuration has been located by previous authors to be the global minimum structure of the phenol-water tetramer [1,3,5]. In that isomer, the four water molecules as well as the phenol OH group act as a proton donor and proton acceptor without further interaction with the phenyl group. After optimisation at ωB97XD/aug-cc-pVDZ level of theory, the isomer converged to PW4_4. The difference between PW4_4 and the isomer located in previous works is that in PW4_4, the water molecules establish further interactions with the phenyl group through CH··· O and OH· · · p interactions. Comparison with previous works shows that all the isomers of Figure 4 are reported in this work for the first time. It is worth mentioning that the location of new isomers in this work has been possible following global optimisation using classical molecular dynamics as implemented in ABCluster [24,25]. This essential step in the exploration of the PESs of clusters has been skipped in previous work, which also explains the difference between current work and previous results. As can be seen in Figure 4, the phenol OH group of low-lying energy structures (from PW4_1 to PW4_6) acts as a proton donor and proton acceptor, establishing therefore two strong hydrogen bonds with the solvating water molecules. However, in less stable structures, the phenol OH group is either a proton donor or proton acceptor, establishing only one hydration bond (see PW4_14, PW4_15 and PW4_16 in Figure 4).
Initially, 20 geometries of the phenol-water hexamer have been optimised at the ωB97XD/aug-cc-pVDZ level of theory. After optimisation, 16 structures are found to be a stable and different one from another. The located structures, covering an energy cutoff of 2.5 kcal/mol, are reported in Figure 5. The most stable structure, PW6_1, exhibits a cage like configuration similar to the case of water heptamer. In addition to strong OH··· O hydrogen bondings, PW6_1 has further bonding interactions resulting from the interaction of the water molecules with the phenyl group (see Figure 8). Exploration of the literature shows that few authors investigated the hydrogen bond networks of phenol-water clusters larger than the pentamer. Guedes et al. [3] reported an isomer with a double cyclic OH configuration as the most stable structure of the phenol-water hexamer. It should be noted that in the configuration located by Guedes et al. [3], the water molecules do not interact with the phenyl group. This is due to the lack of correlation (dispersion interactions) in the HF method used by the authors. On the other hand, after taking into account the dispersion corrections in this work (through the ωB97XD DFT functional), all the located isomers interact with the phenyl group in addition to the OH··· O hydrogen bondings (see Figure 5). In addition, examination of the structures shows that the hydrogen bond network of the phenolwater hexamer is similar to that of neutral water heptamer. However, the phenol-water hexamer structures establish further interactions with the phenyl group, generating lesser free OH stretching.
The number of possible isomers on the PESs of the phenolwater clusters increases with the cluster size. For the phenolwater octamer, among the 25 optimised configurations, 21 are retained and reported in Figure 6. The isomers are located within the cutoff energy of 3.9 kcal/mol. The most stable configuration of the octamer, PW8_1, has a folded cage structure. The folding originates from the interaction of one water molecule with the phenyl group. Thus, in addition to the OH··· O hydrogen bondings, PW8_1 has one CH··· O and one OH· · · p bonding interactions (see Figure 8). As can be seen in Figure 6, almost all the isomers have cage like structures, where most of them are folded cage. Previously, Janzen et al. [14] investigated the experimental spectroscopy of PhOH (H 2 O) 7,8 and reported some ab-initio calculations. Their located stable configurations are similar to PW8_4 and PW8_5, lying 0.4 and 0.5 kcal/mol in this work, respectively (see Figure 6). Similarly, Roth et al. [13]   similar to PW8_5 and PW8_14, for their spectroscopic study of the solvated phenol in water. Apart from these three isomers that have been reported previously, the remaining isomers of Figure 6 are located in this work for the first time. Examination of the configurations shows that the configurations of phenolwater octamer are similar to the structures of neutral water nonamer, where the OH of phenol is considered to be the ninth water molecule [33]. In addition, we noted that the hydrogen bond network governs the stability of the isomers. In the low-lying energy structures, the phenol OH group establishes two strong hydrogen bondings with the surrounding water molecules. However, for the less stable structures, the phenol OH group is only involved in one hydrogen bonding (see Figure 6). Moreover, low-lying energy structures have strong interaction with the phenyl group through CH··· O and OH· · · p interactions (see Figure 8).

located isomers
Regarding the structures of the solvated phenol decamer and dodecamer, PhOH(H 2 O) 10 and PhOH(H 2 O) 12 , we have reported in Figure 7 only the three most stable structures to avoid cumbersome results. For the decamer, 20 isomers have been located within a relative energy cutoff of 3.6 kcal/mol, while for the dodecamer, 25 different isomers have been located on its PES within a relative energy cutoff of 4.9 kcal/ mol. The complete list of the located isomers of the decamer and the dodecamer are presented in the supporting information. Similar to the case of lower size clusters, it is noted that the most stable configurations of the decamer and dodecamer are the isomers where the water molecules establish at least two strong hydrogen bonds and one CH· · · p bonding interaction with the phenyl group (see Figures 7 and 8). For both decamer and dodecamer, the low-lying energy level on their PESs is degenerated. The difference between the two isomers lies in the spatial orientation of the solvating water molecules. The results show that in the less stable structures, the solvating water molecules establish only one hydrogen bond with phenol. Overall, it is noted that the higher the interaction of the water molecules with phenol, the higher the stability of the corresponding isomer, no matter the internal configuration of the solvating water molecules.

QTAIM analysis of non-covalent bondings
To provide more insights about the interaction between the phenol and the solvating water molecules, we performed a QTAIM analysis on the most stable structures of the studied clusters. The located critical points and bond paths of the   Figure 8. For the specific case of the phenol-water monomer, we have additionally reported the electron density (ρ) two-dimensional contour and the atomic basins in the phenyl plane (see Figure 9). The bond paths, the critical points, the electron density 2D contour, and the the atomic basins of PW1 indicate the presence of two non-covalent intermolecular bondings, OH··· O and CH··· O hydrogen bondings. Quantitative data at all bond critical points (BCPs) are provided in the supporting information.

low-lying energy isomers are reported in
The positive values of the Laplacian of the electron density, ∇ 2 r (0.1329 au and 0.0260 au), at the bond critical points (of OH··· O and CH··· O hydrogen bondings, respectively) are indicative of non-covalent interactions. In addition, for the two BCPs, the electron density is evaluated to be 0.0364 au and 0.0071 au, respectively. These values indicate that the OH··· O hydrogen bonding is considerably stronger than the CH··· O hydrogen bonding in PW1. We have reported in Table 1 the interval of ρ and ∇ 2 r at bond critical points of the most stable configurations of the phenol-water clusters.
It can be seen in Figure 8 that the number of non-covalent interactions in phenol-water clusters is dominated by OH··· O hydrogen bondings. The OH··· O hydrogen bondings are also found to be the strongest non-covalent interactions of PhOH (H 2 O) n based on the values of the electron density at BCPs (see Table 1). Besides, it is also reflected in Table 1 that the CH··· O hydrogen bondings are less stronger than the OH··· O hydrogen bondings in all the studied clusters. In addition to the OH··· O and the CH··· O hydrogen bondings, two other non-covalent interactions can be identified in phenolwater clusters: the OH· · · p and the O··· C bonding interactions (see Figure 8 and Table 1). Among all the non-covalent bondings, the O··· C bonding interactions are the weakest interactions in phenol-water clusters. We have noted that only a few of the OH· · · p and the O··· C bonding interactions are identified in this work. Parthasarathi et al. [8] have reported the QTAIM study of phenol-water monomer, dimer and trimer based on HF/6-31G electron density. The authors found that only the phenol-water monomer and phenol-water trimer have each one CH··· O hydrogen bonding. For all the studied clusters, they have only reported OH··· O hydrogen bondings [8]. This could be attributed to the quality   of their electron density which is calculated at a relatively poor computational level of theory. In addition, their results could also be attributed to the missing isomers of the phenol-water clusters which are located in this work for the first time.

Temperature-dependent isomers' probabilities
Temperature effects on the stability of the investigated clusters are assessed using the Boltzmann distribution. The probabilities (relative populations) of different isomers of the phenolwater clusters are calculated for temperatures ranging from 20 to 400 K. As only one isomer of the monomer is located, there is no need to determine the probability. The calculated probabilities of the phenol-water dimer are plotted as function of temperature in Figure 10. The relative population or probability is the Boltzmann distribution given in Equation (5). The results show that the most stable phenol-water dimer, PW2_1, is the most favoured isomer for all the investigated temperatures. The isomers PW2_2 and PW2_3 contribute in trace to the population of the phenol-water dimer at high temperatures. The population of the dimer shows that the isomers PW2_4 and PW2_5 have a negligible contribution to the cluster's population. As noted in our previous works, the isomers with relative energies of higher than ∼1.0 kcal/mol have a negligible contribution to the cluster's population.
In addition to the phenol-water dimer, we have plotted in Figure 11 the temperature-dependence probabilities of the isomers of the studied clusters. For all the studied clusters, the results show that the most stable configuration is the most favoured at all temperatures. Nevertheless, at high temperatures, other isomers contribute significantly to the population of the clusters. For most of the clusters, there is a competition between the isomers at high temperatures. This result is found be to in agreement with previous works on furan clusters and dimethylsulphoxide clusters [53,54]. We have learned from the study of the temperature effects that the isomers that significantly contribute to the population of the clusters have their relative energies within ∼1.0 kcal/mol. This has been also remarked in our previous work on dimethylformamide clusters [55]. Although the structures of phenol-water clusters have similar configurations to those of neutral water clusters, their temperature-dependence follows different trends. It has been found that several isomers compete with the population of neutral water clusters, with some favoured at low temperatures and others more favoured at high temperatures [51,52,[56][57][58]. However, in the case of the phenol-water cluster, the most stable configurations dominate the population of the clusters. This indicates that the phenyl group has a significant influence on the temperature-dependence trend.

Hydration free energies and enthalpies of phenol
Before presenting the calculated hydration free energies and enthalpies of phenol, we would like to properly situate the reader about these properties/energies. The hydration free energies are important to understand the chemical reactivity and chemical kinetics of molecules in solvents. They are also involved in several chemical and biological processes as outlined in our previous work [42]. In this work, we use the cluster continuum 2 solvation model (CCM) to compute the absolute hydration free energy and enthalpy of phenol. As pointed out in the methodology, the CCM has been widely used to calculate the absolute hydration free energy with high accuracy as compared to experiments [34,38,39]. In the cluster continuum model, a few explicit solvent molecules around the solute are treated quantum mechanically, while the remaining is considered as a continuum medium. This model is quantum mechanically representing the solvation of a solute in a liquid.
Using the cluster continuum solvation model as outlined in the methodology section, we have calculated the absolute solvation free energy and the absolute solvation enthalpy of phenol in water. The solvation free energy and enthalpy are calculated using Equations (2) and (3), respectively. For each cluster size, we used all the possible configurations as presented in Section 3.1 weighted by their canonical probabilities. The structures of neutral water clusters are retrieved from our previous works and re-optimised at the ωB97XD/aug-cc-pVDZ level of theory. For the neutral water clusters, only the most stable configurations have been used. To be used in Equations (2) and (3) the gas phase geometry of phenol has been optimised at the ωB97XD/aug-cc-pVDZ level of theory. The solvation free energy and enthalpy of phenol in water are calculated for temperatures ranging from 20 to 400 K. We reported in Figure 12 the variation of the solvation free energy and enthalpy of phenol in water as function of the cluster size n at room temperature (298.15 K).
As can be seen in Figure 12, the calculated solvation free energies and enthalpies for different cluster sizes can be regarded as values oscillating around an average. The fact that the values oscillate around an average indicates that the explicit solvation of phenol has a negligible effect on the calculated hydration free energy and enthalpy. The estimated average values of the hydration free energy and enthalpy of phenol are −1.7 kcal/mol and −12.1 kcal/mol, respectively. Previously, an estimate of −6.6 kcal/mol has been reported by Gallicchio et al. [19] for the hydration free energy of phenol. The same value has been reported earlier by Cabani and coworkers [59]. As compared to theirs, our estimated hydration free energy is underestimated (less negative). Our most negative estimate of the hydration free energy (−3.4 kcal/mol found for n = 3) is also slightly underestimated as compared to the value estimated by Gallicchio et al. [19].
When no explicit water molecule is used, the hydration free energy of phenol is evaluated to be −4.3 kcal/mol. It is worth noting that the method used by Gallicchio et al. [19] is among the experimental techniques used to estimate the hydration free energy. Therefore, our calculated hydration free energy is slightly underestimated.
Regarding the hydration enthalpy of phenol, a previous experimental estimate of −13.6 kcal/mol has been reported by Cabani and coworkers [59]. Although our estimated hydration free energy is found to be slightly underestimated, our estimated hydration enthalpy is in good agreement with experiment. However, an improvement in the estimates would be interesting. To improve the agreement between our estimates and the experiment, one could use a high-level ab-initio method which could be expensive for such a large system. Another way is to benchmark several DFT functionals and select the one providing better agreement with the experiment. In addition, Guedes et al. [60] have estimated several values of the hydration enthalpy of phenol depending on the potential model used in their Monte Carlo simulations. Their estimations vary from −14.9 to −17.8 kcal/ mol [60].
To get further insights into the accuracy of the level of theory and the implicit solvation model, we tried a different level of theory and a different solvation model. Thus, we performed a single-point calculation on the most stable configurations for each value of n. The single-point calculations are performed using the SCS-MP2/aug-cc-pVTZ level of theory, and the SMD solvation model. The free energies and enthalpies at the SCS-MP2/aug-cc-pVTZ are obtained using the frequencies calculated at the ωB97XD/aug-cc-pVDZ level of theory. The calculated solvation free energy and enthalpy, for each cluster size n, are reported in the supporting information (Table S1). Similar to the case of ωB97XD/aug-cc-pVDZ, the results at the SCS-MP2/aug-cc-pVTZ level of theory show that the solvation free energy and enthalpy oscillate around an average. The averaged values of the solvation free energy and enthalpy are evaluated to be −0.8 and −11.3 kcal/ mol, respectively. It comes out from this investigation that the difference between the predictions at these two levels of theory is less than 0.9 kcal/mol.
After examining the effects of explicit solvation, we have investigated the effects of temperature on the hydration free energy and enthalpy for temperatures between 20 and 400 K. The estimated hydration free energy and enthalpy of phenol as function of temperature for n = 1−12 are reported in Figures 13 and 14, respectively. The oscillation of the estimated hydration free energy and enthalpy can be seen in Figures 13  and 14. The results show that the hydration enthalpy for a given cluster size is slowly varying with the change in temperature. This could indicate that the hydration enthalpy is temperature independent. Elsewhere, the hydration free energy is varying almost linearly as a function of temperature (see Figure 13). This indicates that the variation of the hydration free energy of phenol depends on the slope of the curves (which is the hydration entropy). Thus, one can state that the variation of the hydration free energy of phenol is entropically driven. It is worth noting that this behaviour of the hydration free energy and enthalpy of the proton has been reported in our previous works [40,61]. We have found that the solvation free energy of the proton in ammonia varies linearly as a function of temperature, while the solvation enthalpy of the proton in ammonia is found to be temperature independent [40]. Recently, we calculated the adsorption free energy of aniline onto coronene as a function of temperature [61]. It has been found that the adsorption free energy is linearly varying with temperature and entropically driven [61].

Conclusions
In this work, we provided the absolute hydration free energy and the absolute hydration enthalpy of phenol in water for temperatures ranging from 20 to 400 K. To undertake this investigation, we generated initial configurations of phenolwater clusters for n = 1 to n = 12. The generated configurations have been optimised at the ωB97XD/aug-cc-pVDZ level of theory. The results show that the configurations of the phenol-water clusters are similar to those of neutral water clusters without the phenol added. However, the configurations phenol-water clusters are folded as compared to neutral water clusters, due to the interaction of water molecules with the phenyl group. To understand the nature of the interaction between the water molecules and phenol, we performed a QTAIM analysis on the most stable structures of the phenolwater clusters. The analysis shows that the structures are stabilised by strong OH··· O hydrogen bondings and weak CH··· O hydrogen bondings. In addition, OH· · · p and O··· C bonding interactions are also identified.
The located structures of phenol-water clusters are used to compute the solvation free energy and enthalpy. Before that, we examined their relative population (probabilities) for temperatures ranging from 20 to 400 K. The results show that the most stable configurations dominate the population of the clusters. In addition, we have found that few structures compete with the population of clusters at high temperatures. After calculating the hydration free energy and enthalpy of phenol, the results show that the explicit solvation has a negligible effect on the estimated values. The estimated average values of the hydration free energy and enthalpy of phenol are −1.7 and −12.1 kcal/mol, respectively. The temperaturedependence of the solvation free energy and enthalpy shows that the hydration enthalpy of phenol is temperature independent, while the hydration free energy varies linearly as a function of temperature.