Combined Diffraction and Density Functional Theory Calculations of Halogen-Bonded Cocrystal Monolayers

This work describes the combined use of synchrotron X-ray diffraction and density functional theory (DFT) calculations to understand the cocrystal formation or phase separation in 2D monolayers capable of halogen bonding. The solid monolayer structure of 1,4-diiodobenzene (DIB) has been determined by X-ray synchrotron diffraction. The mixing behavior of DIB with 4,4′-bipyridyl (BPY) has also been studied and interestingly is found to phase-separate rather than form a cocrystal, as observed in the bulk. DFT calculations are used to establish the underlying origin of this interesting behavior. The DFT calculations are demonstrated to agree well with the recently proposed monolayer structure for the cocrystal of BPY and 1,4-diiodotetrafluorobenzene (DITFB) (the perfluorinated analogue of DIB), where halogen bonding has also been identified by diffraction. Here we have calculated an estimate of the halogen bond strength by DFT calculations for the DITFB/BPY cocrystal monolayer, which is found to be ∼20 kJ/mol. Computationally, we find that the nonfluorinated DIB and BPY are not expected to form a halogen-bonded cocrystal in a 2D layer; for this pair of species, phase separation of the components is calculated to be lower energy, in good agreement with the diffraction results.


■ INTRODUCTION
The "halogen bond" is a noncovalent interaction between a halogen atom (typically Br or I) and a Lewis base (typically N, S, or O atoms). This interaction has been reported for a broad range of cocrystal combinations in the bulk. 1−6 Important parallels are often drawn between halogen bonding and hydrogen bonding, as both are strong, robust, and directional interactions. 7,8 It has been reported that the halogen bond can be just as strong as the hydrogen bond and in certain cases can even dominate over hydrogen bonding in the molecular recognition processes. 9 This makes the halogen bond a powerful tool in crystal engineering and explains its increasing use in materials chemistry. In this work, we continue our investigation into the role of the halogen bond in 2D supramolecular networks and its potential to control selfassembly in physisorbed layers. 10−13 The halogen bond can be thought of as arising from the electrostatic interaction between a lone pair of electrons and a region of positive molecular electrostatic potential that forms at the tip of a halogen atom. 14−17 This region of positive potential is termed the "σ-hole". For organic halides, the σ-hole is affected by substituents attached to the carbon backbone of the molecule. Electron-withdrawing substituents such as fluorine are considered to result in a more positive σ-hole and hence a greater interaction strength. 16 We have recently reported the formation of a halogenbonded 1:1 stoichiometry cocrystal between 4,4′-bipyridine (BPY) and 1,4-diiodotetrafluorobenzene (DITFB) in a monolayer physisorbed on a graphite surface 10 (molecular structures are illustrated in Figure 1). In the BPY/DITFB colayer structure (see figure 2b in ref 10), the molecules form extended chains of alternating DITFB and BPY. There is evidence of a strong halogen bond between the iodine atoms of DITFB and the nitrogen atoms of BPY, as deduced from the internuclear separation (2.84 Å), which is shorter than the sum of the van der Waals radii of the species. This is in good agreement with the bulk behavior, where halogen bond formation is also observed. 18 Interestingly, there are reports of a similar halogen bonded cocrystal of BPY with the nonfluorinated DITFB analogue (DIB) in the bulk. 18 This halogen bond is reported to be weaker than for the fluorinated case, which is consistent with the current understanding of halogen bonding. 16 The experimental evidence of a weaker bond is supported by DFT calculations on closely related systems. 19 Despite the weaker strength of the bond formed by the nonfluorinated DIB, this interaction still appears to be moderately robust: a halogen bond has also been reported to form between DIB and 1,2bis(4-pyridyl)ethylene, a pyridine-containing molecule very similar to BPY. 18 The study of physisorbed supramolecular networks is of interest due to their applicability in a range of industrial processes. In particular, much work has been devoted to the study of noncovalent interactions that can control self-assembly in physisorbed overlayers. 20,21 In this work, we seek to understand the details of the DITFB/BPY halogen bond in the monolayer and to consider the similar interaction of the nonfluorinated analogue, DIB, with BPY, and identify any differences between the monolayer and bulk behavior.

