The Advantages of Flexibility: The Role of Entropy in Crystal Structures Containing C–H···F Interactions

Molecular crystal structures are often interpreted in terms of strong, structure directing, intermolecular interactions, especially those with distinct geometric signatures such as H-bonds or π-stacking interactions. Other interactions can be overlooked, perhaps because they are weak or lack a characteristic geometry. We show that although the cumulative effect of weak interactions is significant, their deformability also leads to occupation of low energy vibrational energy levels, which provides an additional stabilizing entropic contribution. The entropies of five fluorobenzene derivatives have been calculated by periodic DFT calculations to assess the entropic influence of C–H···F interactions in stabilizing their crystal structures. Calculations reproduce inelastic neutron scattering data and experimental entropies from heat capacity measurements. C–H···F contacts are shown to have force constants which are around half of those of more familiar interactions such as hydrogen bonds, halogen bonds, and C–H···π interactions. This feature, in combination with the relatively high mass of F, means that the lowest energy vibrations in crystalline fluorobenzenes are dominated by C–H···F contributions. C–H···F contacts occur much more frequently than would be expected from their enthalpic contributions alone, but at 150 K, the stabilizing contribution of entropy provides, at −10 to −15 kJ mol–1, a similar level of stabilization to the N–H···N hydrogen bond in ammonia and O–H···O hydrogen bond in water.


Cell optimisation
To provide a suitable level of geometrical convergence for the phonon calculations, experimental crystal structures were geometry-optimized as described in the main text.The optimized unit cell dimensions are compared to experimental values in Table S1.The structures are all monoclinic except for 1F, which is tetragonal.

Comparison of aromatic H…F interactions
Searches of the Cambridge Structural Database 1 were conducted to investigate whether the formation of C-H⋯F interactions is simply a consequence of the optimization of aromatic interactions.Relevant interactions were defined as those shorter than the sum of the contributing van der Waals radii.Data without 3D coordinates, R-factors above 5%, containing disorder, errors, powder data or ions were filtered out of searches and hydrogen bond lengths were normalised to neutron values.Searches were conducted on the CSD 2023.2.0 dataset.The average distances shown in Table S2 no significant systematic difference between the lengths of C-H⋯F contacts in aromatic and non-aromatic systems.

Convergence of phonon calculations
Phonon frequencies for all structures were calculated in the harmonic approximation using the linear response method implemented in CASTEP (as opposed to the finite displacement method based on supercells) using a grid of q-points with a spacing of 0.04 Å −1 .The results were then interpolated onto a finer grid of 16 x 16 x 16 q-points.[Note: q-points define the positions in reciprocal space where phonon frequencies are calculated.]Our approach was based on the following convergence tests.
Calculations of increased computational expense were applied to monofluorobenzene, 1F.An increase in the density of q-points from 0.04 Å -1 to 0.03 Å -1 , yielding an increase from 3 to 12 q-points, was applied.Levels of interpolation from 1 x 1 x 1 to 16 x 16 x 16 were then applied to the results obtained at 0.04 and 0.03 Å −1 ; for even-numbered interpolations (e.g. 2 x 2 x 2), a grid shift was applied to include the Γ-point in the interpolation grid.The value of TS at 150 K was calculated for each interpolated q-point set.The results are tabulated in Table S3.
The value of TS converged with a grid of 16 × 16 × 16 to within 0.01 kJ mol −1 for calculations based on the 0.04 and 0.03 Å −1 grids.At this level of interpolation, the difference between the values of TS for the two sets of calculations was of the order of 0.1 kJ mol −1 .
A comparison of the DoS and calculated INS spectra for the 0.04 and 0.03 Å −1 q-point sets at the highest level of interpolation are provided in Figure S1 and Figure S2 respectively; note that the integrals of both DoS plots are normalised to three times the number of atoms in one unit cell.The differences are very small.
The minor impacts on calculated entropies, DoS plots and INS spectra led us to apply a q-point spacing of 0.04 Å −1 interpolated to 16 x 16 x 16 for all further calculations.
Table S3: Effect of q-point density on the value of TS at 150 K for 1F.

SAPT Calculation convergence
Dimer energies were calculated for mapping potentials using symmetry adapted perturbation theory (SAPT) through the program Psi4. 2 Convergence testing of the most appropriate level of theory for these calculations was completed using a dimer of 1,2-difluorobenzene, 1,2F, contact C, a bridged C-H⋯F interaction (below).Combinations of the double, triple and quadruple zeta basis sets with the augmentation levels of August, July and June from the 'calendar' set (reducing backwards in month for each diffuse shell removed from the calculation) were tested at all levels of SAPT truncation. 3The results are shown in Table S4.
The wall times refer to the SAPT2+3 calculations, which were given 150 GB of memory and 30 processors.Convergence of SAPT energies against time taken for the calculations is shown in Figure S3.
SAPT2+3 was considered the most reliable for all computations.From Table S4 it is clear that, especially the for the lower basis sets, the augmentation level of the calculation has a comparatively modest impact on wall times whilst giving significant improvements to calculated energies.Increases in basis set quality however are very computationally expensive.For the august augmentation level, the shift from a double zeta to quadruple zeta basis set would change the calculation time of a 21-point potential well from ~1 day to ~85 days.From Figure S3 an initial significant improvement to the energy quickly flattens to an area where modest energy improvements require multiple hour increases in wall time.It is from the beginning of this flattening that the aug-cc-pVDZ basis was chosen.

