Halogen Bonding and Cooperative Effects in Chlorine Clathrate: Ab

Jun 27, 2019 - The possibility of forming halogen bonds inside the clathrate cages is an intriguing question. On the basis of X-ray measurements, a pr...
1 downloads 0 Views 13MB Size
Download PDF
Article Cite This: J. Phys. Chem. C XXXX, XXX, XXX−XXX

pubs.acs.org/JPCC

Halogen Bonding and Cooperative Effects in Chlorine Clathrate: Ab Initio Periodic Study Cristina Cuautli and Ramoń Hernań dez-Lamoneda* Centro de Investigaciones Químicas CIQ-IICBA, Universidad Autónoma del Estado de Morelos, Av. Universidad 1001, Cuernavaca 62209, Morelos, Mexico

Downloaded via UNIV OF SOUTHERN INDIANA on July 31, 2019 at 06:53:27 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: The possibility of forming halogen bonds inside the clathrate cages is an intriguing question. On the basis of X-ray measurements, a proposition of “multidirectional halogen bonding” was made, and subsequently, ab initio local correlation studies for dichlorine in a single dodecahedral cage showed evidence of its presence. In this work, the structure and energetics of dichlorine clathrate in its cubic structure I (sI) is analyzed by means of first-principles ab initio periodic calculations based on the density functional theory. We use functionals that are able to properly represent noncovalent interactions particularly halogen bonding. The structural characteristics of the crystal are very well reproduced. Compared with single cage calculations, one important effect of treating the whole crystal is a reduction of the regions with negative values of the electrostatic potential, which reflects in weaker inclusion energies for the crystal. Applying a variety of methods to characterize the electron density, we analyze the guest−host interactions in detail confirming the clear presence of halogen bonding in the dodecahedral cage. The case of the tetrakaidecahedral cage, the largest one, is more subtle, and we show that its quantitative treatment depends on allowing occupancy of several cages in the unit cell giving rise to a strengthening of the halogen bond through cooperative effects.

■

INTRODUCTION In 2013, a detailed X-ray study of the chlorine and bromine hydrate clathrates was performed.1 In the case of chlorine, the crystal structure is of the cubic type usually referred to as structure I (sI), with a unit cell comprising six tetrakaidecahedral cages (51262 or large) and two dodecahedral cages (512 or small). It proved that dichlorine could occupy the small dodecahedral cage and that the percentage occupancy was 32.5%. This result clearly indicated that the guest−host interactions should go beyond the weak van der Waals regime since the free diameter of the cages is estimated at 5.0 Å and the van der Waals size of dichlorine is 5.5 Å.1 This point was further confirmed by the shortest Cl−O distances measured, that is, 2.9−2.97 Å, which are below the estimate from van der Waals radii: 3.27 Å. The larger cage has an estimated free diameter of 5.8 Å and thus a looser fit for the dichlorine guest. The measured shortest Cl−O distances fall in the range of 3.1−3.2 Å, still below the estimated van der Waals radii but not by much, thus requiring a more thorough analysis to clarify the nature of the guest−host interactions. The dichlorine internuclear distance could give valuable information regarding the intermolecular forces; however, the authors report large librational shortening effects, which preclude a reliable estimate of this quantity. On the basis of their data, the authors conclude that halogen bonding (XB) interactions are present and refer to them as “multidirectional halogen bonding” to explain the delocalization observed mainly for © XXXX American Chemical Society

the dibromine case drawing an analogy to previous studies on guests that form hydrogen bonding interactions with the cages.2−9 The dihalogen clathrates have also been studied using a variety of spectroscopic techniques, which have largely increased our understanding of their electronic and vibrational characteristics.10−13 Among the spectroscopic signatures stand out large blue-shifts (400 − 900 cm−1) carefully characterized in the case of bromine clathrate for different phases.14 Here, again, the blue-shifts could be associated with XB interactions, but the observed values are significantly smaller than those calculated with high-level ab initio methods15 for the X2−H2O (X = Cl, Br) complexes: 1600 and 2000 cm−1, respectively. This is not surprising since the water lone pairs required for the XB are engaged in forming hydrogen bonds to hold the cage structure, and thus, their participation in XB is compromised. One option is to disrupt the hydrogen bonding pattern that holds the cages, but it does not seem energetically feasible. So, the question of whether halogen bonding can be present is difficult to decide exclusively on the basis of the experimental measurements performed so far. In 2016, our group performed ab initio local correlation calculations to analyze the nature of guest−host interactions Received: April 8, 2019 Revised: June 26, 2019 Published: June 27, 2019 A

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C for dihalogens in clathrate cages.16 Using energy partitioning analysis,17,18 we could quantify the most relevant components: electrostatic, exchange repulsion, and dispersion. Additionally, the so-called ionic term was seen to correlate with halogen bonding reflected in the dihalogen polarization and small charge transfer from the cage. Combined with estimates of the electronic shifts, it was possible to conclude that, for dichlorine in the small cage, there is clear evidence of XB interactions. More recently, applying real space analysis of the electron density through the interacting quantum atoms approach,19 it has been further confirmed that XB is present for dichlorine in the dodecahedral cage.20 It was shown that the oxygen lone pairs involved in hydrogen bonding can simultaneously participate in halogen bonding so that it is not necessary to disrupt the hydrogen bonding pattern of the cage for this type of guest−host interaction. In this work, we perform ab initio periodic calculations on the dichlorine clathrate to analyze the guest−host interactions and overall energetic stability of the system. Contrary to all previous studies,16,21,22 we take into account the effect of the whole crystal as well as the possible guest occupation of multiple cages. These allows us to conclude the validity of using single cage models. We apply electron density analysis techniques to characterize the nature of the guest−host interactions and to conclude the presence of halogen bonding interactions. Additionally, we shed light on previous experimental results and interpretations.