■ EXPERIMENTAL METHOD
The experimental synchrotron X-ray scattering was performed on Beamline I11 at Diamond Light Source, U.K. 22 The X-ray wavelength used was 1.054700 Å with a detector offset of 0.05899°, as determined by Rietveld refinement of a silicon standard (NIST SRM 640c).
The graphite substrate used was Papyex (>99.6% carbon), a recompressed exfoliated graphite with a specific surface area of 27.8 m 2 g −1 , as determined by nitrogen adsorption. The graphite was cleaned prior to use by outgassing any adsorbed molecules at a pressure of 0.01 mbar and a temperature of 625 K for 3 h. The adsorbates were purchased from Sigma-Aldrich with purities stated on the certificate of analysis as follows: DIB, purity 99.8% by HPLC; BPY, purity 99.9% by HPLC; DITFB, and purity 99.3% by GC. These were used without further purification.
The graphite was dosed with the adsorbates of interest through the vapor phase. The appropriate amount of adsorbate was weighed into a glass tube containing the cleaned graphite. The tube was then evacuated to a pressure of ca. 0.1 mbar and sealed closed. The sample was annealed by heating the sealed glass tube to 470 K for 2 h 30 min, before allowing it to cool slowly over the course of ∼10 h. The dosed graphite was then cut into 3 mm diameter discs that were stacked into a Lindeman capillary and sealed with a flame. The graphite crystallites in Papyex have a preferred orientation, which was used to optimize the collection of scattering from the monolayer by stacking the graphite so that the plane of preferred orientation is aligned with the scattering plane.
The capillaries were rotated on the diffractometer to enhance powder averaging, and the pattern was recorded simultaneously over the angular range 1−90°using a position-sensitive multidetector. 23 Diffraction patterns were collected from the monolayer adsorbed on the graphite substrate and from the graphite substrate alone. Subtraction of one from the other gives the scattering from the monolayer. The sample temperature was controlled with a cryostream (Oxford Cryostreams, U.K.).

