Hydration dynamics and IR spectroscopy of 4-fluorophenol

Halogenated groups are relevant in pharmaceutical applications and potentially useful spectroscopic probes for infrared spectroscopy. In this work, the structural dynamics and infrared spectroscopy of para-fluorophenol (F-PhOH) and phenol (PhOH) is investigated in the gas phase and in water using a combination of experiment and molecular dynamics (MD) simulations. The gas phase and solvent dynamics around F-PhOH and PhOH is characterized from atomistic simulations using empirical energy functions with point charges or multipoles for the electrostatics, Machine Learning (ML) based parametrizations and with full ab initio (QM) and mixed Quantum Mechanical/Molecular Mechanics (QM/MM) simulations with a particular focus on the CF- and OH-stretch region. The CF-stretch band is heavily mixed with other modes whereas the OH-stretch in solution displays a characteristic high-frequency peak around 3600 cm−1 most likely associated with the –OH group of PhOH and F-PhOH together with a characteristic progression below 3000 cm−1 due to coupling with water modes which is also reproduced by several of the simulations. Solvent and radial distribution functions indicate that the CF-site is largely hydrophobic except for simulations using point charges which renders them unsuited for correctly describing hydration and dynamics around fluorinated sites. The hydrophobic character of the CF-group is particularly relevant for applications in pharmaceutical chemistry with a focus on local hydration and interaction with the surrounding protein.


Introduction
Fluorination -and halogenation in general -are common chemical modifications for pharmaceuticals. Approximately 20 % of all small molecule drugs used in medicinal chemistry contain X = F, Cl, Br, or I or a combination thereof. Among these compounds halogenated phenyl rings constitute an important class. 1 Because of the directionality of the interaction along the C-X bond due to the sigma hole, halogenation has emerged as one of the essential chemical modifications in medicinal materials, 2-4 and supramolecular chemistry. 5,6 By changing the halogen atom, the interactions with the environment can be tuned and the hydrophobicity around the modification site can be modulated. 3,[7][8][9][10][11][12] The importance of halogenation as a fundamental concept in medicinal chemistry is highlighted by the improved binding affinities of several ligands towards their receptors. 13,14 Recently, halogenation has also been employed in the context of protein modifications, such as for insulin, to fine-tune thermodynamic stability and affinity to the insulin receptor. 15 A halogen bond "[..]occurs when there is evidence of a net attractive interaction between an electrophilic region associated with a halogen atom in a molecular entity and a nucleophilic region in another, or the same, molecular entity." 16 Hence, halogen atoms act as electrophiles and can form an attractive interaction with a nucleophilic counterpart. Based on the analysis of the molecular surface electrostatic potential (ESP), 17 the "halogen bond" was also associated with a "σ-hole bond" 18 which is a noncovalent interaction between a covalentlybonded halogen atom X and a negative site, e.g. a lone pair of a Lewis base or an anion. 17 Such a "bond" involves a region of positive electrostatic potential, i.e. the σ-hole, and as an extension of one of the covalent bonds to the atom. The σ-hole arises as a consequence of the anisotropy of the ESP around the halogen atom. The strengths of the interactions generally correlate well with the magnitudes of the positive and negative electrostatic potentials of the σ-hole and the negative site. As fluorine has the largest electronegativity and the lowest polarizability, for some time it was in fact assumed that there is no σ-hole and that therefore fluorine is not involved in halogen bonding at all. 17,19,20 However it is now well established that it can have a positive σ−hole and form halogen bonds when it is linked to strongly electron-withdrawing groups including Cl, Br, and I. 21,22 Moreover, there is also experimental evidence for fluorine engaging in halogen bonding. 23 Introducing a fluorine atom into organic molecules can cause major changes in the physicochemical properties such as solubility, chemical reactivity and biological activity compared to non-fluorinated analogues. 24 In particular, fluorine often replaces hydrogen in organic molecules but the size and stereoelectronic influences of the two atoms (hydrogen vs. fluorine) are quite different albeit it is often regarded as isosteric substitution. 4 In bio-inorganic and medicinal chemistry, the formation of intermolecular O-H/F-C and N-H/F-C hydrogen bridges was assumed to be important in binding fluorinated compounds to enzyme active sites. 25 Such interactions affect enzyme ligand binding affinity, selectively coupled with the changes in pharmaco-kinetic properties by fluorine substitution. 10,24 The effects of fluorine substitution on the related pharmaco-kinetic properties like lipophilicity, volatility, solubility, hydrogen bonding and steric effects affect the resulting compound binding, absorption, transport and hence the related biological activity. 26 A variety of functional groups, including C-H, C-OH, C=O, and C≡N, have utilized the C-F bond as a bioisostere. 27 However, it is difficult to generalize the relative ability of fluorine to act similar to a hydrogen or hydroxy group, and different factors must be considered in each case. The van der Waals radius of fluorine (1.47Å) lies between that of oxygen (1.57 A) and hydrogen (1.2Å) and as it is the element with the highest electronegativity, the C-F bond is almost identical to C-OH in terms of bond length and polarity. Despite its three electron pairs, the C-F bond interacts more weakly with the environment compared to an oxygen atom and is better described as "weakly polar" rather than "hydrogen bonding". 27,28 In pharmacological applications the replacement H→F is often considered to avoid metabolic transformation due to the high stability of the CF bond. Examples are drugs interacting with P450 for which fluorination has been widely used to block metabolic transformations. 27 Given the different qualitative characterizations outlined so far, a more molecularly refined picture of the energetics and dynamics involving fluorinated model compounds will be valuable. The present work considers hydrated fluoro-phenol (F-PhOH) as a typical representative. Using linear infrared spectroscopy together with computational characterizations at different levels of theory the structural dynamics and spectroscopy of F-PhOH is charac-terized. The computations use advanced empirical force fields including multipolar interactions, a machine-learned, neural network-based representation of a full-dimensional PES, mixed quantum mechanics/molecular mechanics and ab initio MD simulation techniques.
The measured data from infrared spectroscopy can be directly compared with the computational results. Also, the solvent dynamics is investigated based on frequency fluctuation correlation and spatial distribution functions. In this article, first, the methods are described followed by the results for the spectroscopy of the -CF and -OH from experiments and simulations. Then, the solvent structure is analyzed from radial distribution functions and from 2-dimensional solvent distributions and finally, conclusions are drawn.

