Solving the enigma of weak fluorine contacts in the solid state: a periodic DFT study of fluorinated organic crystals.

The nature and strength of weak interactions with organic fluorine in the solid state are revealed by periodic density functional theory (periodic DFT) calculations coupled with experimental data on the structure and sublimation thermodynamics of crystalline organofluorine compounds. To minimize other intermolecular interactions, several sets of crystals of perfluorinated and partially fluorinated organic molecules are considered. This allows us to establish the theoretical levels providing an adequate description of the metric and electron-density parameters of the C-F⋯F-C interactions and the sublimation enthalpy of crystalline perfluorinated compounds. A detailed comparison of the C-F⋯F-C and C-H⋯F-C interactions is performed using the relaxed molecular geometry in the studied crystals. The change in the crystalline packing of aromatic compounds during their partial fluorination points to the structure-directing role of C-H⋯F-C interactions due to the dominant electrostatic contribution to these contacts. C-H⋯F-C and C-H⋯O interactions are found to be identical in nature and comparable in energy. The factors that determine the contribution of these interactions to the crystal packing are revealed. The reliability of the results is confirmed by considering the superposition of the electrostatic potential and electron density gradient fields in the area of the investigated intermolecular interactions. This journal is © The Royal Society of Chemistry.

Introduction

Non-covalent interactions involving halogens have gained increasing interest of researchers in the last decade because of their significance in chemistry,[1,2] biology,[3] drug design[4,5] and materials science.[6-8] There is ample literature on the X-ray and theoretical studies of organic crystals with intermolecular halogen⋯halogen contacts (“halogen bonds”).[9-12] The F⋯F interactions are often considered together with other C–Hal⋯Hal–C interactions, where Hal = Cl, Br and I.[13,14] However, due to its low polarizability[15] and tightly contracted lone pairs fluorine tends to fall out of the dependences typical of halogen bonds.[16] This leads to a number of specific properties of the C–F⋯F–C interactions in comparison with C–Hal⋯Hal–C. (i) The anisotropic electronic distribution is negligible on fluorine in comparison with other halogens;[17] as a result, the directionality of the C–F⋯F–C interactions is insignificant compared to the C–Hal⋯Hal–C contacts.[18] (ii) The energy of the C–F⋯F–C interactions is less than that of the C–Hal⋯Hal–C interactions in molecular crystals. This follows from the experimental values of the sublimation enthalpy ΔHsub, which equals ∼16, ∼39 and ∼52 kJ mol−1 for crystalline CF4, CCl4 and CBr4, respectively.[19] (iii) Halogen bonds with Cl, Br and I are widely used as building blocks in crystal engineering,[20] while the C–F⋯F–C interactions do not play any structure-directing role. (iv) Fluorine tends to form C–H⋯F–C interactions rather than C–F⋯F–C contacts, whereas heavier halogens seem to prefer the formation of halogen⋯halogen interactions.[21]

Intermolecular C–H⋯F–C interactions usually occur in the crystal structures of organofluorine compounds.[22] In contrast to C–F⋯F–C, these interactions play an important role in directing supramolecular assembly in chemical and biological systems.[21,23] As far as we know, there has been no detailed comparison of the intermolecular C–F⋯F–C and C–H⋯F–C interactions in crystals of organofluorine compounds. Crystallographic analyses have led to controversial conclusions on the considered problem[24] due to the use of area and shape corrections for C–F⋯F–C contacts[18] and the cone correction for C–H⋯F–C contacts.[25] Comparison of C–F⋯F–C and C–H⋯F–C contacts obtained from the X-ray experiment is not straightforward, since normalization of the C–H bonds leads to a shortening of the contacts H⋯F by ∼0.15 Å.[26] In addition, the presence of a relatively short halogen⋯halogen contact in the crystal does not mean that there is an interaction between the atoms in question.[27] All this makes it necessary to refine the hydrogen positions in molecular crystals through geometry optimization by periodic DFT[28] or additional neutron diffraction analysis.[29] However, the number of papers devoted to the neutron diffraction analysis or periodic DFT computations of crystals of organofluorine compounds is very limited.[30,31]

The intermolecular C–F⋯F and C–H⋯F interactions were theoretically studied in detail in the gas phase.[32-34] Complexes of fluorinated methane, HF and F2 with simple molecules (H2O, C2H4etc.) were considered. All the computations point to dominating dispersive contributions, in particular for very weak C–H⋯F interactions with H⋯F distances larger than 2.5 Å.[33] The approaches based on semi-empirical methods, such as the PIXEL scheme,[35-37] gave conflicting results on C–H⋯F–C interactions in crystals. In some cases, these interactions were mainly of a dispersive nature,[36] in other cases they were characterized by a substantial coulombic contribution.[35] There are several reasons explaining the current situation: (i) limited applicability of non-periodic models for describing non-covalent interactions in crystals;[38] (ii) using experimental geometry in PIXEL calculations of crystals of organofluorine compounds; (iii) lack of systematic theoretical studies of C–F⋯F–C and C–H⋯F–C contacts in organic crystals. Indeed, the X23,[39] C60 (ref. 40) and Z20 (ref. 41) benchmark sets for noncovalent interactions in solids do not contain these contacts.

C–F⋯F–C and C–H⋯F–C contacts in molecular crystals are characterized by energies ∼2 and ∼5 kJ mol−1, respectively,[35,36,42-44]i.e. they are weak intermolecular interactions. Firmer conclusions on such interactions are expected from extended studies of complexes in crystals under the conditions in which other intermolecular interactions are minimized.[33] In particular the dependence of the properties of C–F⋯F–C C–H⋯F–C interactions on the underlying donor C(sp, sp2, sp3)–H and acceptor (primary, secondary or tertiary F–C moieties) parameters[36] can be elucidated from thermodynamic experiments.[33]

