Molecular electrostatic potential as a guide to intermolecular interactions: challenge of nucleophilic interaction sites

Introduction

This year (2025) is the International Year of Quantum Science and Technology (IYQ) in recognition of the 100-year anniversary of the initial development of quantum mechanics. Already in 1929 Paul Dirac, one of the pioneers in quantum theory development, stated: The underlying physical laws necessary for the mathematical theory of a large part of physics and the whole of chemistry are thus completely known, and the difficulty is only that the application of these laws leads to equations much too complicated to be soluble. ¹

In 1998 the Nobel committee of chemistry, when awarding the Nobel Prize of Chemistry to Walter Kohn and John Pople, found that Dirac’s statement had been overly pessimistic and stated: Things began to move at the beginning of the 1960s when computers came into use for solving these equations and quantum chemistry (the application of quantum mechanics to chemical problems) emerged as a new branch of chemistry. As we approach the end of the 1990s we are seeing the result of an enormous theoretical and computational development, and the consequences are revolutionising the whole of chemistry. ²

This development has continued and today we are able to compute highly accurate wavefunctions and energies for most chemical systems. Instead it has become apparent that the main problem no longer lies in solving the complicated equations defined by quantum theory but rather to interpret the increasingly complex wavefunctions that result from the solving in terms relevant for the understanding of chemistry. It is illustrative that the two most commonly used quantum based explanatory models in chemistry still are the periodic table of the elements and the drawing of Lewis structures, ^{3 , 4} including the latter’s extensions to resonance structures and curly arrow pushing. ^{5 , 6} These models, which are simple to use but convey important chemical information, have solid foundations in quantum theory although being introduced before the advent of quantum mechanics.

In the field of non-covalent interactions, the surface electrostatic potential (V _S (r)) has emerged as a widely used tool for extracting information of complex quantum chemical calculations to understand and characterize interaction sites of molecules. It was first introduced in the 1990s in the research group of Peter Politzer, ⁷ and some of the first applications were biological recognition interactions and hydrogen bonding. ^{8 , 9} In 1992 Brinck et al. demonstrated that the capacity of halocarbons to interact with Lewis bases can be explained by a tip of positive V _S (r) at the end of the halogen atom opposite to the carbon halogen bond. ¹⁰ This finding was ignored until the boom of halogen bonding in the beginning of the 21st century when it was realized that the surface electrostatic potential can be an effective tool for characterizing halogen bond donating sites. ^{11 , 12 , 13 , 14 , 15 , 16 , 17 , 18 , 19 , 20}

The relevance of V _S (r) was cemented in 2007 when Clark et al. provided an orbital-based rationale for the electron deficiency at the tip of the halogen atom and the resulting maximum in V _S (r) (V _S,max ); according to this rationale the positive tip can be called a σ-hole because of the polarization of the σ-orbital density away from this region towards the C–X bond. ²¹ The σ-hole concept was later expanded to other main-group elements, and has given rise to a wide range of new interactions named after the corresponding group of the periodic table, e.g. chalcogen, tetrel, and pnictogen bonding. ^{18 , 22 , 23 , 24 , 25 , 26} The term π-hole was introduced to describe surface electrostatic potential maxima at positions that are located perpendicular to a planar molecular framework. ^{25 , 27} More recently the concepts of σ-hole and π-hole bonding and the relevance of positive V _S (r) regions for interactions with Lewis bases has been extended to transition metal-based compounds, as pioneered by Brinck and coworkers. ^{28 , 29 , 30 , 31 , 32 , 33}

Even though the surface electrostatic potential (V _S (r)) originally was developed to characterize sites that interact with electrophiles and Lewis acids as well as sites that interact with nucleophiles and Lewis bases, it is clearly the latter application that has been dominating. This is apparent from the extended use of V _S (r) maxima (V _S,max ) as indicators of σ-holes and π-holes. The minima in V _S (r) (V _S,min) should in a similar manner be indicative of sites that interact with Lewis acidic sites, e.g. σ-holes and π-holes, but such examples are less common in the literature. The reason is probably that Lewis basic (nucleophilic) sites are more challenging to characterize and their location depends on the Lewis acid/electrophile. As an example, acetone has two equal minima in the surface electrostatic potential (V _S,min ) close to the top of the O, at angle of around 160° with the C–O bond, while hydrogen bond donors and most other σ- and π-hole donors bind at an angle of 120–130°. This angle is in better agreement with the angle of the spatial minima of the electrostatic potential (V _min ).