Classical Molecular Dynamics Simulations
All classical MD simulations were performed with CHARMM. 29 The bonded parameters are based on CGenFF 30 except for the CF and OH bond for which a Morse potential was used to describe their anharmonicity. To that end, a scan along the CF bond was performed at the MP2/aug-cc-pVTZ level starting from an optimized structure of F-PhOH at this level of theory. The energy of 49 points was computed on a grid ranging from r = 0.75Å to r = 5.55Å in increments of 0.1Å. Then, the energies were fitted to a Morse potential 2 which yields parameters D 0 = 136.316 kcal/mol, r 0 = 1.349 A, and β = 1.603Å −1 . For the OH bond the calculated Morse parameters are D 0 = 120.234 kcal/mol, r 0 = 0.971Å, and β = 2.088Å −1 . To realistically describe the electrostatic interactions, a multipolar (MTP) 31-34 model was also used with MTPs on all heavy atoms up to quadrupoles and point charges for all hydrogen atoms. These parameters were fitted to the electrostatic potential using a fitting environment, 35 see Tables S1 to S3.
Simulations of F-PhOH and PhOH were carried out in a cubic box of 30 3Å3 (28 3Å3 for simulations with the NN-PES, see below) using TIP3P 36 water molecules. Minimization, heating, and equilibration procedures for 40 ps were employed to prepare the system. Production simulations of 5 ns were run in the N V E ensemble at 300 K using a Velocity Verlet. 37 The time step was ∆t = 1 fs and every 5 snapshot was recorded. Lennard-Jones interactions were computed with a 12Å cutoff switched at 10Å. 38 The electrostatic interactions for the monopoles (point charges) are treated using Particle-Mesh Ewald 39 (PME) with grid size spacing of 1Å, characteristic reciprocal length κ = 0.32Å −1 , and interpolation order 6. All bonds involving hydrogen atoms are constrained via the SHAKE algorithm. 40,41 Additional MD simulations were also carried out for PhOH in water with the same setup that was used for F-PhOH in order to directly compare their spectroscopy and solvent structure.
For the simulations with the NN-based PES (see below) the atomic simulation environment (ASE) was used. 42 The van-der-Waals interactions were those from the CGenFF 30 parametrization and the fluctuating charges are from the PhysNet representation, see below.
For both terms interactions the cutoff distance is at 14Å and switched between 13 to 14Å.
To avoid artifacts of the electrostatic Coulomb force in the cutoff range, the Coulomb force at the distance of 14Å was shifted to zero in accordance to the shifted forces method. 43 In the gas phase 1000 trajectories, each 200 ps in length, are run to obtain an ensemble average.
The N V E simulations are run at 300 K initialized from random momenta corresponding to a Maxwell-Boltzmann distribution, with a time step of 0.5 fs, equilibrated for 50 ps and propagated for 200 ps. Simulations in solution N V T using Langevin 44 thermostat at 300 K are performed for 20 trajectories of 100 ps each with a time step of 0.2 fs to obtain a total of 2 ns for PhOH and F-PhOH in solution, respectively. The IR spectra are then calculated from the dipole-dipole moment autocorrelation function [45][46][47] and averaged over all 1000 trajectories.