Figure 1. Left: coordinate system used for the rigid scan. Right: sI unit cell. The red circles indicate the centers of the 512 and 51262 cages where the dichlorine is inserted. Color code: oxygen atoms in red and hydrogen atoms in gray.

the Cl2 were allowed to move with a convergence criteria of 9.52 × 10−4 Å (1.8 × 10−3 au) for the largest atomic displacement and a threshold on the energy change between optimization steps of 10−7 au. In each system, the inclusion energy (I) was calculated as I = EsI − EsIe − ECl 2

(1)

where EsI is the energy of the sI structure with the chlorine molecule in the 512 or 51262 cage, EsIe is the energy of the sI structure where all cages are empty, and ECl2 is the energy of the chlorine molecule. After characterizing the stationary points in the potential energy surface, the guest−host interactions were analyzed through the electrostatic potential, NCI index, ELF, Hirshfeld charges, and charge density differences as detailed below. Computational Details. All calculations were performed with periodic conditions, DFT was used in the Kohn−Sham approximation with the PBE0 hybrid functional24 + D2 (Grimme) dispersion correction25 and POB-TZVP26 localized basis set as implemented in the CRYSTAL 14 code.27 The PBE0-D2 functional was chosen because it reproduces both the angular dependence and strength of the inclusion energies of the potential energy surface (PES) of the isolated 512 cage calculated with the local DF-LMP2 method and an AVDZ/ AVTZ basis set for the water/chlorine molecules.16 Here, only the Γ point was used. The inclusion energy was calculated using the counterpoise method in order to correct the basis set superposition error. The NCI analysis was made with the NCIplot code.28 NCI. The NCI index29,30 has been used in order to analyze noncovalent interactions such as hydrogen and halogen bonds.31 The NCI is based on the calculation of the reduced density gradient (RDG)

■

METHODOLOGY The cell parameters and the atomic position of the oxygen atoms were obtained from a recent X-ray study.1 The hydrogen positions were obtained from a previous periodic study23 where the orientation of the protons was selected since they give rise to a nearly zero net moment dipole of the unit cell and have the lowest potential energy from among 685,686,200 geometries. In the structure thus obtained, the hydrogen atom positions were optimized without the halogen in any cage maintaining the oxygen atoms and the cell parameters fixed. Afterward, two calculation series were performed: One chlorine molecule was placed in the center of (a) one 512 cage and (b) one 51262 cage. The other cages in the unit cell were empty. First, a rigid angular scan of the potential energy was made with a constant Cl−Cl distance in order to obtain an overall view of the potential energy surface efficiently; afterward, selected relevant structures were partially optimized in order to analyze the halogen−water interactions in detail. Rigid Angular Scan. The dihalogen was positioned in the center of the cage under analysis, and the distance Cl−Cl was fixed at 2.00 Å. The rotation of the dihalogen molecule was carried out by the use of spherical coordinates: φ is the azimuthal angle, and θ is the polar angle. The origin of coordinates is in the center of the cage and coincides with the center of mass of the dihalogen molecule. As the unit cell of the sI crystal is cubic, the x, y, and z axes coincide with the a, b, and c unit cell parameters. The internal coordinates are 0, 0, 0 and 0, 0.5, 0.25 for the center of the 512 and 51262 cages, respectively (Figure 1). The angle of variation for each scan was made 15°. Partial Optimizations. The global minima of the potential energy surface obtained with each rigid angular scan was partially relaxed: the cell parameters and oxygen atoms were kept frozen during the optimization. The hydrogen atoms and