Inelastic Neutron Scattering plots
Inelastic neutron scattering plots were made using AbINS through the Mantid workbench GUI. 4 Temperatures used for simulations were taken from the experimental data log and were: 1F-24 K, 1,2F-26 K, 1,4F-20 K, 1,3,5F-30 K and 1,2,4,5F-20 K. AbINS applies a semi-empirical powder averaging model to the relative atomic displacements and frequencies from the phonon calculations.In this model every phonon is treated as an independent quantum harmonic oscillator. 5Atoms are weighted by their neutron scattering cross-sections.Fundamental quantum modes up to 10 were considered, each increasing order having diminishing contribution.The specific instrument resolution function for the TOSCA instrument is built into AbINS and this was used to convolute the theoretical spectra.AbINS also removes experimentally observed overtones and accounts for the lower scattering atom types in the experimental data.Scaling of intensities was completed visually for peaks closest 1000 cm -1 .Data for averages of the front and back detector were used.Comparisons of experimental and simulated INS for all structures over smaller and larger spectral ranges are given below.

Density of states plots
Phonon density of states plots for all structures over ranges of 0-1500 and 0-3250 cm −1 are presented below.These plots were produced through Gaussian smearing in Materials Studio with a defined smear width c of 0.1 THz (~3.3 cm −1 ).Intensities of all plots are normalised such that the total integral is equal to 3N, where N is the number of atoms.

Partial density of states plots
Partial density of states plots for selected atom(s) were derived from phonons at the Γ-point.
The output of the phonon calculation contains N, where N is the number of atoms in the model, eigenvector descriptions for x, y and z perturbations in atom positions for each of the 3N vibrational modes.For each eigenvector the fraction of the total motion attributed to the selected atom(s) is calculated and used to provide a total fractional contribution towards each vibrational mode.The partial density of states plots were produced via the same Gaussian smearing regime described above.The height of the resulting peaks represents the fraction of motion attributed to the selected atom type.Partial density of states plots for all structures at the Γ-point are presented below in the same ranges as the total density of states plots above.

Example: Calculation of the contributions to internal energy and entropy by occupation of energy levels at 150 K in 1,2F
The vibrational frequencies ν at the Γ-point were obtained by periodic DFT as described in Section 2.3 of the main text.Calculation of the contribution to the entropy (S) and thermal energy (U) from each vibration at a temperature T = 150 K is accomplished using the following equations, where k, h and R are Boltzmann's constant, Plank's constant and the gas constant and the value of ν in Hz is obtained from the mode wavenumber (ω in cm −1 ) as equal to 100cω, where c is the speed of light: The values of U and TS were calculated for each of the 141 non-acoustic modes as shown in Table S6.The magnitude of the total entropy calculated at the Γ-point (Table 2 in the main text) is obtained by summing the values in the TS column.

Justification for the neglect of U and zero-point energy
The equations for U and TS given in Section 8 can be plotted as a function of ω for T = 150 K to allow the influence of these quantities on the free energy to be compared (Fig. S24).The plot shows that TS is much more sensitive to the values of the frequencies of low-energy vibrational modes than is U. Small differences in external mode vibrational frequencies make little difference to the contribution of mode to U, and alternative packing arrangements will have similar values of U.This reflects the observation made by Nyman and Day 6 that the difference in heat capacity between polymorphs are very small.In the harmonic approximation the zero-point energy is given by ½Σhν and is dominated by the highest energy internal vibrational modes.These modes can vary between polymorphs, but usually not by much and differences are a minor contributor to the free energy differences between polymorphs, 6 though they can assume significance in the cases where strong intermolecular interactions such as H-bonds are present as these can perturb internal bonds more substantially than the weak interactions that are considered in the present work.

Figure S24 .
Figure S24.Comparison of the contributions of low-energy vibrational modes to entropy and internal energy with T = 150 K.

Table S1 :
Comparison between the experimental and optimised unit cells for all structures studied.

Table S2 :
A comparison of the average bond length of intermolecular H⋯F interactions in the CSD and the aromatic and non-aromatic subsets of this.

Table S4 :
Wall times and calculated energies for different levels of SAPT theory.
FigureS3: Convergence of relative SAPT energies.The y-axis has been extended towards negative % differences for clarity of values which would otherwise lie very close to the x-axis.

Table S5 :
Symmetry unique interactions within the first coordination spheres of fluorobenzenes.All energies are given in kJ mol −1 .U = lattice energy.

Table S6 :
Calculation of vibrational contributions to internal energy (U) and entropy (S) at temperature T = 150 K.The zero-frequency acoustic modes (1-3) have been omitted from the calculation.