Similar observations have also been made for the interactions of weak halogen bond donors, i.e. dihalogens, in that the angular direction of interactions with Lewis bases is better described by V _min than V _S,min . ^{34 , 35} However, metal cations, e.g. Li⁺, bind acetone at an angle of nearly 180°, in relatively good agreement with the angle of the V _S,min . ³⁵ In a similar way as for the prediction of the positions of hydrogen bond accepting sites, the V _min has also been shown to more accurate than the V _S,min for predicting the strength of hydrogen bonding; it has been shown for families of hydrogen bond acceptors that there are better correlations between V _min and the hydrogen bond strength than between V _S,min and hydrogen bond strength. ^{9 , 36}

In this work we have analyzed a number of intermolecular interactions between electrophiles (Lewis acids) and nucleophiles (Lewis bases) and attempted to understand how and why the location of the site for the interaction with the electrophile is dependent on the character of the electrophile.

Methodological background

The electrostatic potential (V(r) that surrounds a molecule is rigorously defined by.

where Z _A is the charge of nucleus A located at R _A. V(r) is a physical observable that can be determined from an experimental or theoretical density distribution function. It is a first order property and thus it is relatively insensitive to the computational method. qV(r) gives the interaction energy between the unperturbed charge distribution of the molecule and a point charge q located at r.

The electrostatic potential of an atom is everywhere positive and decrease asymptotically towards zero with increasing distance from the nucleus. When a molecule is formed from atoms the electron density is redistributed leading to regions of negative V(r) with the minima (V _min ) typically located above π-bonds or at lone pair regions. As an example, V(r) of acetone depicted in the plane of the heavy atoms is shown in Fig. 1. There are two equal V _min at positions corresponding to the lone pairs. These also reflect sites susceptible toward interaction with electrophiles or Lewis acids. ^{9 , 37} There are no corresponding maxima in V(r) as the most positive regions are found at the nuclei where V(r) goes toward infinity. The nuclear positions are not sites for interaction with electrophiles, as the electron density of the approaching electrophile will be repelled by the electron density of the nucleophile before reaching the nucleus.

Fig. 1:

On the left, V(r) in the molecular plane of acetone with V _min labeled by _*. On the right, V _S (r) on the 0.001 au isodensity surface with V _S,min labeled by _*. The V _min positions, which are located inside of the 0.001 au surface, are here indicated blue envelopes that encompasses the V _min . Similar coloring schemes are used for both V(r) and V _S (r) with the most negative areas being deep blue whereas the most positive are read. Note that 1 kcal/mol is equal to 4.184 kJ/mol.

In 1989, Sjoberg and Politzer recognized that the conundrum of identifying electrophilic sites from the electrostatic potential could be solved by computing and plotting V(r) on a molecular surface defined by constant electron density, resulting in a surface electrostatic potential (V _S (r)). ⁷ Typically the 0.001 or 0.002 a.u. contours are used and the corresponding surfaces are at distances from the nuclei that are close to the van der Waals radii of the corresponding elements. ³⁵ Such a surface is expected to have a nearly constant repulsive potential, and therefore the maxima in V _S (r) (V _S,max ) can be used to locate the sites most susceptible to interact with nucleophiles in interactions dominated by electrostatics. ³² This approach has proven highly successful as evidenced by the general adoption of V _S,max as indicators of σ-hole and π-hole donating sites.

Sjoberg and Politzer intended the surface minima (V _S,min ) to serve a similar purpose for identifying sites susceptible to interact with electrophiles, or in today’s vocabular to identify σ-hole accepting sites. The left part of Fig. 1 shows the V _S (r) at the 0.001 density contour for acetone. Similarly to the V _min , there are two V _S,min outside the oxygen in the same plane as the heavy atoms. However, each V _S,min is located at an angle of 162° with the C-O bond, and the V _S (r) value at the top of O (at 180° angle) is less than 0.1 kcal/mol higher, showing the small variation of V _S (r) close to top of the O. However as already mentioned, this location is not representative of the interaction site of many electrophiles, such as hydrogen bond and halogen bond donors, but is close to the interaction site of metal cations. In the figure of V _S (r) we also show the locations of the V _min , which are found well inside the volume of the 0.001 contour at an angle of 128° with the C–O bond.

Results and discussion

Interactions of acetone

We will begin by analyzing complexes between acetone and few electrophiles, i.e. Lewis acids and/or σ-hole donors, and try to relate their geometries to the spatial distribution of the electrostatic potential as well as the surface electrostatic potential (Fig. 2). Beginning with the HF complex, we note that the hydrogen of HF interacts with the lone pair of O of acetone at a C–O–H angle of 118°, which is slightly smaller than the C–O–V _min angle (128°) and much smaller than the C–O–V _S,min angle of 162°. Furthermore, the F atom is slightly tilted towards the CH₃ group of acetone at an O–H–F angle of 170°. It can also be noted that the H atom is located inside of the van der Waals volume as defined by the 0.001 au contour, but it is further away from the O atom than the V _min . Thus, it is not obvious why the C–O–H angle is much closer to the C–O–V _min angle than it is to the C–O–V _S,min  angle.

Turning towards acetone-BF₃ complex, we find that the C–O–B angle (128°) is in perfect agreement with the C–O–V _min (128°), and thus much smaller than C–O–V _S,min angle (162°). The short O–B distance of 1.65 Å is indicative of a dative covalent bond and the B is located inside the 0.001 au contour but the bond is significantly longer than the O–V _min distance.

Fig. 2:

Optimized structures of complexes between acetone and different electrophiles. All neutral complexes interact with similar angular directions as the C–O–V _min  angle.

The two halogen bonded complexes acetone-BrF and acetone-Br₂ are very similar in their angular orientation. The C–O–Br angles are 128° and 126°, respectively, in both cases in close agreement with the C–O–V _min (128°). The O–Br–X angle is for both almost 180°, in agreement with earlier observations that the halogen bond is more directional in character than the halogen bond (compare with the corresponding 170° angle in the acetone-HF complex). It is noteworthy that the C–O–Br angle is close to the C–O–V _min angle in both complexes, even though the O–Br distance is much longer in the Ace-Br₂ complex_. This confirms the previous observation that there is no correlation between how close the intramolecular angle, e.g. C–O–Br, is to the angle with V _min and how close the σ-hole is to the V _min position.

The Li⁺ forms a complex that is linear with the C–O bond, which means that the C–O–Li angle of 180° is larger than the C–O–V _S,min angle of 162°. However, as previously discussed the V _S (r) has a flat value distribution close to the top of O, and the top value at 180° deviates less than 0.1 kcal/mol from V _S,min value. Kollman and Rothenberg already in 1977 identified that Li⁺ forms a nearly linear complex with many Lewis bases whereas the corresponding interaction with hydrogen bond donors is at an angle that agrees with the V _min . ³⁸ They argued that the reason is that exchange repulsion is important for determining the geometry of the Li⁺ complex, whereas its importance is smaller for the hydrogen bonded complexes, as the electron of hydrogen is polarized toward the electronegative atom. This could explain the better agreement with V _S,min for the Li⁺ complex since V _S (r) indirectly takes exchange repulsion into account by using a surface of constant electron density. However, if this rationale is correct then LiF should also form a complex with acetone with a linear C–O–Li bond/interaction, as there is no reason for the exchange repulsion to be of lower magnitude for complexes of LiF compared to Li⁺. Instead we find a C–O–Li angle of 126° in close agreement with the V _min angle and the angles of the other complexes with neutral nucleophiles. An indication of why monatomic cations behave differently from neutral molecules is given by the position of F in the complex of acetone with LiF; the negative F is strongly tilted towards the region of positive V _S (r) above the CH₃ group. Thus, LiF is positioned to maximize the total electrostatic interaction with acetone, as represented by V _S (r), and not only the interaction with the positive end, i.e. Li. In the case of a monoatomic cation, the V _S,min position should be the most favorable position on the surface as V _S (r) gives the electrostatic interaction energy with a unitary point charge. Neutral molecules, however, can often be approximated by a dipole or a collection of dipoles, and the interaction energy with a dipole is proportional to the electric field in the direction of the dipole. This can explain why the interaction angle with a σ-hole donor often coincides with the angle of the V _min ; the electric field perpendicular to the surface is likely to be strongest outside of the V _min .

Interactions of thioacetone

We have also analyzed how the directional binding preferences of thioacetone relates to the spatial and surface distributions of its electrostatic potential (Fig. 3). Looking first at V _S (r) in a plane through the heavy atoms, we note that the overall picture is similar to that of acetone. However, the two lone pair V _min are located at a smaller angle, 107° instead of 128°, and are less negative, −31.5 kcal/mol instead of −51.9 kcal/mol, in comparison to acetone. These differences are translated into V _S (r) with the V _S,min of thioacetone appearing at a much smaller angle compared to that of acetone, 114° instead of 162° and with value of −23.6 instead of −37.5 kcal/mol. In addition, the flat V _S (r) of acetone at the angles close to 180° has no equivalence for thioacetone and V _S (r) at 180° is as high as −6.8 kcal/mol.

Fig. 3:

Corresponding to Fig. 1 but for thioacetone instead of acetone. Similar coloring schemes have been used as in Fig. 1. Note that thioacetone has less negative V _min and V _S,min than those of acetone and are positioned at smaller angles relative the C–S/O bond.

In general, the binding directions of the different electrophiles relates to V _min and V _S,min in a similar manner as for acetone (Fig. 4). The HF molecule binds with a C–S–H angle (100°) that is slightly smaller than the C–S–V _min angle and with the F slightly tilted towards the CH₃ group in a similar manner as for acetone, i.e. he S–H–F angle is 168°. BF₃ binds at an angle that nearly coincides with the C–S–V _min angle, and that also holds for Br₂. The Br₂ is again found to interact with nearly linear conformation with respect to electrophilic site; the S–Br–Br angle is close to 180°, confirming the more localized σ-hole of halogen bond donors compared to the σ-hole of hydrogen bond donors, such as HF.

Fig. 4:

Optimized structures of complexes between thioacetone and different electrophiles. All neutral complexes interact with similar angular directions as the C–O–V _min angle. Li⁺ interacts at an angle closer to the C–O–V _S,min  angle.

In the case of Li⁺, the position is again found to be determined by V _S (r) with the C–S–Li angle (112°) in good agreement with the C–S–V _S,min angle (114°). The interaction of LiF is again dramatically different from that of Li⁺ with LiF having a C–S–Li angle of 101°, which is almost equal to that of HF and smaller than the C–S–V _min angle. In a similar manner as for the acetone•••LiF complex, the F is strongly tilted towards the positive V(r) above the CH₃ and the O–Li–F angle is 129°. The corresponding angle on acetone•••LiF is 127°. These results show that the σ-hole on LiF is even more delocalized than the σ-holes of hydrogen bond donors. Furthermore, the strong tilt of LiF is in agreement with our conclusion that the geometry of the complex is optimized to maximize the electrostatic interaction between the entire charge distributions of the two interacting molecules; thus the negative F of LiF plays an important role in the interaction.

Interactions of hydrogen halides

To investigate if the observations that were made regarding the different angular preferences of cationic and neutral electrophiles are generally valid, we have also investigated some complexes where hydrogen fluoride and hydrogen chloride serve as nucleophiles or Lewis bases (Fig. 5). We will first consider V(r) and V _S (r) of HF (Fig. 6). V(r) in a plane through the heavy atoms has two minima of −28.8 at 124° (H–F–V _min ), but because of rotationally symmetry there is actually a continuous circle of minimum V(r) at 124° and a distance of 1.31 Å from F. V _S (r) also has a continuous circle of minima but at 140° and at a value of −23.8 kcal/mol. Notably there is also a V _S,max at the top of F, in other words a σ-hole. However, the V _S,max at −22.7 kcal/mol is only 1.1 kcal/mol higher than the V _S,min .

Fig. 5:

V(r) (left) and V _S (r) (right) for HF (top) and HCl (bottom). The depiction is similar to Fig. 1 and Fig. 3 but the coloring scheme is shifted as the magnitude of the negative potential range is smaller than for acetone. Note also that the V _min and V _S,min each form a circle of minimum potential at constant angle with the H–X bond.

Fig. 6:

Optimized structures of complexes between HF and HCl as Lewis bases and ClF and Li⁺.

HF forms a weak halogen bonded complex with ClF that has been estimated from microwave spectroscopy to have an H–F–Cl angle of 125° (±5°) angle and an F–Cl distance of 2.76 Å. ³⁹ Our computed structure has an angle of 121° and a distance of 2.62 Å. Considering that the experimental structure includes the lengthening due to anharmonicity of the weak intermolecular Cl–F interaction, the agreement is good. Notably a previous computation at the MP2/6–311++G (2d,p) level resulted in an angle of 125° and a distance of 2.75 Å in almost perfect agreement with experiment. ³⁴ The experimental estimate as well as the M06-2X and MP2 predictions for the angle are in very good agreement with the V _min angle, confirming that V _min accurately predicts the angle of interactions with neutral electrophiles.

The lithium cation forms a strong complex with HF at a Li–FH angle of 180°. Thus, this can technically be viewed as an interaction between two σ-holes, one on each reactant. However, as already mentioned, the V _S,max at the top of F is negative and only 1.1 kcal/mol higher than the V _S,min at 140°. The polarizability is expected to be significantly higher along the H–F bond than perpendicular to it, and thus the preference for a linear complex is most likely due to a larger polarization contribution to the intermolecular interaction at the 180° angle; this compensates for the slightly weaker electrostatic interaction. We also investigated the interaction between HF and LiF, but it results in the reorganization of the structure into a complex between Li⁺ and the symmetric (FHF)^-.

We finally turn towards the HCl molecule. When comparing V(r) and V _S (r) of HCl with those of HF, we see similar changes in the locations of the minima as when going from acetone to thioacetone. The V _min of HCl is at a smaller angle (106°) than that of HF and the V _S,min is at an angle (108°) that is much closer to the angle of V _min than is the case for HF. There is also a much larger difference between the V _S,top and V _S,min value for HCl than for HF. HCl, in contrast to HF, has a well-developed σ-hole at the end position with a positive V _S,max of 5.1 kcal/mol. ClF forms a weak complex with a similar but slightly shorter angle (101°) than the H–F–V _min angle. Li+, however, forms a strong complex with an angle (108°) that similarly to thioacetone coincides with the V _S,min  angle.

Conclusions

Neutral σ-hole and π-hole donors (electrophiles) generally interact with a Lewis basic (nucleophilic) site at an angle that closely corresponds to the angle of the V _min . We have shown that this holds for the full spectrum of polar interactions ranging from weak halogen bond interactions via hydrogen bonds to strong dative covalent interactions. The angle of the V _min is always smaller than the corresponding angle of the V _S,min and it is close to the angular direction where one intuitively would expect to find a lone pair. The rationale for the better agreement with the angle of V _min than V _S,min is that the charge distribution of the Lewis base, as represented by its electrostatic potential, interacts electrostatically with the entire charge distribution of the electrophile and not only with its electrophilic site. The charge distribution of most neutral electrophilic molecules can be approximated by a dipole or a superposition of dipoles, and a dipole interacting with the Lewis base has its lowest electrostatic energy when it is pointed along the direction of strongest electrostatic field, which generally coincides with the direction of the V _min . Further support for the rationale of strongest total electrostatic interaction can be seen in the intermolecular interactions of molecules such as LiF that has a wide and delocalized σ-hole on the electrophilic site, e.g. Li. They form complexes where the negative end (e.g. F) of the molecule is tilted toward a region of positive potential outside the other molecule, as such a geometry strengthens the favorable electrostatic interaction between the two molecules. On the other hand, monoatomic cations, like Li⁺, generally interact with Lewis bases in a direction that is aligned with the V _S,min . This is not surprising as the V _S,min corresponds to the position on the surface that has the strongest electrostatic interaction with a positive point charge. In some cases where two V _S,min are close in space, i.e. at angles close to 180°, it can be more favorable for the electrophilic site to interact at angle of 180° as this leads to a stronger polarization contribution to the interaction. Overall this study confirms the importance of electrostatics for determining the geometry of intermolecular interactions, even for interactions where other energy terms provide the main stabilizing effect, such as dative covalent bonds and some charge shift bonds.

Computational methods

The geometrical structures of molecules and molecular complexes have been computed at the M06-2X/6-31+G (d,p) level using the Gaussian 16 suite of programs. ⁴⁰ The M06-2X functional is highly accurate for main-group chemistry, including noncovalent interactions, and it has been parameterized to include short and medium range London dispersion. Di Labio et al. evaluated non-counterpoise-corrected DFT interaction energies against a revision of the HB23 data set, and found the M06-2X functional to perform well with a mean absolute deviation of only 0.21 kcal mol^{−1 41} A polarized valence-double zeta basis set augmented with diffuse functions is generally sufficient for computing geometries of intermolecular complexes when used together with a hybrid DFT functional, and we found minimal changes in geometrical parameters going from the 6-31+G (d,p) basis set to the much larger jun-cc-pVTZ basis set.

Electron densities for the analysis of V(r) and V _S (r) were also computed at the M06-2X/6-31+G (d,p) level of theory. V(r) is a first order property and as such less sensitive to the computational level than geometries or energies. As an example, V(r) for several of the molecules of this study had previously been analyzed using Hartree–Fock and a smaller basis set, but the angular directions of the V _min are in almost perfect agreement between the studies. ³⁴ V(r) was computed using Gaussian16 and V _S (r) using hs95 (ver.230419) of Brinck. ³² Visualization of V(r) and V _S (r) was done using Vesta and Gaussview, respectively. ^{42 , 43}