RDG =

|∇ρ| 1 2 1/3 4/3 2(3π ) ρ

(2)

that accounts for the local inhomogeneity of the electron density ρ. The RDG takes small values for small gradients, typical of bound regions. Plotting RDG versus sign(λ2)ρ (the second eigenvalue sign of the electron density Hessian matrix multiplied by ρ) for low values of ρ shows peaks where there are noncovalent interactions. The nature of the interactions can be differentiated: the attractive interactions appear at negative values, and the repulsive interactions appear at B

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 2. (a)Angular dependence of the PES for Cl2 @ 512, (b) electrostatic potential mapped into a charge density isosurface (isovalue 0.006 au), (c) cut of the electrostatic potential in ϕ = 135°, and (d) atomic structure of the 512 cage.

the potential is anisotropic and there are preferred orientations. The energies of the global minimum and maximum are −3.89 and −1.19 kcal/mol, respectively. The electrostatic potential for the sI structure without the guest was calculated to examine the electrostatic component of the halogen−water interaction. Figure 2b shows the electrostatic potential of the sI structure mapped into a charge density isosurface (isovalue of 0.006 au). A cut in one plane that splits the 512 cage in half was made to observe the electrostatic potential inside the cage (this is delimited by the dotted line). The other side of the potential is symmetrical by an inversion symmetry operation. Figure 2b shows that the electrostatic potential is mostly positive with neutral and slightly negative regions around the oxygen atoms. A correlation between the electrostatic potential and the geometry of minima and maxima of the PES is shown in Figure 2a,b by labeling with numbers and letters; in particular, one can correlate the inclusion energies with the place where the dichlorine points to in the electrostatic potential. In Table 1, the inclusion energies for the different stationary points are shown. In the structures with the more negative inclusion energies, the dihalogen is oriented toward negative areas of the electrostatic potential allowing an attractive electrostatic interaction with the positive region of the halogen molecule

positive values of sign(λ2)ρ. The stronger the interaction, the higher the value of ρ. By plotting RDG isosurfaces at low values (∼0.4) around the peaks, the atoms involved in each interaction can be identified. ELF. The electron localization function (ELF)32 is a measure of the spatial localization of the electron pairs. It is expressed as 1

ELF = 1+

2

( ) Dσ

Dσ0

(3)

where Dσ is the electronic localization of the system under study and D0σ is the electronic localization of the homogeneous electron gas. The ELF can take values between 0 and 1. The ELF value is near to 1 within the pair region, is small when close to the border between two pair regions, and is near to 0.5 in regions with electron localization like the homogeneous electron gas. The ELF has been successfully used in order to distinguish single, double, and triple covalent bonds in molecules and ionic and metallic bonds in solids33 and to analyze the hydrogen bond.34 Results. 512 Cage. The potential energy surface (PES), that is, the inclusion energy as a function of the guest spherical coordinates (ϕ and θ angles), obtained in a rigid angular scan, is shown in Figure 2a. Given the estimates of the Cl2 van der Waals size (5.5 Å) and of the cage effective diameter (5.0 Å), it is interesting to notice that the inclusion energy is negative for all orientations. As a matter of fact, in the original X-ray determination by Pauling and Marsh,35 the authors assumed that the dodecahedral cages would be empty. The occupation of the 512 cages by Cl2 has been proven more recently through X-ray analysis.1 Notice the presence of several minima with low energy barriers between them (roughly 2 kcal/mol). Although the minima cover a large range of angular space, it is clear that

Table 1. Inclusion Energies (kcal/mol) for the Stationary Points Shown in Figure 2a,b point I point I C

a −3.89 1 −1.19

b −3.80 2 −1.41

c −3.78 3 −1.56

d −3.38 4 −1.61

e −3.36

f −3.36

g −3.33

h −3.21

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 3. Geometries of the analyzed structures for the 512 cage: (a) trimer Cl2 − 2H2O, (b) partially optimized minimum, (c) rigid scan minimum, and (d) rigid scan maximum.

Table 2. Inclusion Energies, Geometrical Parameters, and Hirshfeld Charges in Cl2 for the Analyzed Geometries of Cl2@512a structure

I

Cl−O

Cl−Cl−O

Cl−Cl

Cl2 Hirshfeld charge

trimer partially optimized minimum rigid scan minimum rigid scan maximum

−5.74 −5.42 −3.89 −1.19

2.79 2.76 2.83 3.00, 3.13, 3.19, 3.25, 3.25

172.12 176.47 166.84 138.61, 139.25, 121.88, 128.32, 127.94

2.125 2.105 2.00 2.00

−0.06 0.13 0.14 0.17

a