In the present work, an in-depth study of the intermolecular C–F⋯F–C and C–H⋯F–C interactions in crystals of organofluorine compounds has been carried out by periodic DFT calculations coupled with the experimental data on the structure and sublimation thermodynamics of the studied crystals. To minimize the other intermolecular interactions, several sets of organic crystals were considered: (1) perfluorinated compounds; (2) molecular crystals with dominant C–F⋯F–C/C–H⋯F–C contacts; (3) crystals with C–F⋯F–C, C–H⋯F–C and C–H⋯O interactions. To avoid the controversial conclusions of the previous theoretical studies and crystallographic analyses of the C–F⋯F–C/C–H⋯F–C interactions[33,35,36,45-47] the following steps were taken. First, the optimal level of periodic DFT computations was established and the relaxed molecular structures obtained by this level of approximation were used in subsequent calculations. Second, different theoretical approaches for evaluation of the ΔHsub value of crystals of perfluorinated molecules were considered. Third, a Bader analysis of theoretical crystalline electron density function was performed. Finally, superpositions of the gradient fields of electrostatic potential and of electron density and deformation density maps in the area of the studied intermolecular interactions were plotted to clarify the difference in the nature of C–F⋯F–C/C–H⋯F–C contacts. To the best of our knowledge, a periodic DFT study of the interplay between intermolecular C–F⋯F–C and C–H⋯F–C interactions in organofluorine crystals has not been done yet. The methodological aspects of different DFT-based approaches (including non-periodic calculations) for studying C–F⋯F–C and C–H⋯F–C interactions in crystals of organofluorine compounds are also discussed.

The competition between C–H⋯O and C–H⋯F–C has been the subject of heated discussions.[21,47-52] To clarify this problem, we combined the Cambridge Structural Database (CSD) screening[53] with the results of periodic DFT computations of the crystals of partly fluorinated molecules with C–H⋯O interactions.

The aim of this work is threefold. (1) To find a functional/basis set combination for periodic DFT calculations, which gives best results for the intermolecular C–F⋯F–C contacts in crystals and to establish a theoretical approach, which gives the reasonable ΔHsub values for crystals of perfluorinated molecules. (2) To determine the specific features of intermolecular C–H⋯F–C/C–F⋯F–C interactions in crystals of organofluorine compounds. (3) To clarify the role of the C–H⋯F–C contacts in crystal packing with/without C–H⋯O contacts. We start by introducing crystal structures containing organic fluorine as well as the used computational methods in Section 2. Section 3.1 describes the search for the level of DFT calculations which gives the best results for the metric and electron density parameters of the C–F⋯F–C contacts. Estimation of the lattice energy of crystals of perfluorinated molecules by different theoretical approaches is discussed in Section 3.2. The role of the C–H⋯F–C contacts in crystal packing is described in Section 3.3. The competition of the C–H⋯F–C interactions with nonconventional C–H⋯O bonds is discussed in Section 3.4. The applicability of Bader's theory to describing the intermolecular interactions and differences between C–H⋯F–C interactions and conventional hydrogen bonds is discussed in Section 4. Finally, we draw conclusions in Section 5.

Methods

Crystals under study

To study the nature of fluorine-based intermolecular interactions and interplay of interactions of particular types, several sets of molecular crystals are considered with the intention to minimize the possible contribution of other intermolecular contacts (Scheme 1). (i) Perfluorinated compounds with only C–F⋯F–C contacts: tetrafluoromethane (CF4),[54] hexafluorobenzene (C6F6)[55] and pentafluoropyridine (C5NF5).[56] (Fluorine and nitrogen atoms are considered to be isolobal with respect to the supramolecular architecture in the crystalline phase[57]). (ii) Molecular crystals with prevalent C–F⋯F–C/C–H⋯F–C contacts: 1,1,2,4,4-pentafluorobuta-1,3-diene (C4HF5),[58] 1,2,3,5-tetrafluorobenzene (m-C6H2F4),[23] 1,2,3,4-tetrafluorobenzene (o-C6H2F4)[59] and 2,3,5,6-tetrafluoropyridine (C5NHF4).[56] (iii) Crystals where a competition occurs between intermolecular C–F⋯F–C, C–H⋯F–C and C–H⋯O contacts: 2,3-difluoro-1,4-naphthoquinone (DFNAPQ),[60] 3-(3,5-difluorophenyl)-1-phenylprop-2-en-1-one (YICBES)[61] and 3,4-difluoro-6-phenyl-2H-pyran-2-one (KEGWEZ)[62]. The designations of the crystals given in bold; they will be used below.

Scheme 1

The molecular structures of the single-component crystals under study, divided into three classes by the type of predominant non-covalent interactions. Refcodes of the crystals are given in parentheses.

At the level of theory validation stage we additionally choose a smaller training set of crystals, where the lattice energies calculated by different approaches are compared against experimental ΔHsub values. For this purpose we select from ref. 19 the fluorine-containing crystals for which both X-ray diffraction data and ΔHsub are available, namely, CF4, C6F6, and pentafluorobenzoic acid (C6F5COOH).[63]

In order to reveal the applicability of non-periodic models to the crystals of partially fluorinated molecules, the DFT computations of non-periodic systems were carried out. For that purpose the dimers which form C–H⋯F–C synthons of types A and D[23] are considered. DFT computations of isolated dimers were performed using Gaussian 09 package.[64]

Periodic DFT calculations

Two different approaches are widely used for the atomic simulation of molecular crystals within periodic DFT calculations: the first approach exploits a linear combination of plane-wave functions as a basis set[28] and the second one adopts a linear combination of local atomic orbitals usually represented as Gaussian-type orbitals.[65] The pros and cons of each approach are discussed elsewhere.[66]

Two codes with plane wave basis sets were tested for computations of molecular crystals with perfluorinated compounds. The pseudopotential method is implemented in Quantum Expresso,[67] while ELK is based on the full-potential linearized augmented-plane-wave method.[68] Due to the difficulties associated with geometry optimization of the crystals with C–F⋯F–C contacts in Quantum Expresso and bad self-consistent field (SCF) convergence in some cases in ELK, the plane wave basis sets were not used below (for more information see ESI†).