Instantaneous Normal Mode Analysis
From the production simulation, 10 6 snapshots were taken as a time-ordered series for computing the frequency fluctuation correlation function (FFCF). The FFCF was determined from instantaneous harmonic vibrational frequencies based on a normal mode analysis. Such instantaneous normal modes (INM) are obtained by minimizing F-PhOH while keeping the surrounding solvent frozen. Next, normal modes were calculated using CHARMM for 5 modes (ν 1 to ν 5 in ascending order) between 1100 to 1400 cm −1 in terms of participation ratio of CF stretch in that particular mode. In a separate analysis step, the participation ratios of the CF, CO, and CH stretch and the COH bending coordinates to these 5 normal modes were determined.

Frequency Fluctuation Correlation Function and Lineshape
From the INMs the frequency trajectory ω i (t) and the FFCF, δω(0)δω(t) was computed.
From the FFCF the line shape function is determined within the cumulant approximation. To compute g(t), the FFCF is numerically integrated using the trapezoidal rule and the 1D-IR spectrum is calculated from 48 where ω is the average transition frequency obtained from the distribution, T 1 = 1.2 ps 49 is the vibrational relaxation time and α = 0.5 is a phenomenological factor to account for lifetime broadening. 48 From the FFCF, the decay time is determined by fitting the FFCF to a general expression 50 The two systems were comprised of the solute (F-PhOH or PhOH), and 331 and 311 water molecules, respectively. The system size was selected so that a direct comparison with the full QM simulation of F-PhOH in water can be made, and assess the impact of quantum description of the solvent including polarization effects on the geometric and spectral prop- The systems were first equilibrated classically, using AMBER18. 66 F-PhOH and PhOH were modelled with the GAFF2 force field, 67,68 while the TIP3P model was used for water. Following an initial minimization, the two systems were equilibrated in the isothermal-isobaric (N P T ) (300 K, 1 atm) ensemble with the Berendsen barostat 69 and Langevin dynamics 44 for pressure and temperature control respectively, followed by simulations in the NVT ensemble with Langevin dynamics, for a total of 100 ns. A time step of 2 fs was employed. In view of the small periodic box size, a reduced real space cutoff of 7Å was used for the nonbonded interactions.
In both systems, the solute was treated at the QM level and the solvent at the classical (MM) level. The QM setup, and the simulation time step were the same as described above for the full QM simulations. The QM/MM MD simulation protocol was also similar to the one described for the full QM simulation, apart from the use of two separate Nosé-Hoover thermostats for the QM and MM parts respectively, during the equilibration with BO and CP MD. For the F-PhOH system, the equilibration and production runs lasted 10.1 ps and 35.9 ps respectively, while for the PhOH system they lasted 12.6 ps and 25.0 ps, respectively.
During the production phase, frames were saved with the same frequency as in the full QM simulations (0.48 fs).