The energies are in kcal/mol, the distances are in Å, and the angles are in degrees.

the crystal. The PES and the electrostatic potential of the extracted 512 cage are given in Figure S2. A deeper analysis of the halogen−water interaction can be made by describing in detail the main stationary points and other regions of the PES. We also focus on a particular fragment: the trimer (Cl2 − 2H2O), which is ubiquitous for the most stable structures. The geometries selected are shown in Figure 3, and their inclusion energies and main geometrical parameters of Cl2 and the closest H2O molecules are shown in Table 2. The minima and maxima geometries are a and 1 from the the rigid scan shown in Figure 2a. The global minimum of the rigid scan was partially optimized as detailed in the Methodology, Partial Optimizations section and will be denoted as partially optimized minimum. Since the inclusion energies of structures b and c (Figure 2a) are very close to the global minimum of the rigid scan PES (a), these structures were also partially optimized. The inclusion energies of the b and c partially optimized structures are −4.05 and −4.35 kcal/mol, respectively, both still higher than that of the partially optimized structure a, which remains the global minimum with an energy of −5.42 kcal/mol. This shows that the main features of the angular dependence can be captured with the rigid scan. The same is shown through the rigid angular scan for the Cl−Cl distance of the partially optimized minimum (2.105 Å) given in Figure S3. The main features are the same in both scans, and the energies are very similar; the main deviations observed in the maxima are most likely due to an increased repulsion. The diameters of the cage, measured as the distance between the opposing oxygen atoms to which the dihalogen points to, are 7.62 Å for the partially optimized structure a and 7.89 Å for partially optimized structures b and c. The corresponding optimized Cl−O distances are 2.76, 2.90, and 2.85 Å for a, b, and c, respectively. Thus, the dichlorine is stabilized in the minor axis of the 512 cage, and the observed Cl−O distances fall in the expected region for a strong halogen bond. One objective of this study is to analyze whether halogen bonds can be formed inside clathrate cages. In this respect, the optimized trimer (Cl2 − 2H2O) represents a reference structure to compare with since the only intermolecular interactions it exhibits are two halogen bonds and, as stated above, this structure is ubiquitous for the most stable stuctures in the cage.

known as the sigma hole. The negative electrostatic potential matches the area where the oxygen atoms are localized. The minimum energy paths between nearby minima of the PES (see Figure 2a) have the dihalogen moving along the hydrogen bonds that link the oxygen atoms associated with the minima (see Figure S1). The energy barriers are less than 2.0 kcal/mol indicating the possibility of rotation inside the cage through paths distributed almost spherically. In other words, the only forbidden regions for movement are along the center of the pentagonal faces. This is in agreement with the predicted angular distributions obtained in molecular dynamics simulations, which are nearly flat for this cage.21 It is interesting that the predicted mobility for dichlorine was not clearly seen in the X-ray studies1 where only a single structure was determined. This is in contrast to the case of bromine in the 51262 cage where delocalization is observed. The structures with Cl2 oriented toward the center of pentagonal faces are the less favorable; in those areas, the electrostatic potential is slightly positive and the charge density is very small. Although the electrostatic potential directs the halogen orientation to a great extent, there are features that cannot be explained solely in terms of the electrostatic potential; for example, the plane of the electrostatic potential for ϕ = 135° (Figure 2c) contains the global minima (a, θ = 45) and one local minima (d, θ = 120) of the PES. According to the electrostatic potential, the point d should be more stable than a due to the more negative electrostatic potential at which the Cl2 points to, but its inclusion energy is less by 0.51 kcal/mol. The relevance of other components of the interaction energy such as exchange repulsion, induction, and charge transfer has been documented on previous studies for isolated cages.16,36 In particular, exchange repulsion can be a determining factor for the most stable orientations. Although dispersion is the main attractive component, it is more isotropic and does not play a significant role in determining the most stable geometries.16,36 In order to observe the effect of the whole crystal in the Cl2cage interaction, the potential energy surface of the 512 cage extracted from the sI structure was calculated. The global minimum and maximum inclusion energies in that case are −6.11 and 0.93 kcal/mol, that is, both stronger attractive and repulsive interactions, respectively. This is most likely due to stronger electrostatic interactions in the isolated cage, which decreases upon saturating all oxygen lone pairs when forming D

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 4. NCI analysis of Cl2 @ 512: (a) isosurface of the partially optimized minimum and (b) rigid scan maximum (isovalue = 0.45). 2D plots of (c) trimer Cl2 − 2H2O, (d) partially optimized minimum, (e) rigid scan minimum, and (f) rigid scan maximum structures.