In the CRYSTAL17 (ref. 69) calculations, the B3LYP and PBE functionals were employed with all-electron Gaussian-type localized orbital basis sets 6-31G**, pob-TZVP, and 6-31(F+)G**. In the latter diffuse functions with exponent factor equal to 0.1076 were added to the fluorine atoms. It is greater than the critical value 0.06 (ref. 70) below which problems in the SCF procedure arise.

London dispersion interactions were taken into account by using the semi-empirical D3 scheme.[71] Tolerances on energies which control the self-consistent field convergence for geometry optimizations and frequency computations were set to 1 × 10−8 and 1 × 10−10 hartree, respectively. The shrinking factor of the reciprocal space net was set to 4. The space groups and the unit cell parameters of the considered crystals obtained in the experimental studies[54-56,58-63] were fixed for the calculations. This is a common approximation for calculating molecular crystals with Hal⋯Hal contacts[27,72,73] as the change in the volume of a cell of molecular crystals appears to be insignificant after full optimization.[40,74]

Periodic DFT computations of molecular crystals sometimes lead to the appearance of imaginary frequencies.[75,76] This problem is usually solved by reducing the space symmetry to a minimum (P1 space group).[74] In the present study, such procedure was applied to the CF4 crystal (C2/c space symmetry group[54]).

Estimation of the lattice energy

The approaches to lattice energy Elatt evaluation can be divided into two groups. The methods of the first group consider the as difference between the total electronic energy of a relaxed bulk crystal (Ebulk) and an isolated molecule (Ei) in its relaxed geometry:Here Z indicates the number of molecules in the unit cell.

In order to obtain physically relevant results, the basis set superposition error needs to be taken into account in the Ei calculation either by Boys–Bernardi counterpoise correction[77] (MOLEBSSE keyword in CRYSTAL17 (ref. 70)) or geometric counterpoise (gCP) scheme recently proposed by Grimme et al.[39] It should be noted that the basis set superposition error evaluation is not straightforward for conformationally flexible molecules and crystals with more than one molecule in the asymmetric unit.[78] According to our computations, the geometry is virtually unchanged, and conformation energies of studied compounds are rather small (∼1–3 kJ mol−1). Therefore, Elatt values evaluated using eqn (1) can be compared with the experimental sublimation enthalpy of the considered crystals.

The second group of methods considers the lattice energy as sum of energies of all pair intermolecular interactions within the asymmetric unit:where Eint, is the energy of a particular intermolecular interaction. Indices j and i denote the atoms belonging to different molecules. In the present study, Eint, was evaluated using the empirical correlation basing on Bader analysis of crystalline electron density (QTAIMC). The calculation details for this scheme are given in Section 2.4. Other methods from this group used in present work include the Coulomb–London–Pauli scheme developed by Gavezzotti in AA-CLP[79] and PIXEL[80] versions, as well as CE-B3LYP approach.[81] In all three methods, the energy of each pair interaction is calculated as a sum of four terms (coulombic, polarization, dispersion and repulsion). While AA-CLP uses predefined force field parameters for different atom types, PIXEL and CE-B3LYP use molecular wavefunctions and orbitals taken from the gas-phase calculations.

To ensure the consistence of obtained results with already existing data on fluoroorganic compounds,[35-37] default calculation parameters were used for AA-CLP, PIXEL and CE-B3LYP. In AA-CLP and PIXEL, C–H/F–H/OH distances (if any) were normalized to standard values prior lattice energy estimation using the built-in function.

Evaluation of energy of intermolecular interactions in crystals using QTAIMC

In QTAIMC,[82] the particular intermolecular interaction is associated with the existence of a bond critical point (BCP) between the pair of atoms. The energy of each specific interaction is considered totally independent of the others. The effects of the crystal environment, long-range electrostatic forces, etc. are taken into account implicitly, via the crystalline electronic wave function, and are coded in the BCP features.

The energy of the particular intermolecular interaction Eint was evaluated according:[83]Gb is the local electronic kinetic energy density at the BCP.[82] It is obtained from the crystalline wave function using TOPOND14.[85]Eqn (3) with k = 0.429 yields reasonable Eint values for conventional H-bonds, C–H⋯O and C–H⋯F interactions in crystals.[84,86,87]

Superposition of gradient fields of electrostatic potential and electron density; deformation density maps

To reveal the nature of the H⋯F interactions in organic crystals, we have plotted the deformation density maps and superposition of the gradient fields of electrostatic potential and electron density for crystalline C4HF5 and KEGWEZ. C4HF5 has a short H⋯F contact, ∼2.35 Å, while KEGWEZ has H⋯F and H⋯O contacts with similar distances, ∼2.5 Å. For this purpose, we computed the theoretical structure factor lists for relaxed crystals using the CRYSTAL17. The theoretical structure factors are calculated for each possible hkl index characterized by sin θ/λ < 1.3 Å. The Hansen & Coppens[88] formalism is used for multipole refinement implemented in the MOLLY program.[89] After that, the multipole parameters are used to calculate gradient fields of the electrostatic potential and the electron density and deformation density maps by the WinXPRO program.[90] The superposition of gradient fields of the electrostatic potential and the electron density is usually studied to reveal the electrostatic contribution to different interactions in crystals.[91] The theoretical background is given in ESI.†

Results

Two terms are widely used in the literature to describe weak non-covalent forces in crystals, namely “contacts” and “interactions”. In contrast to “contact”, the term “interaction” implies the existence of a BCP,[27] which is localized by QTAIMC.