■ COMPUTATIONAL METHOD
The periodic boundary conditions DFT code CASTEP 24 was used to optimize the lattice parameters for the pure BPY and DIB overlayer crystals and for the colayers. Given the relative chemical inertness of the graphitic substrate and the flatness of the potential energy surface suggested by the experimental results, we have modeled the three selfassembled systems without explicitly considering the surface− adsorbate interactions. We used the Perdew Burke Ernzerhof 25 exchange-correlation functional with a kinetic energy cutoff for the plane-wave basis of 340 eV. Long-range intermolecular interactions are not accounted for by traditional DFT in either the LDA or GGA formalism, so to estimate the total binding energy between the aromatic molecules forming the supramolecular network, we applied the TS dispersion force correction method developed by Tkatchenko and Scheffler, 26 in which the C 6 interatomic coefficients are derived from the electron density of the molecular system. The quality of the TS corrections for surface calculations and supramolecular selfassembly calculations is generally robust. 13,27,28 Dispersion force corrections, rather than nonlocal functionals, were used because in the plane-wave DFT code employed in this work (CASTEP), cell structure optimizations cannot be performed with fully nonlocal functionals. This is because there is no known analytical expression for the stress tensor components in the functional formalism. Furthermore, in CASTEP, nonlocal functionals are orders of magnitude more expensive in terms of computational time.
During the geometry optimizations, the forces are converged with a tolerance of 0.05 eV/Å, while the electronic energy tolerance is set to 10 −6 eV. In these calculations, the molecular structure, the unit cell dimensions, and the lattice symmetry have not been constrained.
■ RESULTS: SYNCHROTRON X-RAY DIFFRACTION DIB Overlayer. The experimental diffraction pattern of DIB overlayer recorded at a coverage of 0.512 equivalent monolayers (MLs), and a temperature of 200−230 K is shown in Figure 2a in gray. Scattering from the graphite substrate and small-angle "Porod" scattering have been subtracted 29 so that the final pattern shows only diffracted intensity from the physisorbed layer. These peaks have the characteristic sawtooth shape of diffraction from 2D layers. Hence we can conclude that DIB forms a solid crystalline overlayer at this temperature and coverage. Imperfect subtraction of the strong (0002) peak from the graphite substrate at 1.9 Å −1 is evident and will obscure any scattering from the overlayer in this region; data from this region have therefore not been included in the Figure. The imperfect subtraction of this graphite feature also means that the relative intensity of the neighboring overlayer diffraction peak at 1.79 Å −1 is less reliable than that in the other peaks.
The process of analyzing ML diffraction data has been described in detail elsewhere. 30 In brief, because of the small number of X-ray reflections present in the accessible range of momentum transfer, any fit to the data must be constrained as much as possible. Hence, the structure of the DIB molecule used in the fitting process was taken unchanged from the 3D crystal structure (Cambridge Crystallographic Database refcode ZZZPRO05), and only rigid body rotations and translations of this molecule were considered. In addition, high-symmetry plane groups with fewer degrees of freedom were considered before lower symmetry plane groups.
The Gaussian, Lorentzian, and Lorentzian-squared peak shape models of Schildberg and Lauter were considered in this work. 31 The Lorentzian-squared model has been used for the final fit, as this provides the closest approximation to the experimental line shape. This model includes terms for the size and preferential orientation of the graphite crystallites, which were fitted to the experimental data. A single temperature factor set to unity was used.
The experimental pattern was indexed with a rectangular unit cell of dimensions a = 16.38 Å, b = 6.81 Å, and ν = 90°. The (h, 0) and (0, k) reflections are systematically absent from the experimental pattern for odd h and k. The rectangular symmetry of the lattice and the presence of (h + k) ∈ 2Ζ+1 (odd) reflections imply that the unit cell has p2gg symmetry. 32 A unit cell of this size can contain two close-packed DIB molecules adsorbed with the ring flat on the graphite substrate. When flat, the DIB molecules possess a two-fold rotation axis. Because of the size and shape of the molecules, for them to pack in a physically reasonable manner without overlap, the two-fold rotation axes of the molecules must coincide with the rotation axes of the unit cell. This means that the molecules have no translational freedom and must be centered at the origin and (1/2, 1/2).
If the molecules are constrained to lie in a plane parallel to the graphite substrate, then the only degree of freedom is rotation about the surface normal (the z axis). Because the two molecules in the unit cell are related by symmetry, the rotation of one molecule determines the rotation of the other. Hence, there is only one freely variable parameter to be fitted. This parameter is also constrained to a reasonably narrow range of angles to avoid overlap with adjacent molecules. A fit to the experimental pattern based on this highly constrained model produces a good match to the experimental data.
The fit can be slightly improved by relaxing the constraint that the molecules must lie parallel to the surface. This introduces two extra parameters for each molecule: rotation about the C 2 axis passing through the two iodine atoms (I−I axis) and rotation about a second C 2 axis perpendicular to the first that passes through the center of mass of the molecule (referred to as the C 2 ′ axis for convenience). The pattern is largely insensitive to the structure in the z direction. As molecules are tilted up from the surface, the majority of the change in the pattern results from the change in the x and y coordinates of atoms projected onto the plane rather than the change in z coordinate of the atom per se. This means that the fitting procedure is relatively insensitive to small rotations, where the x and y coordinates do not change significantly.
The best-fit structure does not show any significant rotation about the I−I axis. However, we stress that the fitting procedure is largely insensitive to rotation about this axis, with an uncertainty of ca. ±15°. The calculated pattern is somewhat more sensitive to rotation about the C 2 ′ axis because rotation about this axis results in a change in x and y coordinates of the strongly scattering iodine atoms. The best-fit structure has a rotation about the C 2 ′ axis of 7 ± 5°. Because the overlayer is now no longer perfectly flat, the unit cell symmetry is strictly no longer p2gg. Instead, the rotation axes become centers of inversion, and the glide lines become either glide planes or screw axes. Because each of the two molecules in the unit cell can be rotated independently about the C 2 ′ axis, this leads to two possible structures, depending on whether the second molecule in the unit cell has the same or opposite sense of rotation to the first. We have named these possibilities "glide" and "screw", respectively (based on the symmetry element parallel to the a axis; see Figure 2b). It is not possible to distinguish between these two structures based on the experimental diffraction data.
The agreement between experimental and calculated patterns can be measured with several parameters. Here we have calculated the "goodness of fit" parameter, R, which is given by: where N is the number of data points and n is the number of fitted parameters. However, as discussed elsewhere, it is not possible to compare these values directly with similar calculation on bulk diffraction patterns. 12 The best fit has R = 0.14 and χ 2 red = 115 and is shown in black in Figure 2a. The corresponding glide/screw structure is shown in Figure 2b. Given the highly constrained nature of the model, the fit to the experimental data is still reasonably good. The small area of this highly symmetrical unit cell allowed us to Langmuir Article dx.doi.org/10.1021/la402910a | Langmuir 2013, 29, 14903−14911 perform a search of all possible molecular rotations assuming the above p2gg symmetry cell with two molecules at fixed positions. From this, we conclude that this is the only p2gg symmetry structure that fits the experimental data.
The DIB diffraction pattern was also collected at lower temperatures down to 100 K. On cooling below 160 K, several of the diffraction peaks are observed to split. This is interpreted as a reduction of symmetry of the overlayer lattice arising from a change of the unit cell angle from 90°. However, to assign all of the split reflections, a unit cell doubled in both the a and b directions is required. This gives rise to a rather large unit cell containing at least eight molecules. The total number of parameters to fit is now much higher (as all of these eight molecules can be varied independently given the low plane group symmetry). Hence we are unable to uniquely determine this new low-temperature structure on the basis of diffraction data alone.
Mixed DIB/BPY Overlayer. The formation of the cocrystal overlayer of DITFB and BPY has been previously reported by us using synchrotron diffraction. 10 Here we address the behavior of the DIB and BPY overlayer using the same approach. In addition, we aim to use DFT calculations to understand the cocrystal formation or phase separation of both these binary mixtures. Figure 3 presents the diffraction data from the pure DIB overlayer, the pure BPY ML, and a 1:1 mixture of DIB and BPY all at the same coverage of 0.5 ML. The graphite background has been removed in preparing these figures, and so the   scattering represents that from the overlayers alone. We note that in all cases the peaks have the very characteristic sawtooth line shape indicative of overlayer diffraction, confirming the formation of solid overlayers of the adsorbates.
All of the peaks in the pattern of the 1:1 mixture are found in the pattern of either DIB or BPY, typical of phase separation of the two components. Significantly, there are no new peaks, which would be expected if a cocrystal was formed. Figure 3 also shows the pattern obtained by summing the patterns of DIB and BPY, which is the expected pattern if the two components phase separate. This pattern agrees very well with the experimentally observed pattern of the 1:1 mixture. Hence we conclude that the 1:1 mixture of DIB/BPY does not form a cocrystal but phase separates on the surface. This is in marked contrast with the DITFB/BPY combination and rather surprising as the DIB/BPY combination forms a cocrystal in the bulk.
The conclusion that phase separation is occurring is further supported by the change in the diffraction pattern with temperature ( Figure 4). ML diffraction patterns of the 1:1 system were recorded between 200 and 300 K. At ∼244 K, several peaks disappear from the pattern of the 1:1 mixture, all of which are found in the pure DIB pattern. Hence at 244 K, we identify melting of a pure DIB phase to leave behind a second solid phase. The remaining peaks at temperatures greater than 244 K, illustrated in Figure 5, are consistent with a solid overlayer phase of only BPY.
■ RESULTS: DFT CALCULATIONS DIB. Figure 6 shows the unit cell and electron density difference obtained by DFT for the pure DIB overlayer. This Figure is obtained by subtracting the calculated electron density for the isolated DIB molecule from the electron density calculated for the molecules in the 2D overlayer.
Calculations indicate that the screw symmetry structure is only marginally more favorable than the glide (+5 meV) and the planar geometry (+26 meV). The optimized cell parameters for pure DIB in the screw symmetry arrangement are very close to the experimental results (+0.47% a, +2.5% b, +0.24% ν), therefore suggesting that the interactions with the graphite surface (neglected in the present computational work) do not influence significantly the equilibrium lattice structure of the DIB overlayer. The van der Waals forces dominate the intermolecular interactions with a total contribution of 0.5 eV per cell, almost 90% of the total binding energy per cell. The covalent and directional (hydrogen-and halogen-bonding) bonding is accounted for by standard DFT; therefore, we calculate the binding energy coming from intermolecular H bonding by performing calculations without vdW corrections (DFT+D). The difference in the total binding energy per cell obtained by DFT+D and standard DFT calculations provides us with the intensity of long-range (vdW) dispersion interaction. 13 The charge accumulation and depletion upon the formation of the overlayer from isolated DIB monomers can be seen in Figure 6. In this Figure, it is evident that there is rather little concentration of negative charge between the iodine of one molecule and closest hydrogen atom of the neighbor DIB. This indicates a low degree of directional bonding.
BPY. Figure 7 presents the electron density difference for BPY. This Figure is obtained by subtracting the calculated electron density for the isolated BPY molecule from the electron density calculated for the molecules in the 2D overlayer.
This simulation was initially constrained to have the unit cell of the experimentally determined overlayer structure. 11 However, removing this initial constraint resulted in <0.6% change (+0.53% a, −0.44% b, +0.01% ν) in lattice parameters. Hence we conclude that the experimentally determined overlayer structure and that determined by DFT are in good agreement. Interestingly, we also conclude that the effect of the graphite on the overlayer structure is again very small, given that the graphite substrate was not included in the simulations.  The high symmetry of the overlayer structure means that the packing is essentially dominated by the H bonding between the nitrogen of one ring interacting with two hydrogen atoms on the other molecule. There are eight of these interactions in the overlayer. The energy of the H-bonding interaction is calculated to be 0.482 eV, and hence we conclude that each bond is ∼60 meV. This is a reasonable value for the weak hydrogen bond expected in this configuration. The total binding energy of the overlayer, accounting for vdW type interactions, is −0.873 eV. (The H bonds account for 55% of the energy gain.) BPY/DITFB. Figure 8a shows the electron density difference map for the DITFB/BPY cocrystal. Again, comparison between the experimentally determined lattice parameters and those determined by DFT calculations agree well, within 2.3% (−1.36% a, −2.24% b, −1.45% ν). This again indicates that the X-ray diffraction and DFT calculations agree well on the overlayer structure and that the role of the graphite is minimal. The structure is characterized by lines of alternating ...BPY-DITFB-BPY-DITFB... molecules with halogen bonds between the species in the chain, as inferred from X-ray diffraction by the short internuclear separation of the N and I atoms.
The energy of this structure is found to have noncovalent (including vdW dispersion forces) interactions of ∼1.078 eV (104 kJ/mol). There are two weak hydrogen bonds on either side between the chains of molecules. By calculating the binding energy of a linear chain of DITFB-BPY molecules separated from any other chain, we can estimate the contribution to the binding energy coming from the halogenbonding alone. This is found to be 0.249 eV per bond typical of a strong halogen bond in the bulk. 33 The energy of the interchain hydrogen bonding is then estimated to correspond to ∼0.062 eV each or a total of 0.248 eV.
The electron density difference map reveals an interesting topological distribution. Around the iodine atoms is a torus of electron density accumulation, while at the top end toward the nitrogen atom there is a region of electron deficiency. This is a clear topological representation of the "σ-hole" previously suggested to be involved in halogen bonding. Furthermore, the nitrogen lone pair is evident as a region of electron density that fits into the iodine sigma hole. (See Figure 8b.) In conclusion, this DFT calculation agrees well with the experimental observation that there is strong halogen bond formation in this binary combination of BPY/DITFB and that a cocrystal is favored.
BPY/DIB. Because there is no experimental evidence of the cocrystal of BPY/DIB, we need to make a guess for the most likely cocrystal structure, if it did form. In doing this, we have taken the BPY/DIB cocrystal to be isomorphic with the BPY/ DITFB cocrystal above.
The electron density difference for the hypothetical BPY/ DIB cocrystal, isomorphic to the BPY/DITFB combination above, is given in Figure 9.
The purpose of these calculations is to understand why the BPY/DIB colayer is not observed in the overlayer, and instead the two species separate into two distinct ML phases (pure DIB and pure BPY overlayers).
The results of the DFT calculations show a substantial decrease in the total directional (non-vdW) bonding from 0.75 to 0.36 eV (52% decrease from BPY/DITFB to BPY/DIB). We can conclude that although DFT predicts a total favorable  interaction between DIB and BPY the strong decrease in hydrogen and halogen bonding in the hypothetical cocrystal hints at a relatively high instability of the overlayer under experimental conditions.