Machine-Learned Potential Energy Surface (PES)
To validate in particular the PC-and MTP-based simulations using an empirical force field a complementary model based on a machine-learned PES was also pursued. For this Phys-Net, 70 a deep neural network (NN) of the message passing type, 71 was used to obtain a representation of the potential energy for both PhOH and F-PhOH. PhysNet uses Cartesian coordinates and nuclear charges to learn an atomic descriptor for the prediction of energies, forces, partial charges and molecular dipole moments to describe chemical systems and their properties, such as infrared spectra.
PhysNet was trained on ab initio energies, forces and dipole moments calculated at the MP2/6-31G(d,p) level of theory using Molpro 72 according to the protocol reported in Ref. 70 The reference data, containing different geometries for both molecules, is generated from MD simulations at 50, 300 and 1000 K using CHARMM force field (5000 geometries each yielding a total of 30000 geometries) and extended with geometries obtained from normal mode sampling 73 at temperatures between 10 and 2000 K (6600 geometries for each molecule). The complete data set thus contains 43200 PhOH and F-PhOH structures. The performance of the PhysNet PES is reported in Figure S1 which shows the correlation between the reference MP2 and the PhysNet energies for a test set of 3700 randomly selected points with

Experimental Infrared Spectroscopy
Experimental spectra of PhOH and F-PhOH in water have been measured in attenuated total reflection (ATR) geometry using a home-constructed ZnSe prism in a Bruker Tensor 27 FT-IR spectrometer. Identical concentrations of 0.76 M have been used for both samples, which is the saturation limit for F-PhOH. As a background spectrum, pure water has been measured as well and subtracted from the PhOH and F-PhOH spectra. The subtraction procedure is deemed reliable in the shown spectra, with the exception of the range between ≈ 3100 cm −1 to 3500 cm −1 , where the very strong OH stretch vibration dominates. In fact, in this range a dispersion-shaped response is observed, which is attributed to the effect the samples have on the water spectrum, i.e., a shift to somewhat higher frequencies due to overall weaker water-water hydrogen bonds.

Gas Phase Spectra
First, the performance of the PC-and MTP-based empirical force fields, of PhysNet, and of the DFT-BLYP/DCACP based ab initio MD simulations was assessed for F-PhOH in the gas phase. For this, MD simulations of F-PhOH in the gas phase were carried out and the power and infrared spectra were determined and compared with experiments. and computed (in the gas phase) spectra for F-PhOH in the range of (1100-1400 cm −1 ) and (2600 to 4000 cm −1 ). Panel A: experimental (black) spectrum in CCl 4 from Ref. 61 extracted using g3data. 76 The arrows in the right hand column refer to the experimental frequencies for the CH and OH stretches. 61 Panels B and C: IR from the Fourier transform of the total dipole moment correlation function and CF/OH power spectra from PC (blue) and MTP (red) simulations for F-PhOH in the gas phase. Panel D: global (solid) and CF/OH (dotted) power spectrum in violet from QM simulations in the gas phase. Panel E: IR spectrum (maroon) obtained from ML/MM MD simulations.
The measured CF and CO stretching modes occur mainly at 1226 and 1262 cm −1 , respectively, while they couple to one another and potentially to other modes. According to the analysis 61,75 (see Table 1), the CF stretch is coupled to the in plane bending of the ring and also the C-H bend while the CO stretch couples to the CC and CF stretching vibrations. Therefore, the CF stretch is intimately coupled with other modes and for that reason it is not possible to assign a local CF stretching mode to one particular frequency. Table 1: Vibrational frequencies in cm −1 for PhOH (gas phase) and F-PhOH (in CCl 4 ) in the range of 1100-1400 cm −1 . 61,74,75 The contributions (in terms of local deformations) to each vibrational mode indicate strong mixing and are those from the literature. 61,74 The assignment of the bands has been made on the basis of the calculated potential energy distribution. 61,77 Symbols ν and δ refer to stretching and bending modes, respectively.
In summary, the gas phase spectrum from finite-T MD simulations find comparable patterns for the frequencies in the 1100 to 1400 cm −1 region when compared with experiment. It is also found that the CF stretch is coupled to other modes in this spectral range and no local mode for this motion can be assigned.