To reveal applicability of non-periodic models to the description of C–H⋯F–C interactions in organic crystals the dimers extracted from crystalline m-C6H2F4 and the cocrystal of C6F6 and anthracene[23] were computed at the PBE/6-31(F+)G** level (Fig. 1). The results obtained by this approach differ from the experimental data. In the m-C6H2F4 dimer (synthon A), the calculated values of the H⋯F/C⋯F distances are shorter than the standard experimental values by ∼0.15 Å (Table S1†). In the C6F6 – anthracene dimer (synthon D) the calculated values of the H⋯F/C⋯F distances agree with the experimental data. However, the spatial orientation of molecules is radically different from the orientation of these molecules in a two-component crystal C6F6 – anthracene (Table S1† and Fig. 1). The reasons for these disagreements are the disregard of the effects of the crystal environment in non-periodic models. We conclude that gas-phase computations have limited applicability for a description of intermolecular C–H⋯F–C interactions in organic crystals. The difference between experimental gas-phase (gas electron diffraction) and crystal-phase structures of perfluorinated species was also emphasis in ref. 92 and 93.

Fig. 1

C–H⋯F–C synthons referred to in this study. C–H⋯F dimer motif in the o-C6H2F4 crystal (synthon A) and C–H⋯F–C dimer motif in the cocrystal of C6F6 and anthracene (synthon D). Spatial orientation of molecules in the C6F6 – anthracene dimer, obtained as a result of optimization at the PBE/6-31(F+)G** approximation (lower panel). H⋯F interactions are denoted by dotted lines.

Selecting the optimal level of periodic computations

To select the optimal level of periodic DFT computations (functional/basis set), calculations were made for five crystals with multiple C–F⋯F–C interactions: CF4, C6F6, C6F5COOH, C5NHF4, and C5NF5. The results of calculations of intermolecular C–F⋯F–C distances in different approximations are given in Table S2.† The C–F⋯F–C contacts with the F⋯F distances significantly less than the sum of their van der Waals radii (∼2.94 Å[94,95]) are better reproduced in calculations than the longer C–F⋯F–C contacts. The inclusion of dispersion correction D3 has little effect on the theoretical geometric characteristics of the crystals under consideration. The best results are achieved using the 6-31(F+)G** basis set, where (F+) denotes diffuse functions added to fluorine atoms (Table S2†). Replacing the basis set with pob-TZVP[96] does not significantly affect the theoretical C–F⋯F–C distances but greatly increases the computational cost. B3LYP/6-31(F+)G** and PBE-D3/6-31(F+)G** give the best results for the F⋯F distances in the considered crystals. The obtained results are in agreement with ref. 97, according to which the B3LYP/6-311+G* level shows excellent agreement with the experimental structures of isomeric fluoro thiourea. The disadvantage of using the B3LYP functional with the 6-31(F+)G** basis set is the poor convergence of the SCF procedure.

The theoretical values of the F⋯F distances, the electron density and its Laplacian at the BCP of the C–F⋯F–C contacts evaluated using PBE-D3/6-31(F+)G** are compared with the experimental data in Fig. 2 and Table S3.† This level provides an adequate description of the electron density distribution for the F⋯F distances less than 2.94 Å.

Fig. 2

The deviations of experimental F⋯F distances ΔR(F⋯F) in crystalline CF4, C6F6, C5NF5, C5NHF4 (black circles) and C6F5COOH (red triangles) from theoretical F⋯F distances R(F⋯F) computed at the PBE-D3/6-31+(F)G** level. The horizontal dotted line corresponds to zero deviation. The vertical dotted line is a sum of van der Waals radii for fluorine atoms (∼2.94 Å (ref. 94 and 95)).

The electron density ρb at the BCP is less than 0.003 a.u. for the contacts with the F⋯F distances larger than 2.94 Å (Table S3†). This BCP value is too small to be determined with certainty by the existing theoretical methods and precise X-ray diffraction experiments.[86,98] Identification of the C–F⋯F–C bonding at the F⋯F distances larger than 2.94 Å involves the use of advanced experimental methods.[30] It should be noted that in the crystals of perfluorinated molecules, a significant part of the C–F⋯F–C contacts is characterized by the F⋯F distances greater than 2.94 Å, see Fig. 2 and ref. 34. The use of the term “interaction” at such distances must be very careful.

The Elatt values of crystals of perfluorinated molecules

The theoretical Elatt values of the crystals of perfluorinated molecules (CF4, C6F6 and C6F5COOH) were evaluated using different theoretical approaches and compared with the experimental ΔHsub values in Table 1. In the case of eqn (1), PBE-D3/6-31(F+)G** is the best level of approximation, as shown in Table S4.† Theoretical Elatt values of crystalline C6F6 and C6F5COOH significantly exceed the experimental values for most DFT methods (Table S4†). This result does not contradict the literature data,[99] according to which root-mean-square errors can reach 20 kJ mol−1 if no attempt is made to correct ΔHsub for finite temperature effects. Using the gCP correction[39] in Elatt calculations does not improve the agreement with the experiment (Table S4†) as this approach was developed for organic crystals containing H, C, N, and O atoms. Neglecting the D3 correction results in negative lattice energies of “van der Waals” crystals CF4 and C6F6 (Table S3†). This can be explained by underestimating weak non-covalent interactions in DFT calculations, e.g. see Table S1 in ref. 39 and 100. The results are consistent with ref. 101. According to it, “⋯ halogen bonds are just too complex for simple dispersion corrections”. Thus, the correct description of the energy properties of perfluorinated crystals involves the use of a more accurate many-body dispersion energy correction.[102]

Comparison of the Elatt values calculated using different approaches with the experimental ΔHsub values of crystals of perfluorinated molecules. The units are kJ mol−1

Crystal	Eqn (2)			Eqn (1)	ΔHsub[19]
	Eqn (3)a	AA-CLP	CE-B3LYP	PBE-D3/6-31(F+)G**
CF4	68.4 (20.5)b	19.6	25.4	17.1	14.0–17.0
C6F6	82.1 (33.4)	44.3	46.8	69.8	46.0–49.8
C6F5COOH	147.1 (92.6)	89.5	111.0	143.6	92.0