In the maxima, there are five oxygen atoms near each chlorine, but the geometrical parameters are not favorable for halogen bond interactions. The NCI index was calculated to search all possible noncovalent interactions. In Figure 4a,b, the isosurface for the RDG of the partially optimized minimum and rigid scan maximum are shown; the isovalue is 0.45. The noncovalent interactions present in the structure are the hydrogen bonds in blue and several others in green color between the halogen and water molecules. NCI isosurfaces of all structures look alike. In order to make a deeper analysis of the interactions, the 2D RDG versus sign(λ2)ρ plots are shown in Figure 4c−f. The interactions shown in Figure 4 can be classified according to their value of sign(λ2)ρ. The strongest appear for values less than −0.03, a set of weak interactions both attractive and repulsive appear for values less than |0.01|, and one peak appears between values −0.02 and −0.01. Inspecting the isosurfaces of each peak, one can deduce the kind of interaction it represents. The strongest interactions are

The similarity between the trimer and the partially optimized global minimum in their geometrical parameters is remarkable such that the dichlorine orientation in the optimized structure allows the formation of two halogen bonds. The only important difference is the distance of Cl−Cl; in the cage, it is shorter in comparison to the trimer. An elongation of the Cl−Cl covalent bond is observed in the presence of a halogen bond,15,37 correlating with the strength of the XB. The Cl−Cl distance of the isolated molecule calculated with the DFT methodology is 2.101 Å, suggesting a weaker halogen bond in the cage than in the trimer. The interaction energy of the trimer and the inclusion energy of the partially optimized minimum are very similar; however, in the cage, dichlorine interacts with the other water molecules, too, and part of the net attractive binding surely comes from this. The geometrical and energetic information for the partially optimized minimum is consistent with the presence of two halogen bonds, but further inspection is necessary to properly characterize it. E

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 5. Isosurface of the ELF in the (a) trimer Cl2 − 2H2O and (b) partially optimized minimum (isovalue = 0.85). (c) ELF along the path Cl− O for the structures analyzed of the 512 cage.

maintain the cage structure, it is natural to expect that they are not available to interact with the halogen molecule. As stated in the Introduction, the nature of the dichlorine-cage interactions is more subtle, and strong evidence for halogen bonding has been given in the case of the dodecahedral cage.16,20 In order to observe the electron localization in the crystal, in particular the oxygen lone pairs, the computation of the ELF was made. Values of the ELF above 0.8 are recommended to observe the covalent bonds and the lone pair electrons.33 The ELF isosurfaces with a 0.85 isovalue for the trimer and partially optimized minimum are shown in Figure 5a,b. In both figures, the lone pairs of oxygen atoms are observed. Also, notice that the halogen lone pairs are localized in a torus around the chlorine atoms, and the hole of the torus at the ends of the Cl−Cl bond corresponds to the sigma hole. The expected lock and key structure characteristic of halogen bonding can be clearly seen in the small cluster and more importantly in the clathrate crystal, confirming previous analysis performed on single cages.20 A more quantitative analysis can be performed by computing the ELF along the path that links the oxygen and chlorine atoms, thus proving that the lone pairs indeed point toward the sigma hole. The 2D plot of the ELF is shown in Figure 5c. Only one Cl−O path in each case is shown. In the x axis, extremes are the Cl and O atomic centers and all distances are normalized. The plot shows the shell structure of the atoms; from left to right, the peaks belong to Cl 1s, Cl 2s2p, Cl 3s3p, O 2s2p, and O 1s. For the electrostatic interaction between O and Cl to be optimal, it is necesary that the sigma hole (electrophilic region) is aligned with the oxygen lone pair. This is reflected when the path linking chlorine and oxygen has lower ELF values in the Cl 3s3p peak due to the sigma hole and high ELF values in the O 2s2p peaks due to the electron pair. The O 2s2p peaks in almost all cases analyzed reach values above 0.8, which means that the lone pair can establish a directional interaction with the Cl atoms. In the case of the geometry of the maxima, the O 2s2p is below 0.8, so the oxygen lone pair is not in the halogen bonding path, as expected. In the case of the chlorine 3s3p peaks, there are two kinds of them: the first ones with maximum values of 0.72 are paths through the sigma hole and appear for the partially optimized minimum, rigid scan minima, and trimer systems, and these cases are consistent with the electrostatic model of the halogen bond. In the second kind, the peak for the rigid scan maximum has a value of 0.9, belonging to a path across the torus of the Cl2 lone pairs, which presumably leads to a repulsive interaction between the oxygen and Cl2 lone pair charge density.