Spectroscopy and Dynamics of F-PhOH in Water
After assessing the energy functions considered in the present work, the spectroscopy of F-PhOH in solution is analyzed, see Figure 2. The experimentally measured spectrum from the present work is the black trace in Figure  For ML/MM MD simulations a broad band is observed between 1220 to 1290 cm −1 in the IR spectrum which is blue shifted compared to the experiment. The CF-power spectrum with peak maximum at 1273 cm −1 (dotted maroon trace in Figure 1E) indicates that part of this broad IR-lineshape is due to CF-stretching motion. Moreover, the CO power spectrum has a first peak at 1250 cm −1 which contributes to the broad IR peak below 1300 cm −1 and a second, prominent peak at 1323 cm −1 , see Figure 1E. The overlapping peaks of the CF-and CO-power spectra indicate coupling between the two types of motion.
The high-frequency region of the experimental spectrum for solvated F-PhOH (Figure 2A) above and below 3000 cm −1 involves a broad absorption extending from ∼ 2700 to 3100 cm −1 , a region between 3100 and 3500 cm −1 that is not reliable due to the dominating background from the OH stretch vibration of bulk water which can not be subtracted off completely (grey area in Figure 2A), and a high-frequency feature at 3643 cm −1 assigned to the free OH vibration presumably originating from F-PhOH but possibly also from H 2 O.
This latter assignment is less likely, though, as the same sharp peak appears for PhOH in solution but shifted to the red by ∼ 40 cm −1 . If the signal was due to water it is expected to occur at closer frequencies given the similarity of the solutes. and PhOH in the gas phase shows that this signal corresponds to the "free OH stretch" vibration, see Figure 1. Features at 3150 cm −1 in the infrared spectrum ( Figure 1B) are due to the CH-stretch vibrations which could be brought into better agreement with experiment by slight reparametrization of the force constants. These features are not present in the OH-power spectra, Figure 1C, as expected which confirms the assignment to the CH-stretch vibration. The structured spectrum below 3000 cm −1 is present in both, the PC and MTP simulations, albeit with lower intensity.
For the full QM and QM/MM simulations the global power spectra ( Figure 2D solid violet and orange) find a signal centered at 3000 cm −1 which is typical for the CH-stretch modes. At higher frequency (∼ 3250 cm −1 ) the OH-stretch vibration is located which is confirmed by the OH-power spectra (dotted violet and orange traces). However, no signal in the 3600 cm −1 region is present which suggests that the "free OH" signature in these simulations, expected around 3500 cm −1 from the QM gas phase simulations ( Figure 1D), is absent.
Simulations with the PhysNet energy function primarily find the high frequency -OH stretch at 3803 cm −1 with broad, largely unstructured undulations below 3000 cm −1 . It is likely that the MP2/6-31G(d,p) level is not sufficient for quantitatively describing the spectroscopy of F-PhOH. Accounting for a frequency scaling of 0.94 for harmonic frequencies 78 shifts all frequencies to the red which is more consistent with the experimentally determined spectra.
Specifically, the 3800 cm  Table S4 for the frequency maxima. For obtaining the instantaneous normal modes the water environment was frozen and the structure of the solute was optimized, followed by a normal mode calculation. The band positions compare well with the frequencies from