E int was computed using eqn (3) for intermolecular interactions with ρb > 0.003 a.u.; the electron-density characteristics at the BCP were evaluated using periodic electron density obtained at the PBE-D3/6-31(F+)G** level.

The Eint values in parenthesis were computed with coefficient value equals to 0.129 in eqn (3) for F⋯F and F⋯O contacts.

The Elatt values evaluated using eqn (2) overestimate ΔHsub regardless of the level of periodic DFT calculations (Tables 1 and S5†). Applicability of this approach to the Elatt assessment of crystalline CF4 is not straightforward, due to a large number of the F⋯F contacts with R(F⋯F) > 2.94 Å (Fig. 2) and ρb < 0.003 a.u. However, this approach works better for crystalline C6F5COOH because the main contribution to the lattice energy is made by OH⋯O contacts between carboxyl groups. Moreover, in C6F5COOH, most of the F⋯F contacts are shorter than 2.94 Å (Fig. 2). The inclusion of vibrational effects would alter the obtained results because the C–F bond length would be slightly elongated. However, the investigation of the temperature effect on the Elatt value is beyond the scope of the present study.

We have tried to find the coefficient value in eqn (3) for the F⋯F contacts, which gives “reasonable” Elatt values (the standard deviation from ΔHsub is less than ∼10 kJ mol−1 (ref. 99)). The coefficient value equals ∼0.129 for perfluorinated molecules. This value can be justified as follows. The ΔHsub value of a CF4 crystal is ∼2.5 times smaller than that of CCl4, and the Cl⋯Cl distances are systematically larger than F⋯F. According to[27], the energy of Cl⋯Cl interactions in crystals is well described by eqn (3) with k = 0.429.

The semi-empirical schemes for Elatt evaluation PIXEL,[80] AA-CLP[79] and CE-B3LYP[81] are widely used due to their robustness and user-friendliness. We have tested their performance on the set of the studied compounds. Eqn (2) with the AA-CLP force field[79,103] gives the best results for the Elatt values of perfluorinated crystals (Tables 1 and S4†) while CE-B3LYP overestimates it for CF4 and C6F5COOH. The PIXEL method, considered to be a more advanced version of AA-CLP, also fails to provide correct lattice energies for CF4 and C6F5COOH (Table S5†). The possible reason for poor performance of PIXEL and CE-B3LYP methods may be caused by the 6-31G** basis set used to estimate the molecular orbital overlap in CE-B3LYP[81] or electron density distribution in the molecule in PIXEL.[80] The basis sets chosen for parameterization of both methods lack diffusion orbital functions, which are essential for the agreement with the experimental distances and sublimation enthalpies (Section 3.1).

C–F⋯F–C vs. C–H⋯F–C in crystals containing only C, H and F atoms

Sections 3.3 and 3.4 will exclusively use molecular parameters in crystals optimized in the PBE/6-31(F+)G** approximation. The metric and electron-density features, as well as the energy of the C–H⋯F–C and C–F⋯F–C interactions in crystals containing C, H and F atoms only, are given in Tables 2 and S3.† Theoretical H⋯F distances range from ∼2.35 to ∼2.70 Å in the studied crystals and are systematically shorter than F⋯F distances (Fig. 2), in accordance with the X-ray data of crystals containing only C, H and F atoms.[21,23,43,104,105] Indeed, the shortest intermolecular C–H⋯F–C and C–F⋯F–C distances equal to ∼2.20 Å[106] and ∼2.60 Å,[107] respectively. The energies of the C–H⋯F–C contacts (5–7 kJ mol−1) are found to be larger than the energies of the C–F⋯F–C contacts (<4 kJ mol−1), in agreement with the literature data.[35,36,42-44]

Metric and topological characteristics of C–H⋯F–C and C–F⋯F–C contacts in crystals containing only C, H and F atoms, evaluated by periodic DFT calculations at the PBE-D3/6-31(F+)G(d,p) level. ρb is the electron density, ∇2ρ(r) is the Laplacian of the electron density; λ1, λ2, λ3 are the electron density curvatures at the BCP of the H⋯F/F⋯F contacts. Eint is the energy of the corresponding contact, evaluated by eqn (3)

Contacta	H⋯F/F⋯F distanceb, Å	ρ(r), a.u.	∇2ρ(r), a.u.	λ 1, a.u.	λ 2, a.u.	|λ1/λ3|,	E int c, kJ mol−1
C4HF5
C9–H1⋯F9	2.335 (2.420)	0.011	0.030	−0.014	−0.013	0.246	7.3
C5–F9⋯F1	2.900 (2.909)	0.005	0.029	−0.005	−0.005	0.132	1.8
C5–F9⋯F13	2.940 (2.909)	0.006	0.030	−0.005	−0.003	0.132	1.3
C1–F5⋯F5	2.880 (2.889)	0.007	0.036	−0.007	−0.002	0.156	1.6
m-C6H2F4
C4–H1⋯F1	2.550 (2.530)	0.0055	0.025	−0.006	−0.006	0.167	5.7
C6–H2⋯F4	2.578 (2.614)	0.0053	0.025	−0.005	−0.005	0.143	5.6
C4–H1⋯F3	2.590 (2.684)	0.0064	0.029	−0.006	−0.005	0.150	6.4
C2–F2⋯F4	2.897 (2.923)	0.0057	0.030	−0.006	−0.005	0.146	2.2
o-C6H2F4
C5–H1⋯F3	2.467 (2.557)	0.0068	0.031	−0.007	−0.007	0.159	7.3
C5–H1⋯F4	2.520 (2.950)	0.0060	0.027	−0.006	−0.006	0.154	6.3
C6–H2⋯F2	2.529 (2.566)	0.0062	0.027	−0.007	−0.006	0.175	6.3
C3–F3⋯F4	2.836 (2.858)	0.0064	0.033	−0.006	−0.006	0.133	2.4
C1–F1⋯F2	2.867 (2.950)	0.0062	0.032	−0.006	−0.006	0.140	2.3
C4–F4⋯F4	2.882 (2.943)	0.0071	0.034	−0.007	−0.007	0.146	2.5