hydrogen bonds. The weakest attractive interactions are between the halogen and the distant oxygen atoms; due to their low density values, these interactions must be dispersion. The repulsive interactions are steric repulsion between the halogen and the hydrogen bonds that hold the cage structure. All previous interactions appear in the NCI analysis of all structures. The isosurface for the peak with density values of sign(λ2)ρ from −0.02 to −0.01 represents interactions between Cl atoms and the closest oxygen atoms; comparing this peak with the one found in the trimer Cl2 − 2H2O (Figure 4c), it strongly suggests the presence of a halogen bond. There are patterns that help differentiate the cases shown. For example, the peaks shift to lower values of sign(λ2)ρ following the sequence trimer, partially optimized minimum, rigid scan global minimum, and rigid scan global maximum with sign(λ2)ρ values of −0.018, −0.018, −0.016, and −0.013, respectively, indicating the strongest noncovalent interaction for the trimer and the weakest for the maximun. It is remarkable that the peaks of the trimer and the partially optimized minimum appear at the same value of sign(λ2)ρ suggesting the same interaction strength, a point that was already suggested on the basis of geometrical and energetic characteristics. Nonetheless, the optimum O−O distance in the trimer is 7.68 Å, and the equivalent, but restricted, O−O distance in the cage is smaller: 7.62 Å. This is consistent with the fact that the Cl−Cl distance in the cage is also slightly shorter. Additionally, for the case of the maximum, the peak shifts to the region of dispersive interactions where the orientation is not favorable to form halogen bonds, but due to the closeness of several oxygen molecules, the total dispersive interaction becomes stronger. Thus, the NCI analysis is very useful for identifying strong halogen bond interaction, but if it becomes weakened, then care must be taken to differentiate it from the case of mainly dispersive interactions such as those present for the maximum. In this respect, it is worth mentioning that the ″multidirectional halogen bonding″1 would be better described as a sampling of geometries that encompass a range of interactions, from halogen bonds with the highest strength to mainly dispersive attraction with lower strength. In addition to the strength of the interactions, the analysis relies on other characteristic features such as geometry and electrostatic component of the interaction. A halogen bond occurs when there is evidence of a net attractive interaction between an electrophilic region associated with a halogen atom in a molecular entity and a nucleophilic region in another, or the same, molecular entity.38 In the case of clathrate hydrates, the nucleophilic region would be associated with the oxygen lone pair electrons; however, since these are engaged in forming two hydrogen bonds to F

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 6. Charge density difference isosurfaces for the (a) partially optimized minimum, (b) rigid scan minimum, (c) rigid scan maximum, and (d− f) dimers/trimer extracted from the partially optimized minimun of Cl2@512 cage (isovalue 0.0006 au).

Figure 7. (a) Angular dependence of the PES for Cl2@51262, (b) electrostatic potential mapped into a charge density isosurface (isovalue 0.006 au), (c) sample of structures with ΔI less than 1 kcal/mol with respect to the global minimum, and (d) molecular structure of the 51262 cage in the crystal sI.

Figure 6a−c shows the electron density difference isosurfaces in the partially optimized, rigid scan minima, and maxima. The blue color designates an increase, and the red is a decrease of charge density. The Hirshfeld charges of the Cl2 are shown in Table 2. The dihalogen charge in the cage is slightly positive, which is opposite to those reported with NBO21 and with QTAIM,20 albeit in these cases the negative charges are fairly small and calculations were performed on single cages. All results are consistent with the fact that charge transfer is small for halogen bonding in these systems. The different predictions may reflect the different methods used to evaluate the charges but could also be related to differences between single cage models and the full crystal. Nevertheless, we will show that all studies agree in overall trends. The charge in the dihalogen inside the 512

In the region between the atoms, the values of the ELF are small, but the subtle differences between the structures are revealing. The partially optimized minimum, rigid scan minima, and trimer systems have values greater than 0.045 ,and the rigid scan maximum has a value of 0.035. Even though the differences are small, the electron localization in that region characterizes the noncovalent interactions and is consistent with what was shown in the NCI analysis. In addition to the electrostatic component, the charge transfer and the polarization are other components of the halogen bond.39 In order to analyze these components, we have computed the rearrangement of charge density of both guest and host in the crystal through the Hirshfeld-I40 charges and the electron density differences. G

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

isosurface. A correlation between the electrostatic potential and the minima and maxima in the PES is shown in Figure 7a,b: the letters are for the minima, and the numbers are for the maxima. The labels in the electrostatic potential indicate the regions where the dichlorine points to. The corresponding inclusion energies are shown in Table 3.