Comparison of the Spectroscopy for Hydrated PhOH and F-PhOH
The experimental spectra for F-PhOH (black) and PhOH (green) in water are reported in Figure 2A. For the experimental spectra between 1100 cm −1 and 1400 cm −1 pronounced differences between the two compounds are found. Most prominently, the single band with maximum at 1242 cm −1 for PhOH, which is one of the "X-sensitive modes" involving considerable motion of the phenol ring and the CO group, 80 is shifted to the red for F-PhOH and split into at least two (at 1201 and 1222 cm −1 ), but possibly several more peaks, some of which overlap with the peak from PhOH. Other features, such as the broader band with peak maximum at 1381 cm −1 for PhOH are also shifted to the red (band maximum at 1368 cm −1 ) for F-PhOH, see Figure 2A.
Consistent with experiment, the number of spectral features for F-PhOH is larger than for PhOH in the 1200 to 1300 cm −1 range. However, none of the computed spectra display the pronounced double-peak structure above 1200 cm −1 for F-PhOH with the broad peak for PhOH to the blue of it. The simulations using MTPs find a single absorption between the low-and high-frequency absorption in F-PhOH, see Figure 2B   which points towards a largely hydrophobic behaviour of the CF site. This is also consistent with notions from pharmaceutical chemistry in that a CF group reduces ligand solvation and increases its hydrophobicity. 27 The number of water molecules within distance r is reported in Figure S4 and shows that for small F-water separations (r ≤ 4Å) the occupation from simulations with MTP and PhysNet is similar but clearly below that of QM/MM and QM simulations (which are identical) whereas for larger separations (r ∼ 5Å) that from QM simulations approaches the MTP and PhysNet simulations. For g H OH −O W from the QM and QM/MM simulations (blue and red traces in Figure 5C) the position of the first maximum is almost identical whereas the height of the first peak differs.
Both pair correlation functions have a first minimum around 2.5Å with an amplitude close to 0 which suggests that during the 25 ps simulation one solvent water molecule is strongly bound to the OH-group of F-PhOH and does not exchange with the surrounding solvent.
As a consequence the "free phenolic -OH stretch" (spectroscopic feature around and above 3600 cm −1 , Figure 2 Figure S7 which manifestly shows symmetric water distribution according to the underlying dynamically averaged spatial symmetry of the solute due to rotation of the OH group around the -COH axis, see Figure S8. Two-dimensional solvent distribution functions were generated from the positions of the water-oxygen atoms around F-PhOH. For that, the structures of the 5000 snapshots were oriented with C1 in the origin, the C1-C4 bond along the x-axis the [C1,C4,H] atoms in the xy−plane. A 2-dimensional histogram of the water positions was generated and then refined from kernel density estimation using Rstudio. 86 The distribution of the solvent water around the CF-part of F-PhOH is comparatively flat for all simulations with PC, MTP and PhysNet when contrasted with the COH-moiety of the solute, see Figure 6. It is also found that for PhOH (top line) the solvent distributions resemble those for F-PhOH (bottom line).
The asymmetry in the solvent distribution around the -COH group is due to the reference atoms chosen for the superposition of all structures. If the phenol-oxygen atom is excluded in reorienting the structures the 2d distribution becomes manifestly symmetric ( Figure S7) and also clarifies that the -OH group of the solute rotates on the time scale of the simulations, see Figure S8. Specifically, rotation of the -COH group is an activated process and occurs on the ∼ 100 ps time scale from MTP and ML/MM simulations.
The weak interaction between solvent and solute around the CF-site of F-PhOH can also be gleaned from the behaviour of the frequency fluctuation correlation functions (FFCFs). For this, the FFCFs were determined from the frequency trajectories ω i (t) for the five modes in the CF-stretch region, i.e. ν 1 to ν 5 , between 1100 cm −1 and 1400 cm −1 and for the OH-stretch vibration, see Figures S9 and S10. The FFCFs contain information about the coupling between a particular mode and the environmental dynamics. The fitting parameters for ν 1 to ν 5 in Table S5 show that the fast correlation is generally τ 1 ∼ 0.1 ps whereas the longer time scale ranges from τ 2 = 0.28 ps to τ 2 = 0.83 ps. Such short correlation times point towards weak solvent/solute interactions close to the CF-site. The static components ∆ 0 are very small, too, which also imply rapid, unspecific dynamics around the CF-group.
For the OH-stretch vibration, the fast correlation time is τ 1 ∼ 0.1 ps while the longer one is τ 1 ∼ 0.5 ps which are comparatively short. Compared with the FFCFs for vibrations involving the CF-stretch the t = 0 amplitude and the static offset ∆ 0 is larger by one to two orders of magnitude which points towards somewhat slower dynamics around the OH bond of F-PhOH.

Discussion and Conclusion
The present work discusses the infrared spectroscopy of and solvent distribution around F-  80 In CCl 4 solution and the pure liquid an additional band is at 3350 cm −1 . 80 The present work finds a sharp peak at 3596 cm −1 which is assigned to the "free" phenol-OH stretch in a non-hydrogen bonded environment. This is supported by the observation that the peak is only marginally shifted from the PhOH OH-stretch in vapor and the fact that the spectroscopic feature is sharp and therefore can not be due to water.
More recently, infrared spectra were recorded for PhOH complexed with variable numbers of water molecules in the gas phase, 87,88 in matrices, 89 and for PhOH at the air/water interface using vibrational sum frequency generation (SFG). 81 The cluster studies all report the phenolic-OH stretch vibration at frequencies above 3000 cm −1 whereas the experiment at the air/water interface assigns a very broad signature in the SFG signal extending from 2550 cm −1 to 3500 cm −1 to the OH-stretch mode. 81 This finding is consistent with the present experiments which report a broad absorption for both, F-PhOH and PhOH extending down to ∼ 2700 cm −1 which is also assigned to the phenol-OH stretch for water-coordinated -COH, see Figures 2 and S2.
Hydration around the -CF group as characterized by 1d-and 2d-solvent distribution functions does not feature any pronounced hydrogen bonding unless a conventional PC model is used which, however, exaggerates the directed interaction along the CF bond. On the other hand, solvent water molecules form intermittent H-bonds with the -OH group of F-PhOH which leads to the broad spectroscopic response below 3000 cm −1 observed experimentally and reproduced by several of the computational models.
From the perspective of halogenic modifications of organic frameworks used for drug design the present work suggests that the local hydration of -CH and -CF groups is comparable (see Figures 4 and S5), supporting the notion that "F behaves like a large H-atom". 27 Although the C-F bond has three isolated electron pairs, it has weaker electrostatic interactions compared to O due to the small size and high electronegativity, which compromises its hydrogen bonding ability and can thus better be described as a weakly polar interaction rather than a hydrogen bond. Such insights are important guidelines for rational drug discovery as they provide a basis for directed modification and evolution of ligands with specific interactions in protein binding sites.

Supporting Information Available
The supporting information provides tables for the force field parametrizations (Tables S1 to S3), frequency maxima for frequency distributions (Table S4), and parameters for the FFCFs (Table S5) together with Figures S1 to S10 the report quality of the PhysNet model, additional solvent distribution functions, dihedral time series and the FFCFs for the CF-and OH-stretch frequencies.     Average frequency ω of the asymmetric stretch in cm −1 , the amplitudes a 1 and a 2 in ps −2 , the decay times τ 1 and τ 2 in ps, and the offset ∆ 0 in ps −2 . The FFCF for the OH-stretch is reported in Figure S10.  Figure S1: Correlation between the ab initio and PhysNet energies for a set of 3700 randomly selected points averaged over 980/982 trajectories with R 2 = 0.9999 and an root mean square error of 0.0037 eV.  Figure S5: The radial distribution function between the water-oxygen atoms and the two carbon atoms flanking the COH group in PhOH from a 5 ns simulation with MTP. It is demonstrated that the solvent distribution is converged. Comparison with the red trace in Figure 4A indicates that hydration around the CF and CH groups is similar with the first maximum at a similar C-O W separation.   Figure S7: 2d-distribution from MTP simulations with TIP3P water. The structures are reoriented with reference to the solute heavy atoms, but excluding the phenolic-oxygen atom. This differs from the selection used in Figure 6 and allows to visualize the symmetric distribution of the solvent water molecules around the -COH group. . With the shorter cutoff there is no water molecule hydrogen bonded to the rotatable COH group for most of the simulation time and occasionally one water is sufficiently close to form an H-bond. For the longer cutoff water presence is frequent and H-bounding (green crosses) is more prevalent. The alteration between H-bonded and non-H-bonded -COH group explains the spectroscopy seen in Figure  S2. A direct correlation between rotation of the CCOH angle and water occupation can not be established convincingly although it is likely that rotation of the COH group requires water to detach from it. The dihedral time series is locally averaged over a time window of 50 ps. Simulations with ML/MM MD yield a comparable time series with transitions on the ∼ 100 ps time scale.  Table S5. Logarithmic scale is chosen for the y-axis.  Table 2. Logarithmic scale is chosen for the y-axis.