Only the contacts shorter than the sum of van der Waals radii are shown, see Fig. 3 for atomic numeration.

The experimental value is given in parentheses (C–H bonds were normalized to 1.09 Å).

In the case of the F⋯F contacts the coefficient value in eqn (3) equals to 0.129 (Subsection 3.2).

A comparison of the pattern of bond paths in the m-C6H2F4 crystal (middle panel of Fig. 3) with the literature data (Fig. 2 in ref. 23) shows the limited applicability of the geometric criterion to the localization of weak intermolecular C–F⋯F–C interactions in crystals of organofluorine compounds. F⋯F distances greater than 2.94 Å (ref. 34) characterize the vast majority of these crystals. Locating BCPs at such distances using QTAIMC is almost impossible.

Fig. 3

Fragments of crystalline C4HF5 (upper panel), m-C6H2F4 (middle panel) and o-C6H2F4 (lower panel). The C–H⋯F–C and C–F⋯F–C interactions are denoted by dotted lines. Two types of C–F⋯F–C contacts occur in C4HF5 crystal, namely, intra- and interlayer ones. The latter are shown with red dotted lines.

The C–H⋯F–C interactions are characterized by a higher electron density value than the C–F⋯F–C interactions (Table 2) in accordance with the literature.[35,108-110] The C–H⋯F–C and C–F⋯F–C contacts can be referred to the pure closed-shell interactions,[32,111] which are characterized by electron density shifted to the interacting atoms, due to the low values of the λ1/λ3 at the BCPs (Table 2).[112] However, consideration of the BCP characteristics does not allow identifying any qualitative differences between the interactions under consideration. To establish the nature of the C–H⋯F–C interactions and to distinguish difference between them and other weak interactions, like C–F⋯F–C, the superposition of the electrostatic potential and electron density gradient fields were considered. This approach identifies the electrostatic contribution to the interactions in solid state.[75,91]

For this purpose the C–H⋯F–C and C–F⋯F–C interactions in the C4HF5 crystal were studied (Fig. 3, upper panel). The superposition of gradient fields of the electrostatic potential and electron density map for C4HF5 crystal is given on the Fig. 4. The C9–H1⋯F9–C5 contact is characterized by significant penetration of v- and ρ-basins of H and F atoms, which reveals the importance of the electrostatic contribution to this interaction. The boundaries of the corresponding basins along the bond path for the intralayer C5–F9⋯F1–C1 interaction (upper panel in Fig. 3) are equivalent, so C5–F9⋯F1–C1 is a typical closed-shell interaction. In case of interlayer F⋯F interactions (C1–F5⋯F5–C1 and C5–F9⋯F13–C13), the boundaries of v-basins nearly coincide with the boundaries of ρ-basins. Therefore, the electrostatic contribution for the interlayer C–F⋯F–C interactions in the C4HF5 crystal is negligible.

Fig. 4

Superposition of gradient fields of the electrostatic potential (pink) and the electron density (blue) in crystalline C4HF5: in the plane of the H1⋯(BCP) ⋯F9 interaction (left), in the plane of the inter-layer F5⋯(BCP)⋯F5 interaction (lower right), in the plane of the intra-layer F1⋯(BCP)⋯F9 (middle right) and the inter-layer F9⋯(BCP)⋯F13 (upper right) interactions. Bond paths are given by black lines. The nuclei positions are given by red circles.

This result is confirmed by a plot of the distribution of the electrostatic potential on a surface of constant density (Fig. 5). The areas corresponded to the interlayer C–F⋯F–C interactions are characterized by a positive electrostatic potential in the isolated molecule. Therefore, these interactions will not be driven by the electrostatic factor in the crystalline structure. However, the fact that the crystalline environment strongly affects the distribution of the electrostatic potential should be taken into account. The distribution of the electrostatic potential for the C–H⋯F–C interaction changes considerably for a molecule removed from the crystal (Fig. 5a and c). This is completely opposite to the interlayer C–F⋯F–C interactions (Fig. 5b and d). Thus, the analysis of the distribution of the electrostatic potential for isolated molecules mainly illustrates the “driving force” in the process of crystal lattice formation.[34] We conclude that this approach should be used with caution when studying the nature of intermolecular interactions in the solid state.

Fig. 5

Electrostatic potential in the C4HF5 molecule, panels (a) and (b), and the molecule removed from the C4HF5 crystal, panels (c) and (d) projected onto surface of constant electron density. The latter value is 0.011 a.u. on panels (a) and (c) and 0.007 a.u. on panel (b) and (d). These values correspond to the electron density at the BCPs for the C9–H1⋯F9 and C1–F5⋯F5 interactions (Table 2). Bond paths are marked with yellow lines. The grey, light grey and yellow circles represent the carbon, hydrogen and fluorine atoms, respectively.

A significant electrostatic contribution to the C–H⋯F–C interaction in crystalline C4HF5 with the H⋯F distance of 2.35 Å, which is shown by the superposition of the gradient fields of the electrostatic potential and the electron density map, is consistent with studies based on the PIXEL scheme. According to these studies, short C–H⋯F–C contacts have a significant Coulomb contribution.[35] In crystals with H⋯F distances around 2.5 Å, the C–H⋯F–C interactions are also driven by the electrostatic factor (Section 3.4). This result differs from the data obtained by gas-phase[33] and PIXEL computations.[36,37] According to these studies, C–H⋯F interactions with H⋯F distances of about 2.5 Å are mainly dispersive. We believe that the reason for the disagreement is the use of non-periodic models and non-relaxed molecular geometry in the PIXEL approach.