cage decreases when the attractive interaction energy increases. To understand the rearrangement of the charge density during the inclusion of the dihalogen in more detail, the charge density differences were calculated for three systems extracted from the partially optimized minima (Figure 6d−f): two trimers and one pentamer Cl2 − 4H2O system. The water molecules considered were those nearest to the Cl atoms (first set, Figure 6d) and others linked to them by a hydrogen bond (second set, Figure 6e). The charge density differences were calculated from the density of the total cluster and those of the water cluster and the chlorine molecule. In Figure 6d, one can appreciate that the first set of water molecules donates charge density to the dihalogen molecule as expected for the halogen bond, Figure 6e shows that the second set of water molecules is capable of receiving charge density from the dihalogen, and Figure 6f shows both phenomena. The Hirshfeld charges calculated for the dichlorine in these systems are −0.07, +0.04, and −0.06 for panels (d−f), respectively. Thus, in the cage, the dichlorine donates charge density to the distant water molecules and receives charge density from the nearest water molecules through a halogen bond. The stronger the halogen bond, the more charge density is received and the lower the chlorine charge. This same trend can be observed in the cages (Figure 6a−c); the charge density loss of the nearest water molecules decreased when the attractive interaction energy diminishes. On the other hand, Figure 6 shows regions of charge density gain between the halogen atoms and the nearest water molecules in the trimer, partially optimized minimum, and rigid scan minima structures. This indicates a covalent part of the interaction as expected for the halogen bond.41 The rigid scan maxima does not show this type of region. The nearest water molecules suffer a charge density loss in the lone pair regions, and thus, the hydrogen bonds, where that water molecule acts like a hydrogen acceptor, are weakened. In the case of the partially optimized structure, the hydrogen bond distance increases from 1.760 to 1.764 Å when the dichlorine is inserted in the cage; this weakening is also reflected in the NCI 2D diagram: the peak below −0.04 sign(λ2)ρ belonging to that interaction is present in the rigid scan minimum but disappears from the partially optimized structure. The other hydrogen bond, which is in one edge of the cage, was weakened, too, and the O−H distance increases from 1.784 to 1.790 Å. 51262 Cage. The PES of the 51262 cage, calculated through the rigid angular scan (see Methodology), is shown in Figure 7. All orientations of the dichlorine have an attractive inclusion energy larger than in the 512 cage. This stabilization explains why all the 51262 cages are occupied as has been observed experimentally.1 One can clearly see that the region defined along θ = 90 is preferred by the Cl2 molecule and the less attractive energies are along θ = 0. These last orientations have the chlorine molecule pointing toward the center of the hexagonal faces of the cage. The average diameters of the cage are 9.0 and 6.0 Å along θ = 90 and θ = 0, respectively. The van der Waals diameter of the chlorine molecule is 5.5 Å indicating the expected repulsion for the θ = 0 orientations. The inclusion energies of the global minima and maxima are −6.78 and −4.67 kcal/mol. In order to analyze the electrostatic component of the dichlorine−cage interaction, the electrostatic potential of the empty cage was calculated and mapped in an electronic density

Table 3. Inclusion Energies (I) in kcal/mol for the Stationary Points Shown in Figure 7a,b point I

a −6.78

b −6.73

c −6.65

d −6.26

1 −4.67

In the global minimum, the chlorine atoms point to negative areas of the electrostatic potential, which indicate that the electrostatic interaction is the greatest for this dihalogen orientation; however, in the range of θ between 60 and 120° where the inclusion energy differences are less than 1 kcal/mol, the dichlorine points to both slightly negative or positive areas. Still, it should be noted that the dihalogen always points to the regions where the electrostatic potential attains its lowest values and thus plays a role in determining the most stable geometries. The small interaction energy differences over a wide region of angular space indicate the possibility of movement of the dihalogen (Figure 7c). This is in agreement with the angular distribution observed for this system with molecular dynamics21 where the angular distribution function shows a maximum around θ = 90. This also coincides with the highest probability positions of the dichlorine determined by X-ray measurements.1 In order to analyze some representative cases, four structures were chosen (Figure 8): the minima and maxima of the PES, the restricted optimization of the minima (see the Partial Optimizations section), and a structure with unfavorable geometric parameters for the halogen bond but with an inclusion energy similar to all of the minima (ϕ = 135, θ = 90), which we will name X. The last structure is representative for several cases contained in Figure 7c. It is important to clarify that all these structures do not correspond to stationary points. The inclusion energy and the geometrical parameters of the chosen structures are shown in Table 4. Since the b orientation (Figure 7) has an inclusion energy very close to the global minimum, it was partially optimized, too. The resulting interaction energy is −6.92 kcal/mol, and the distances Cl−O are 3.44 and 3.41 Å. In this case, because the halogen is along the 9 Å axis, both Cl atoms cannot get close enough to the opposite oxygen atoms. The inclusion energy for the partially optimized global minima is −7.57 kcal/ mol, only 0.7 kcal/mol below the partially optimized b geometry even though this last one does not satisfy the geometrical criteria that favor the formation of a halogen bond. Comparing the partially optimized global minima of the 51262 and 512 cages, we notice that, for the former, the Cl−O distance is larger and the Cl−Cl−O angle deviates more from colinearity; this indicates that the XB is weaker in this cage. Another important difference is that the two closest Cl−O distances are different indicating that one of the XB interactions is stronger than the other. This was not the case for the smaller cage where the distances are equivalent. This result agrees with the X-ray measurements1 where it is observed that the halogen molecule is not located in the center of the cage. H

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 8. Geometries of the analyzed structures for the 51262 cage: (a) partially optimized minimum, (b) rigid scan minimum, (c) X, and (d) rigid scan maximum.