■ DISCUSSION AND CONCLUSIONS
For the pure DIB overlayer, the DFT calculations suggest that of the possible structures determined by X-ray diffraction the lowest energy structure is the screw-symmetry structure. However, with such small energy differences separating the "screw" structure from the "glide" and "flat" arrangements (5 and 26 meV respectively), it is likely that the experimental structure is a statistical average of these structures at the experimental temperature of 200 K. It is often the case that the 2D overlayer structure resembles a plane from the 3D "bulk" structure. 12,13 Interestingly, the two proposed lowest energy structures also closely resemble the same plane from the two observed bulk polymorphs. The "glide" structure is very similar to the (010) plane of the low-temperature alpha polymorph, and the "screw" structure is very similar to the (010) plane of the high-temperature beta polymorph, which is shown in Figure  10. 34 The zigzag arrangement of iodine atoms closely resembles the placement of iodine atoms in iodoalkane overlayers. 35,36 In the flat arrangement, where the iodine atoms approach each other most closely, the I−I separation is 4.36 Å, somewhat larger than the sum of the van der Waals radii, 3.96 Å. This implies that any interaction between iodine atoms is rather weak, in contrast with the strong iodine−iodine interactions that have been observed for iodoalkanes, where iodine atoms on neighboring molecules are positioned closer than the sum of their van der Waals radii. 35 The lack of strong interactions is supported by the DFT calculations, which indicate that the vast majority of the binding energy of the layer is due to dispersion interactions. However, the molecules are oriented so that the σhole of one iodine atom would point at the belt of electron density surrounding the circumference of the neighboring iodine atom. The DFT calculations ( Figure 6) suggest a modest change in electron density around the iodine atoms, indicative of some interaction between iodine atoms in addition to dispersion alone.
The combined results of synchrotron X-ray diffraction and DFT calculations have been used to address halogen bond formation in overlayer cocrystals of BPY and iodobenzenes. When the iodobenzene species is fluorinated there is evidence of strong halogen bond formation in 2D layers (DITFB/BPY). However, when the iodobenzene species is not fluorinated no halogen bond is formed and the two species (DIB/BPY) phase separate on the surface. This is very different to the bulk behavior of the DIB/BPY mixture, which forms a cocrystal. 18 The DFT results suggest that this phase separation arises due to the much lower binding energy between the DIB and the BPY molecules in the BPY/DIB colayer. Because DFT provides us with an accurate estimate of the amount of hydrogen and halogen bonding energy per unit cell, as well as the dispersion (vdW) interaction between neighboring molecules, we can also derive from these quantities the total energy change (driving force) for forming (or breaking) the colayer. The binding energy per cell needs to be corrected for the change in the surface density due to compression or expansion of the surface cell parameters when going from the mixed phase (BPY/DIB) to phase-separated DIB + BPY. This can be accomplished by calculating the binding energy per unit area, γ (obtained by dividing the total binding energy by the area of the unit cell). The results are summarized in Table 1.
The result of the calculations shows that the BPY/DIB colayer formation from the two separated DIB and BPY phases is not favorable (Δγ > 0). We therefore conclude that the lack of formation of a cocrystal phase of BPY/DIB is that the total intermolecular binding energy is maximized when the two components are separated rather than in the cocrystal. The cocrystal is destabilized in the DIB with respect to the DITFB cocrystal due to both weaker halogen bonding and weaker interchain interactions.
It would be convenient to perform a similar calculation for the overlayer DITFB/BPY combination. Unfortunately, the structure of the pure DITFB overlayer crystal on graphite is not available. However, the bulk structures of DITFB, BPY and the cocrystal are known and hence provide some verification of the validity of the first-principle methods applied in this work. Hence we have calculated the formation energies of the bulk cocrystals BPY/DIB and BPY/DITFB from the respective bulk crystal BPY, DIB, and DITFB. The pure bulk crystal structures of BPY, 37 DIB, 34 and DITFB 38 contains two, four, and two molecules per unit cell respectively, while the BPY/DIB 18 and