To identify the role of C–H⋯F–C interactions in the formation of a crystal architecture, a series of organic crystals consisting of molecules with varying degrees of fluorine substitution were examined. In benzene – partially fluorinated benzenes – hexafluorobenzene series the packing motif changes as follows (Fig. 6): herringbone packing with edge-to-face orientation[113] – layers[23] – edge-to-face orientation.[55] This phenomenon is also observed for naphthalene, anthracene, pyridine (Fig. S1†). This points to the structure-directing role of C–H⋯F–C interactions due to the dominant electrostatic contribution to these contacts. The C–H⋯H–C and C–F⋯F–C interactions do not have that feature. Their energy in crystals is lower than that of C–H⋯F–C, which is manifested in higher values of vaporization enthalpy of partially fluorinated alkanes in comparison with the corresponding fluoroalkanes (Fig. S2 and S4†). The vaporization enthalpies of these crystals are directly proportional to the sublimation enthalpies (Fig. S3†).

Fig. 6

Transformation of packing type in the row benzene – partially fluorinated benzenes – hexafluorobenzene. Refcodes of the crystals are given in parentheses. The C–H⋯F–C interactions are denoted by dotted lines.

We conclude that: (i) the energies of the C–H⋯F–C contacts (5–7 kJ mol−1) are found to be larger than the energies of the C–F⋯F–C contacts (<4 kJ mol−1); (ii) C–H⋯F–C interactions are driven by the electrostatic factor and have the structure-directing role in crystals containing C, H and F atoms.

The competition between C–H⋯O and C–H⋯F–C interactions

One of the CSD analysis tools is the probability density function (PDF). It has a dimension inverse to the measurement of the physical quantity being analyzed, e.g. the H⋯X distance and the CAr–H⋯X angle, here CAr–H denotes aromatic hydrogen, and X = O or F. The integral of the PDF over the entire area of the analyzed physical quantity is equal to one. The maximum value on the PDF curve corresponds to the most frequently encountered H⋯X distance (Fig. 7a) or CAr–H⋯X angle (Fig. 7b). The data on Fig. 7a leads to the following conclusions. (i) The H⋯F/H⋯O distances are determined by the proton acceptor group polarity of the fragment in question. The shortest H⋯O distances are formed by the OC group. (ii) The H⋯F distances of aromatic fluorine (F–CAr) are similar to the H⋯O distances of the OH–C and OC–C groups (Fig. 7a). (iii) Aliphatic fluorine (F–Csp) forms the longest H⋯F distances. It follows from Fig. 7b that C–H⋯F–C/C–H⋯OC contacts are directional if the H⋯F and H⋯O distances are shorter than 2.35 Å (Fig. 7b). If the H⋯F and H⋯O distances are greater than the specified value (2.35 Å < R(H⋯X) < 2.70 Å), the directional character of the studied contacts gradually disappears. 2.70 Å is nearly equal to the sum of van der Waals radii of H and F (∼2.67 Å) or H and O (∼2.75 Å).

Fig. 7

CSD-based estimation of the probability density function (PDF) of the CArH⋯X distance (a) and CAr–H⋯X angle (b) distribution, where X = O, F. The C–H bond lengths were normalized to the neutron-diffraction value (1.089 Å); the cone correction was also considered.[14]

The energies of the ArC–H⋯O and ArC–H⋯F interactions are comparable, if the H⋯F/H⋯O distances are close to each other (Table S6†). Analysis of the nature of the C–H⋯O and C–H⋯F interactions shows their similarity (Fig. 8 and 9). Both interactions are characterized by a significant electrostatic component (Fig. 8). The C–H group is non-polar, as evidenced by the absence of charge depletion on the deformation density maps in the region of this bond (Fig. 9).

Fig. 8

Left panel: superposition of gradient fields of the electrostatic potential (pink) and the electron density (blue) for the KEGWEZ crystal. The bond paths are given by black lines. The nuclei positions are given by red circles. Right panel: the fragment of the KEGWEZ crystal. Atoms forming the C–H⋯F–C, C–F⋯F–C, C–H⋯O interactions are labeled.

Fig. 9

Deformation density maps for the C7–H2⋯O2 contact (left) and C8–H3⋯F1 contact (right) in the KEGWEZ crystal. The solid red lines represent positive contours and the broken blue lines represent negative contours. The values of charge depletion on the hydrogen atoms are given. All the contour lines are drawn at the intervals of ± 0.1 eÅ−3, except the additional contours with a step width of 0.01 eÅ−3 at the values close to those given on the map (0.03–0.1 eÅ−3). For atomic labels see Fig. 8.

The theoretical Elatt values of crystalline KEGWEZ, DFNAPQ and YICBES are given in Table 3. They are in reasonable agreement with the literature data.[52]Eqn (2), which was used to estimate the Elatt values, allows the lattice energy to be divided into different terms.[114,115] The contributions of the C–H⋯O and C–H⋯F–C interactions are comparable and responsible for more than 2/3 of Elatt (Table 3). This means that both interactions can play a significant structure-directing role in organic crystals without conventional H-bonds or other relatively strong intermolecular interactions.

The theoretical Elatt values of crystalline KEGWEZ, DFNAPQ and YICBES evaluated using eqn (2) and relative contributions of various terms to the lattice energy

	KEGWEZ	DFNAPQ	YICBES
E latt, kJ mol−1	70.4	71.6	69.2
Contribution
C–H⋯OC	35%	38%	31%
C–H⋯F–C	34%	28%	42%
C–H⋯pi, C–F⋯pi etc.	31%	34%	28%

The occurrence of the C–F⋯F–C and C–H⋯F–C interactions in isomeric species were analysed for 1-(4-fluorobenzoyl)-3-(isomeric fluorophenyl)thioureas.[97] Their crystal packing shows H⋯F distances in the 2.46–2.82 Å range (Table S7†). In accordance with the results obtained in this subsection, almost no fluorine–fluorine contacts are realized in crystals with the C–H⋯F/C–H⋯O/C–H⋯S interactions.