Table 4. Inclusion Energies, Geometrical Parameters, and Hirshfeld Charges in Cl2 for the Analyzed Structures of Cl2@51262 Cagea structure

I

Cl−O

Cl−Cl−O

Cl−Cl

Cl2 Hirshfeld charge

partially optimized minima rigid scan minima X rigid scan maxima

−7.57 −6.78 −6.66 −4.67

2.97, 3.01 3.05 3.51 3.60 (2O), 3.35 (4O)

171.38, 163.39 167.02 156.98 123.66, 126.45

2.112 2.00 2.00 2.00

0.07 0.08 0.11 0.13

a

The energies are in kcal/mol, the distances are in Å, and the angles are in degrees.

Figure 9. NCI analysis of Cl2 @ 51262. 2D plots of (a) partially optimized minimum, (b) rigid scan minimum, (c) X, and (d) rigid scan maximum structures.

time, the smaller size of the 512 cage should lead to a larger

It is important to clarify that there is no contradiction in the fact that XB interactions are stronger in the 512 cage than in the 51262 cage even though the total inclusion energies are larger for the latter. Other than electrostatic interactions, the main attractive component comes from dispersion, which is largely isotropic,16,36 and given the larger number of water molecules in the 51262 cage, it should lead to larger attraction. At the same

repulsion. So, as stated at the beginning of this section, we find that the inclusion energies are more attractive for the larger cage for all orientations. The best conditions for XB are found in the smaller cage, but these represent a partial contribution to the overall inclusion energy. I

DOI: 10.1021/acs.jpcc.9b03268 J. Phys. Chem. C XXXX, XXX, XXX−XXX

Article

The Journal of Physical Chemistry C

Figure 10. ELF isosurface of the (a) partially optimized minimum and (b) along the path Cl−O for the structures analyzed of the 51262 cage.

Figure 11. Isosurfaces of the ρ difference of the (a) partially optimized minimum, (b) rigid scan minimum, (c) X, and (d) rigid scan maximum structures of the 51262 cage.

The NCI index was calculated for these structures in order to further characterize the existing noncovalent interactions. The isosurfaces are similar to those of the 512 cage (Figure 4a). Figure 9 shows the 2D RDG versus sign(λ2)ρ plots. The peaks in the plots belong to hydrogen bonds in values of sign(λ2)ρ lower than −0.03, the dispersion and the steric repulsion interactions are between −0.01 and 0.01. In the plots for the partially optimized global minimum and rigid scan minima, there are peaks between the hydrogen bonds and the dispersion interactions; the isosurfaces of these peaks are between the chlorine and the oxygen atoms of the nearest water molecule, and this is consistent with a weakened XB as previously mentioned, in fact weaker than the interaction in the maxima of the 512 cage. On the other hand, the plots for the maxima and X structure do not show that peak; that is, there is no NCI evidence for the presence of the halogen bond in those structures, as expected. If a halogen bond is present, the lone pair electron of the water molecules must point to the sigma hole of the halogen. The ELF was calculated to establish if the lone pairs of the oxygen atoms are available to interact with the dichlorine in this cage, and it is shown in Figure 10. In general, the form of the 2D ELF plot is very similar to that of the 512 cage. The peaks of the states O 2s2p are above 0.80, and the lone pairs of oxygen are capable of establishing an electrostatic interaction with the dichlorine. On the other hand, there are three groups of 3s3p Cl peaks through the path Cl−O: the states of the partially optimized and minima with maximum values below 0.7, the states of the X structure with maximum values between 0.7 and 0.8, and the states of the maximum with maximum values of 0.9. This indicates that, for the most stable structures, the path points to the sigma hole region, as expected. In the other systems, the path passes through the lone pair electron region of Cl, which

probably induces electrostatic repulsion with the lone pair of oxygen. In the region between the atoms, there are subtle differences in the values of the ELF. Although these are small, there is a trend. The minimum values of the ELF are 0.035, 0.03, and