*3.2. Calculation of Interaction Energy*

To determine the energy of the studied intermolecular interactions with fluorine, all complexes obtained from the PDB were divided into subgroups based on the following: (i) angle of HB (ranged from 90◦ to 180◦ with a step of 10◦), (ii) distance of HB (ranged from 2.5 to 4 Å with a step 0.1 Å), (iii) donor type of HB (OH, NH, +NH, CH), and (iv) whether fluorine is bonded to an aromatic (Far) or aliphatic (Fal) carbon. Then, one representative PDB complex was randomly selected from each rectangle defined by unit distance and

angle change. In the next step, all the selected systems were visually inspected to identify those in which the HB with fluorine was not the main stabilizing interaction and the number of adjacent supporting interactions was the smallest. These identified complexes were used to calculate the interaction energy of the fluorine-containing HB for the given geometric parameters. Figure S1 illustrates the distribution of complexes in a given HB distance–angle interval.

The appropriate ionization states at pH 7.4 ± 0.5 were determined using Epik v3.4 [28,29]. Using an in-house script, the structure of the ligand and the amino acid participating in the interaction was extracted, the missing atoms (included in the peptide bond) were added, and their positions were optimized (OPLS3 force field). The interaction energy of the identified complexes with fluorine was calculated using three commonly used approaches as follows.

The first method, named difference approach or Diff, works based on the assumption that the total interaction energy equals the energy required to separate two interacting molecules. Thus, the energy between the HB donor (X) and the acceptor (Y) is calculated as the difference between the total energy of the X···Y complex and the sum of the total energies of its frozen components [18,30]:

$$E\_{\text{int}} = E(X \cdot Y) - \left[ E(X) + E(Y) \right]^2$$

The energies of the separated molecules, as well as that of the complex, were calculated in Gaussian G16 software [31], using the Minnesota functional M06-2X [32,33] and Karlsruhe basis set def2-TZVP [34]. The polarizable continuum model (PCM) (solvent = water) was used for the calculation [35,36].

The second approach works based on Bader's quantum theory of "atoms in molecules" (QTAIM). In this approach, the topological analysis of electron density was carried out in AIMAll program [37]. The electron density calculated in Gaussian G16 at the M06-2X/def2- TZVP level and the PCM (solvent = water) were used in the analysis. The energy of the noncovalent bonds detected in the crystal structures was calculated using the Espinosa equation as follows:

$$E\_{int} = \frac{1}{2}v(r)$$

where *Eint* is the energy of the interatomic interaction and *v*(*r*) is the kinetic energy at the bond critical point (BCP). The above equation can be used for all types of HBs, van der Waals interactions, and weak interactions such as H···H and C–H···O [38].

The third approach works based on the energy decomposition analysis. Bonding was analyzed using the extended transition state (ETS) method [39], with the natural orbitals for the chemical valence (NOCV) scheme [40–42]. In this approach, the total energy of bonding between the interacting molecules (ΔEint) is divided into different components as follows:

$$
\Delta E\_{\rm int} = \Delta E\_{\rm dist} + \Delta E\_{\rm el} + \Delta E\_{\rm Pauli} + \Delta E\_{\rm orb} \tag{1}
$$

where Δ*Edist* is the energy required to promote the separated fragments from their equilibrium geometry to the structure they will take up in the complex, Δ*Ee*<sup>l</sup> is the energy of the electrostatic interaction occurring between the two fragments in the supermolecule geometry, Δ*EPaul*<sup>i</sup> is the energy of repulsion between the occupied orbitals of the two fragments, and Δ*Eorb* or the orbital interaction term refers to the energy of the stabilizing component due to the final orbital relaxation. All calculations were performed using the Amsterdam Density Functional (ADF) program [43–46], using the ETS-NOCV scheme. The Becke, Lee, Yang, and Parr exchange-correlation functional with the Grimme dispersion correction (B3LYP-D3) was used. A standard double-ζ STO basis containing one set of polarization functions was adopted for all the electrons (TZP). The total (electronic) bonding enthalpies (ΔE = Δ*Eint*) did not include the zero-point energy (ZPE) additions, finite temperature contributions or basis set superposition error corrections (BSSE).