Discussion

The applicability of Bader's theory to describing the intermolecular interactions was recently questioned.[116] It should be noted that the H-bond energy/enthalpy may be estimated from the spectroscopic[117] and geometric[118] characteristics of H-bonds in crystals. All these approaches give comparable results for molecular crystals with weak or moderate conventional H-bonds.[38,83,87] For example, the H-bond energy in water ices calculated using eqn (3) and the Iogansen approach[117] are consistent with the experimental value of the H-bond enthalpy.[119]Eqn (3) with k = 0.429 yields reasonable energies for intermolecular interactions with significant electrostatic contribution. A change in the nature of the interaction leads to serious differences in the energies calculated by eqn (3) from the estimates made by other approaches. For example, the appearance of a significant covalent component in short (strong) H-bonds, causes a serious variation in the energies calculated using different empirical approaches, see Table 2 in ref. 79 and Table 2 in ref. 120.

The value of the coefficient k in eqn (3) is ∼0.67,[72] 0.429 (ref. 27) and 0.129 for the I⋯I, Cl⋯Cl and F⋯F interactions in crystals, respectively. This trend is consistent with the results of ref. 121.

The fundamental difference between the C–H⋯F–C interactions and the conventional O–H⋯O bonds can be demonstrated by the superposition of the electrostatic potential and electron density gradient fields (Fig. 10). The boundaries of v- and ρ-basins along the C–H bond path overlap perfectly for the C–H⋯F–C contacts, but for O–H, the penetration of v- and ρ-basins is significant. Thus, the O–H covalent interaction is mainly driven by the electrostatic factor. This result is supported by the deformation density maps for C–H⋯F–C and O–H⋯O (Fig. 11). There are regions of charge accumulation and charge depletion in the O–H direction, unlike the C–H case. The lack of charge accumulation near the oxygen atom demonstrates that O–H is a typical polar interaction as opposed to the C–H bond. The presence of the polar O–H bond causes more intensive charge depletion on the hydrogen atom, which leads to stronger O–H⋯O interactions in comparison with C–H⋯F–C and C–H⋯O (Tables S8 and S9†).

Fig. 10

Superposition of gradient fields of the electrostatic potential (pink) and the electron density (blue) for the TISQER crystal. The bond paths are given by black lines. The nuclei positions are given by red circles. For atomic labels see Fig. S6.†

Fig. 11

Deformation density maps for the C8–H3⋯F1 interaction (left) and O1–H8⋯O5 H-bond (right) in the TISQER crystal. The solid red lines represent positive contours and the broken blue lines represent negative contours. The values of charge depletion on the hydrogen atoms are given. The contours are drawn at the intervals of ± 0.05 eÅ−3. For atomic labels see Fig. S6.†

A similar situation arises when we compare the H⋯F energy in C–H⋯F and F–H⋯F fragments. The energy of the F–H⋯F contact (∼23 kJ mol−1 (ref. 122)) is much higher than that of the C–H⋯F ones. This is due to the high polarity of the F–H covalent bond compared to the polarity of the C–H bond. It should be noted, that the X23 (ref. 39) and C60 (ref. 40) benchmark sets for noncovalent interactions in solids do not include molecules with fluorine atoms, while the Z20 set[41] includes only F2 and HF molecules. Obviously, these sets are hardly applicable to adequate description of C–H⋯F interactions.

Conclusions

The PBE-D3 level with the 6-31G** basis set in which diffuse functions are added to the fluorine atoms provides an adequate description of the metric and electron-density parameters of C–F⋯F–C interactions in crystals of perfluorinated compounds, while the AA-CLP approach gives the best values of the sublimation enthalpy of these crystals. It has been found that plane wave periodic DFT computations with the most popular pseudopotentials may fail to describe the crystal structure of fluorine-rich compounds. A quantitative description of the lattice energy of perfluorinated organic molecules implies the use of advanced dispersion energy corrections in the periodic DFT computations or specially parameterized force fields in MD simulations.

Periodic DFT computations followed by Bader analysis of crystalline electron density in combination with thermochemical data enable us to estimate the C–F⋯F–C interaction energy Eint using the following relation: Eint = kGb. Here, Gb is the local electronic kinetic density at the BCP and k = 0.129. This value is much smaller than the coefficient value of 0.429, which is used for evaluating the C–H⋯F–C energy. As a result, the average energy of intermolecular C–F⋯F–C interactions in crystals of neutral organic molecules is <4 kJ mol−1, which is less than the energy of the C–H⋯F–C interactions varying from 5 to 7 kJ mol−1. The presence of C–F⋯F–C interactions in crystals should be verified by Bader analysis of crystalline electron density. This approach gives reliable results for the F⋯F distances less than ∼2.94 Å.

The nature of the considered intermolecular interactions was established through the superposition of the gradient fields of electrostatic potential and electron density. In contrast to weak C–F⋯F–C interactions, intermolecular C–H⋯F–C interactions are driven by the electrostatic factor. The analysis of the Cambridge Structural Database shows the directed nature of the C–H⋯F–C interactions with H⋯F distances less than 2.35 Å. The electrostatic character and partial directionality of C–H⋯F interactions stipulate the structure-directing role of these interactions in crystals containing C, H and F atoms. This is manifested in the differences in the crystalline packing of partially fluorinated aromatic molecules as compared with non-fluorinated or perfluorinated molecules.

Intermolecular C–H⋯F–C and C–H⋯O interactions are identical in nature and comparable in energy. Different factors determine their contribution to the crystal packing: the underlying donor and acceptor properties of the C–H/C–F/CO moieties, the number of C–H⋯F–C and C–H⋯O contacts, etc. C–H⋯O and C–H⋯F interactions play a significant structure-directing role in crystals of neutral organic molecules without conventional H-bonds.

Conflicts of interest

There are no conflicts to declare.