Julia Contreras-García

Roberto A Boto

Fernando Izquierdo-Ruiz

Igor Reva

Tatiana Woller

Mercedes Alonso

A benchmark for the Non-Covalent Interaction index (NCI) or...is it really all in the geometry?

Keywords:

Describing non-covalent interactions has shown to be of paramount importance in many areas of theoretical chemistry and related disciplines, such as biochemistry and material science. However, non covalent interactions are subtle effects, very difficult to reproduce from most common computational approaches. Electron density studies have shown to provide a good semi-quantitative visual approach to such interactions, which are much less prone to method dependency. But to which extent? This is the question addressed in this contribution. The reduced density gradient approach is given the third degree so as to provide the user with a benchmark on how it is affected by the computational method and the basis set of choice. We have assesed the dependence of the NCI results on the geometry. This last question is addressed in detail to dissect how, why and when the NCI method can be used to understand dispersion interactions. Along various examples we will show that the NCI index is very little dependent on the method and basis set used in the calculation of the electron density as long as the geometry is kept fixed. Indeed, the biggest variations in NCI come from changes in the geometry. Thus, methods which provide descriptions of a given interaction type of different accuracies, will yield different electron density organizations. This gives no qualitative variations in the NCI 3D picture. But it is reflected in quantitative NCI measures even in very subtle cases. Moreover, in the case of a failure of the calculation method, NCI can also reveal the sources of its error. NCI volumes are able to locate the energetic ordering in various conformational situations, but always in a relative manner. Absolute values should not be used in comparisons, nor between compounds that do not belong to the same family.

I. INTRODUCTION

Non-covalent interactions (NCIs) are of paramount importance in chemistry, biology and material science; 1,2 playing for example a key role in the interaction of a protein and a drug, as well as in the formation of molecular clusters. These classes of interactions span a wide range of binding energies, and traditionally encompass hydrogen bonding, dipole-dipole interaction and London dispersion. More specifically, NCIs set up the force field scenario through which chemical species interact with each other without a significant electron sharing between them. They represent the machinery through which molecules interact at long distance and establish how they will approach and pack together.

It is thus not surprising that during the last decade NCIs have been subject of a vast number studies from both theoretical and practical considerations. In general, many efforts are currently invested to reconcile the traditional chemical pictures with the outcome of first principles quantum-mechanical methodologies. 3,4 Quantitative and rigorous formalisms based on the topological and visual analysis of scalar fields related to the electron density, ρ( r), have shown to provide very useful information. Within this framework, some of us have developed an index for the analysis of non-covalent interactions derived from the electron density and its gradient. 5 Although mainly applied to non-covalent interactions, it can also be used to detect covalent interactions. 6,7 What is very interesting about its application to non-covalent interactions is the fact that promolecular approaches (sum of non-relaxed atomic densities) becomes a reasonable approach for the electron density and its derivatives. This means that NCI can be calculated in a very fast (and extremely parallellizable) manner for big systems. 5 This raises two theoretical issues:

• What is the influence of the method and the basis set on the NCI index? How trustworthy are promolecular results? In other words, a benchmark on the method dependency of NCI was conspicuously missing.

• Where does the energetic variability come in? Some weak interactions, such as dispersion, are due to correlation (and thus not contained in a promolecular approximation, nor in many DFT calculations). Thus, where does energetics play the role that leads to the stabilization of non-covalent interactions? Very simple approaches based on the try alone, so 3D space seems a reasonable approach to characterize these interactions too. How and when NCI derived from non-dispersive calculations can be a reasonable approach? and what is its interpretation?

We will first introduce the theoretical background and the computational details in sections II and III, respectively. In sections sections IV and V we will try to address the two issues presented above. The promolecular approach and the influence of computational methods (DFT and wavefunction methods) and basis sets when visualizing and quantifying NCI results will be tackled in section IV, and the uttermost importance of geometry in NCI calculations, but also its ability to disentagle subtle situations and methodological failures.

The contribution ends with a summary of the conclusions we have drawn on the use and limitations of NCI.

II. THEORETICAL BACKGROUND A. Reduced density gradient

The origin of NCI can be traced back to the exchange contribution of the generalized gradient approximation (E GGA X ) with respect to that of the local density approximation (E LDA X ):

E GGA X -E LDA X = - F (s)ρ 4/3 ( r)d r (1)
where F (s) is a function of the reduced density gradient, s, for a given spin. This parameter accounts for the local inhomogeneity:

s = 1 c s |∇ρ| ρ 4/3 (2)
where c s = 2(3π 2 ) 1/3 (note that c s is not the Fermi constant) and the 4/3 exponent in the density ensures that s is a dimensionless quantity.

It assumes large values not only for large gradients, but also in regions of small densities, such as the exponential tails far from the nuclei. Likewise, small values of s occur for small gradients, typical for bonding regions, but also for regions of large density. The homogeneous electron gas is characterized by s = 0 everywhere. When interference between atomic clouds occurs, peaks appear in the s(ρ) plot that reveal the presence of an interaction. Plotting s isosurfaces around these minima, non-covalent interactions appear as very intuitive isosurfaces (see Figure 1). In the NCI method, the sign of the second eigenvalue of the electron density Hessian matrix, λ 2 ,is used to distinguish between atractive and repulsive interactions thanks to a color code: blue for attractive such as hydrogen bonds (HBs), green for very weak interactions such as van der Waals (vdW) and red for steric repulsion. The first two are illustrated in Figure 1.

B. Local analysis: NCI critical points

When plotted in terms of the density, the graph s(ρ) assumes the form of a decaying function f (x) = ax -1/3 , where a is a constant. This can be easily proved from an STO model density. For ψ = e -αr we have ρ = e -2αr and ∇ρ = -2αρ so that:

s(ρ) = 1 c s 2αρ ρ 4/3 = 2α c s ρ -1/3 (3)
Figure 2a shows such behavior for an atomic hydrogen density. 5 When two atoms come However, NCI features do not only appear associated with AIM critical points. NCI interactions also appear in weak closed shell interactions, which do not show AIM critical point (e.g. intramolecular interactions 8,9 ). Indeed, we may differentiate two types of s critical points: 7

1. AIM critical points: critical points of ρ(r), for which ∇ρ(r) = 0.

2. Non-AIM critical points: Points where:

∇ 2 ρ(r) ρ(r) - 4 3 (∇ρ(r)) 2 ρ(r) 2 = 0 (4)
This means that the topology of s recovers all AIM critical points, but also those related to shells 7 (like the One Electron Potential or OEP 10,11 ) and to very weak interactions (for which ∇ 2 ρ( r) >0). The latter cases are usually related to very weak intramolecular interactions, where the balance between geometric constrains and the weak interaction leads to an absence of critical point. To avoid misleading with AIM critical points, we will call "Interaction Critical Points" (ICPs) the points where ∇s = 0, so that we will have AIM-ICPs and non-AIM-ICPs.

C. Global analysis: NCI integrations

In order to integrate NCI volumes, it is necessary to establish a unique definition of the NCI region. Since the difference between the interacting and non-interacting densities is directly reflected in the s(ρ) diagram, it is possible to define the NCI region, Ω N CI , as the domain(s) in 3D space with (ρ, s) values lying in the s(ρ) peak. To identify this region, both the monomer and the dimer densities must be computed and compared. The lowest edge of the monomer s(ρ) curve is splined (green in Fig. 3a) and all points of the dimer s(ρ) plot lying below the splined curve (red in Fig. 3a) are localized in real space (Fig. 3a). 12 Density properties can be integrated within this space to obtain the volume (V N CI ) of the NCI region, or the charge (q N CI ) enclosed within it.

V N CI = Ω NCI d r (5)
q N CI = Ω NCI ρ( r)d r (6)
III. COMPUTATIONAL DETAILS All the geometries in section IV were taken from the s22 set 13 and not reoptimized.

Wavefunctions were obtained with the Gaussian09 package [START_REF] Frisch | Gaussian 09[END_REF] . In order to compare the results from different SCF methods, we have included in our study Hartree-Fock (HF), post-HF methods (MP2, CISD, CCSD) and density functional approaches.

Second-order Möller-Plesset perturbation theory (MP2) is the most affordable wave function-based method beyond the HF approximation that provides an approximate description of all relevant vdW interactions, electrostatics, induction, and dispersion. Furthermore, MP2 is free from spurious electron self-interaction, which leads to noticeable improvements for hydrogen bond description. However, one of the serious shortcomings of MP2 theory is a noticeable overestimation of the dispersion interaction energy. As a reference, we have included the quantum chemistry gold standard: coupled cluster theory with single and double excitations. Perturbative tripe excitations have not been included since they do not affect the underlying electron density and the output of wfn files is not yet implemented in Gaussian. CISD has been included for completeness.

Since nowadays calculations of big systems, where NCIs are extremely relevant, are usually performed with density functional theory (DFT) methods, we have also included two of the most common functionals: B3LYP [START_REF] Becke | [END_REF]16 and B97D. 17 B3LYP is by far the most popular density functional in chemistry, but there is growing evidence showing that B3LYP degrades as the system becomes larger since it fails to correctly describe van-der-Waals interactions. 18 Last years have been marked by an intense research of new DFT variants accounting for dispersive effects. From the comparative benchmark studies emerges that the M06-2X, ω-B97X-D, and B97-D functionals significantly outdo B3LYP in estimating interaction energies in π-complexes. Since these studies are done at fixed geometry, a posteriori corrections become a fair test for the results. Pople basis sets have been used to test the dependence with the basis, including separately polarization and diffusse functions to independently analyze their effect.

The analysis of pseudopotentials was carried out in a methane hydrate of structure sI with were carried out with Quantum-Espresso. 20 The XDM (aka. exchange-hole dipole moment) model was used to correct for dispersion with parameters obtained from Refs. 21,22: a 1 = 0.136 and a 2 = 3.178 Å. Pseudopotentials and projector augmented wave all-electron description of the electron-ion-core interaction 23 were tested. The plane wave energy cutoff was selected to 60 Ry (checking that the total energy is converged at that value) and the reciprocal space was divided into a Monkhorst-Pack 24 grid of 2 × 2 × 2 points. The convergence threshold for the SCF procedure was set to 10 -8 Ry. Geometry optimization was performed using the Broyden-Fletcher-Goldfarb-Shanno 25 algorithm included in the Quantum Espresso package, 20 relaxing lattice parameters and all crystallographic positions.

The system with various phenyl rings in section V (4"'-Methyl, 5"'-Phenyl-[1,1';2',1";3",1"';2"',1"";3"",1""']-sexiphenyl -"sexiphenyl", from now on-see Fig. 7a) was chosen to illustrate the release of π interactions. For the purpose of this work, only two pairs of conformers were selected (not the most stable ones), differing by orientation of a single dihedral angle in each case. Change in this particular dihedral makes either possible or impossible a T-shaped non-covalent interaction. The geometries of these two pairs of conformers were fully optimized at the B3LYP [START_REF] Becke | [END_REF]16 and M06 26 levels of theory with the 6-31G(d,p) basis set. These geometries are provided in the Supporting information.

The Houk-List data were taken from the metadata contained in Ref. 27, which had been calculated at the B3LYP/TZVP/SCRF=DMSO level of theory with and without D3 28 dispersion correction. The availability of these data are greatly acknowledged. [START_REF] Rzepa | data repositories in chemistry and their integration with journals and electronic laboratory notebooks[END_REF] NCI isosurfaces were obtained with a modified version of the NCIPLOT code 30,31 with grids of 0.1 Å along each axis. Electron density at the AIM critical points was obtained with promolden. 32 Non-AIM critical points were obtained from 2D graphs with increments 0.05 bohr along each axis. A new version of the program, which will be soon released (NCIPLOT-4.0), 31 was developed to perform the volume and charge integrations. 30,31 Instead of computing the monomer self-consistently, it computes pseudo-densities (approximated monomer densities from the dimer wavefunction) for each of the monomers and splines the resulting s(ρ) curve. We used the following default values: grids of 0.09 bohr along each axis, and a threshold of 0.05 (0.1 for Houk-List results to include the forming bond) and 0.5 a.u. in the density and s, respectively. In the solid case, Critic2 was used for the NCI calculations. 33,34 In all cases, isosurfaces were visualised with VMD version 1.9.2. 35

IV. THE EFFECTS OF THE LEVEL OF CALCULATION

In order to understand the effects of the method we will analyze a set or representative molecular systems:

• a hydrogen bond (water dimer), a dispersive interaction (CH 4 dimer) and a repulsive clash (bicyclo[2.2.2]octene, from now on bicyclooctene) showing AIM critical points

• formic acid dimer as an example of mixed interactions of different strength (hydrogen bonds and dispersion)

• ethanediol where the hydrogen bond does not have an associated AIM critical point (i.e. it is a non-AIM-ICP).

A stand-alone section is included on the analysis of promolecular (non-relaxed densities)

at the end of the method section. In the analysis of basis sets, a special section will be devoted to pseudopotentials due to their relevance in solid state calculations. Since the main objective of this section is to analyze the effects of the method on NCI results, we will use fixed geometries in all cases. The influence of the geometry on NCI results will be highlighted in section V. 4b,g). Non bonding interactions in bicyclooctene result in density depletion, so that λ 2 > 0 (Figs. 4c,h). It can be observed that both hydrogen bonds and steric clashes appear at greater densities (although different λ 2 sign) than van der Waals (see Figure 4a and4c and Table I). Both strong and weak interactions are present in formic acid dimer (Figures 4d and4i). AIM critical points can be found in Table S1. The case of non-AIM-ICPs is illustrated with ethanediol (Figs. 4e,j). As highlighted in the previous section, it corresponds to a very weak interaction where both positive and negative eigenvalues of the electron-density Hessian matrix are present in the absence of AIM critical points.

A. SCF method

Figure S1 (Supporting Information) shows the s(ρ) behavior of our test systems at the 6-311G level (3D images remain qualitatively the same as in Figure 4). At first sight, the location of the peaks is very little altered, so we have also collected the electron density at the ICP (Table I).

Except for HF, the influence of the electronic structure method is rather small. This is in agreement with previous results for QTAIM. 36,37 In general, peaks at the HF level appear at smaller densities. This is due to the localizing effect of HF. 38,39 Hartree-Fock densities localize the electrons within molecular units, 39 so that the electron density at the ICP diminishes in all the interaction types. Inclusion of correlation, either directly (MP2, CISD)

or by parametrization (functionals) re-establishes values similar to the CCSD reference. It is interesting to note how the method failures are translated into the electron density. MP2 leads to overbinding, with slightly overestimation of densities at the ICP, specially (as noted in the literature) in dispersion interactions. As far as DFT calculations are concerned, they have an effect contrary to that of HF, since functionals tend to delocalize. 38 Thus, the electron density at the ICP is always larger than for CCSD in attractive interactions. The effect however is not as well behaved due to the parametrizations. B97D seems to specially overbind the strong hydrogen bonds, whereas it behaves better in methane dimer (but not in formic acid dimer). The case of bicyclooctene is interesting: both functionals find smaller densities at the ICP than CCSD, which seems to indicate that they underestimate steric repulsion probably due to self-interaction errors.

It is interesting to note that the calculations have been carried out at the s22 fixed geometry, so that the over(under)binding is not translated into the geometry, but only into the electron density. These effects are small and do not yield qualitative (not even quantitative) differences in the NCI picture of the systems. Two main conclusions can be extracted.

On the one hand, global quantities are more affected by the method (see Table II), but the qualitative behavior remains. This is related to the stability of NCI results, and backs up the fact that NCI can be applied in spite of the method used for evaluating the density.

Even non dispersive methods, such as HF or B3LYP give global results that qualitatively agree with CCSD calculations. Changes in the global quantities confirm that NCI does change with the calculation method, but the changes are not very important (even for the volume) and conserve the trends, as long as the geometry is fixed.

On the other hand, the same behavior in ICPs and NCI volumes is observed in all cases.

This is a very interesting observation since the calculation of electron densities at the ICP is much faster than the integration of volumes, so the former can be directly taken as an indicator of the quantitative deviations that should be expected.

From the analysis of volumes, we can also draw conclusions on the visual approach and the choice of isosurfaces. The same s value should be used when comparing various NCIs, both in the same or in different systems, provided the same method was employed to obtain their respective electron densities. This not being the case, different s values are seemingly required to compare NCI diagrams on similar grounds (e.g. from wavefunction, multipolar).

In other words, a shift of cut-offs is needed to obtain comparable images. This is related to the fact that s roughly behaves like ρ -1/3 (see Eq. 3), so that the effect of the method on the density is directly followed by the s-value.

B. No relaxation: promolecular densities

Features of the electron density in the weak interaction region are very stable with respect to the calculation method, to such an extent that these features are already contained in the sum of atomic densities, ρ at i , in what is also known as "promolecular density", ρ pro : 40,41

ρ pro = i ρ at i (7)
A promolecular density obtained from simple exponential atomic pieces is able to predict low-density, low-reduced-gradient regions qualitatively similar to SCF results:

ρ at i = j c j e -r/η j (8)
where c j and η j are adjusted to fit closely spherically-averaged, density-functional atomic densities. 5 Promolecular densities lack the relaxation introduced in a SCF Hartree-Fock or DFT calculation; but are very fast to compute. Thus, qualitative NCI analysis is applicable to large systems, including biosystems, where the description of the interplay of structure and reactivity is crucial. 42 Because the calculation of the electron density in these systems becomes extremely computationally expensive, the promolecular density becomes an attractive option.

Of course, it is important in these cases to understand how the use of promolecular densities affects the NCI calculation. When relaxed densities are compared to promolecular ones, a shift in the peaks is observed (see Figure S2). In general, the promolecular approximation yields bigger values for strong interactions and slight variations in either direction for van der Waals interactions (Table III). The largest shift is observed in the non-bonded overlap cases, whereas it is moderate in the case of stabilizing interactions. This can be understood as the process of relaxing densities at a fixed geometry: steric clashes play the major role, which affects the whole interacting region, and to a smaller extent the stabilizing ones. As an example, the electron density at the peak in bicyclooctene changes from 0.053 to 0.023 upon convergence, whereas it only changes from 0.029 to 0.027 and from 0.018 to 0.014 in water dimer and ethanediol, respectively. Non bonded dispersive interactions also show big relative changes (e.g. from 0.006 to 0.002 in methane dimer), but since these densities are smaller in absolute value, the absolute difference stays negligeble from a visual point of view (see Figure S2). This same pattern is reproduced when the interaction pattern becomes more complicated, as is the case of formic acid dimer.

Since in the promolecular case, the changes are not negligeble, cutoff values need to be changed if we want to obtain similar pictures to SCF ones: the electron density cutoff needs to be bigger and the s isovalue, smaller. The isosurface value needs to be diminished because reduced densities are contracted around the ICPs upon relaxation. This is reflected in Figure 5 for all interaction types. As already observed in the analysis of ICPs, the greatest change is observed in steric clashes, where the repulsive region is greatly reduced upon convergence.

In the same line, all volumes are bigger in the promolecular approximation (Table IV).

Indeed, values start to look of the same order when s=0.4 in the promolecular calculation (except in bicyclooctene), whereas s=0.5 is used in the SCF ones. As a guiding line, the default values for visualization in the NCIPLOT 30 code go from ρ=0.05 and s=0.5 a.u. in SCF calculations to ρ=0.07 and s=0.3 in the promolecular case.

C. The basis set

Figure S3 shows the variation of NCI 2D plots with the basis set on going from double to triple zeta, and on adding diffuse and polarization functions. CCSD has been used in all cases. Linear dependency problems were found with the convergence for Dunning basis sets (and also for 6-311G** calculations of bicyclooctene), so we have restricted ourselves to Pople basis sets as a proof of principle of basis set dependency. As observed in Figure S3, variability upon basis set change is even smaller than upon method variation (Figure S1).

This had already been observed in the study of benzene dimers with Dunning basis sets, where the influence of adding diffuse functions was found to be negligible, since an almost identical distribution of low-gradient spikes was obtained for the aug-cc-VTZ and cc-VTZ basis sets. 43 As in the previous section, 2D and 3D plots were qualitative equivalent (Figure S3), so we have focused on electron densities at the ICP and NCI volumes to quantify changes.

It is important to separate attractive and repulsive interactions here, and to make a joint analysis of ICP density and NCI volumes. As a general trend, in attractive interactions, the increase of the basis set leads to a diminution of the electron density at the ICP. This is clearly observed in water and methane dimers. In both cases, polarization functions have the greatest effect. This can probably be attributed to the match between the method and the base: since CCSD allows for correlation to take place, the use of higher order angular momenta functions is favored. Correlation diminishes bond order in general and reduces the electron density at the critical point, but has different effects on the total NCI volume.

A localized interaction, such as the one in water dimer leads to smaller NCI volumes along with the diminution of ICP density. However, the effect in methane dimer volume is opposite. Due to the nature of dispersion, the main effect of the basis set in methane dimer is a delocalization of the electron density, an increase in its planarity, leading to bigger NCI volumes. So the bigger basis variability reduces the density at the ICP and increases the NCI volume in vdW interactions. A similar effect is observed in ethanediol, although in this case the biggest effect is by inclusion of diffuse functions. Since these two effects are opposite, they are more difficult to analyze when interactions are mixed in the same system, as in the formic acid dimer example. However, the bigger NCI volumes due to an electronic delocalization in the interaction region remain. Finally, repulsive interactions like the one in bicyclooctene yield both an increase in density and volume as the basis set is increased.

All in all, this highlights the need for a compensated basis set, which can have an important role in providing the variability necessary for electronic delocalization. However, this is only relevant for a quantitative description, and does not yield important qualitative differences in the NCI picture (see Figure S3).

Pseudopotentials

The use of pseudopotentials (PP) is widespread in 3D periodic calculations because of the popularity of the plane-wave (PW) approximation, which needs pseudopotentials to avoid the difficult description of the nuclei with plane-waves. If a topological study is going to be carried out on-top of a PP calculation, it is important to take into account some problems that routinely arise. The lack of the core electrons has been reported to produce spurious maxima (specially in light atoms) which affect the AIM topology.

To override this problem the projector augmented wavefunction (PAW) method was developed. Theoretically, with the PAW pseudopotentials it is possible to reproduce the electron density at the nuclei. However, a very fine grid is needed to appropriately plot the steep function in the core region.

Figure 6 shows a comparison of PP and PAW in a methane hydrate with sI structure (178 atoms -8 molecules of methane and 46 molecules of water) taken from 19. We can divide the effect in two regions: core and valence. Most commonly, NCI is used in the region far from the nuclei for identifying non-covalent interactions. As can be seen in Figure 6, all of the above approaches provide a coherent s(ρ) description of the non-covalent interaction region.

However, in some cases, we might be interested in the whole s-topology (e.g. chemical ). This issue can be solved based on the assumption that the core electrons are barely affected by the electronic relaxation so that the core density can be added from a database (also making the calculation lighter). This methodology is included in the CRITIC2 program, 33,34 which we have used to calculate the PP+core and PAW+core curves. This simple treatment solves the issues related with the core-valence region. Moreover, the core shells are also recovered (curves at sign(λ 2 )ρ <0.3, s ∈0.6-1.0). 44

V. OPTIMIZATION OF THE GEOMETRY OR "HOW CAN DISPERSION BE REFLECTED IF THE CALCULATION DOES NOT INCLUDE DISPERSION?"

One of the main questions related to the visualization of non-covalent interactions is related to the visualization of vdW interactions from electron densities that have been obtained from a method which does not contain long range correlation.

In contrast to steric repulsion, it can be seen in Figure 5 that qualitative changes from promolecular to SCF results are very small. This is related to the nature itself of vdW interactions: since dispersion interactions yield a very small perturbation in the system density, it is very well described even at the promolecular level. In the same line of thought, it can also be analyzed to a very good approximation from non-dispersion corrected electron densities. And here a strong "if" comes in: if the geometry is good enough. Hence, in systems where dispersion can become a leading term in the optimization, using a dispersion corrected method can become a must for a good NCI qualitative picture, since it can change the geometry of the system. We have chosen two examples to illustrate this effect:

• Case study 1: a conformational study of a system with various phenyl rings where T-shape interactions can be turned on and off. In this case the relative energies from dispersion-corrected and non-dispersion-corrected functionals agree (though with a big shift).

• Case study 2: two Houk-List diastereoisomers, where including dispersion changes the preferred stereoselectivity.

It should be noted that in both cases, we have chosen conformational studies to make sure that we can compare results. Since gauges appear in between families, conclusions are more difficult to draw in different systems.

A. Case study 1: Dispersion does not change the relative energies

We have chosen a system with various phenyl rings (see Fig. 7a). It has a central toluene ring (ring G) with the other six rings placed in three branches, one with one ring (ring B),

another one with two rings (rings C and V) and the last one with three rings (rings Y, O and R). We have analyzed the effect of rotating two dihedral angles in order to release T-shape A first look at the entry M06-2X in Table VII shows good agreement between relative volumes and the stability of the system, in such a way that the energetic ordering corresponds to the same ordering in NCI volumes (the more stable, the bigger the volume). This is a very interesting example where we can have a first cast at relative energies from a visual approach (even without the integrations -see Figures in S.I.). Moreover, NCI enables to detect the presence of the T-shape interactions in 1 and 3 that explain their greater stability over rotamers 2 and 4, respectively. It should be noted that this kind of structural-stability information is a necessary step in inverse design and in the rationalization of chemistry.

It is also interesting to check the results from different functionals. The lack of dispersion in B3LYP is not able to correctly describe dispersion and T-shape interactions. This gives rise to very small energy differences between the four conformers. In this case, although the ordering among pairs is correct, a deviation appears between the energies and the volumes of compounds 2 and 4 (see Table VII). It is interesting here to determine which accuracy we might expect from this type of estimations. From the NCI volumes point of view, the absence of dispersion in the optimized electron densities is reflected in smaller volumes and smaller differences from one conformer to the other. 2 kJ mol -1 differences in B3LYP-sexiphenyl (see pairs 1-2 and 3-4 in Table VII) are translated into ca. 10 a.u.

difference in the NCI integration. This can give us an idea of the accuracy that one can expect from NCI integrations, and restores the slight mismatch found. In general, whether dispersion is included or not in this example, NCI volumes are able to (i) predict the correct conformational ordering (ii) identify main changes in energetics within chemical accuracy.

Another point that should be raised is that B3LYP volumes do not differ significatly from M06-2X ones, whereas stabilization energies go from ca. 2 to 10 kJ mol -1 within a pair. This is highlighted in Figure S4, where both interaction energies and NCI voluemes are plotted together. Whereas a nice correlation is found for M06-2X, B3LYP energies cluster together in pairs (to T-shape or not to T-shape). This absence of correlation points at an inherent limitation of the NCI approach that should always be considered: NCI volumes should only be used as a relative measure within a family and always following the same calculation method if the geometry is not fixed.

Whereas volumes had a coherent behavior for fixed geometries, this is not the case when the geometry is changed. This is related to the fact that the abilities and deficiencies of a given calculation method for describing interactions will be reflect in the resulting electron density. Section V B serves to this connection in a subtle case.

B. Case study 2: Dispersion changes the relative energies and the geometries

This section illustrates the relevance of including dispersion in the geometry optimization in those systems where dispersion plays a major role in the stabilization energy, such as big non polar systems and transition states.

The Houk-List transition states have been recently revisited with current computational capabilities. It was found that the stability of conformers was given not only by the electrostatics in the NCH δ+ • • •O δ-interaction, but also by the T-shape dispersive interactions between the proline and the phenyk group, specially important in the syn conformers (see Ref. 27 for more details). In this revisit, several models 27 were devised to boost non-covalent interactions altering the diastereoselectivity ratio. The model in Figure 8 was found to even alter the initial ratio of (S,R)/(S,S) -favoring the (S,S) instead of the (S/R)-but only when dispersion was taken into account.

This inversion is related to a change in the geometry, which also leads to differential

NCIs. The relevance of dispersion interactions is only retrieved when dispersion is directly included in the optimization. Indeed, for the B3LYP/TZVP calculations, Table VIII shows the D3 Grimme's correction and the NCI volume are bigger in the least stable case -the (S,S). When the D3 term is included in the optimization (B3LYP-D3/TZVP), dispersion terms are boosted both in (S,S) and (S,R) compounds, as can be seen from the D3 and V N CI without dispersion), subtle effects arise that (i) explain relative stabilities, and (ii) highlight the abilities/deficiencies of constructing a proper electron density with a given method.

VI. CONCLUSIONS

We have first analyzed the effect of the method and the basis set in NCI calculations on representative systems at fixed geometry. When compared to CCSD results, only HF shows a noticeable deviation both in the ICP electron density and the volume of the NCI region.

Other methods yield very similar positions for the ICP, independently of the method and basis set of calculation. The bigger variations are observed for the NCI volumes in complex systems, but they still provide a correct qualitative and semi-quantitative picture. It is important to note that the deviations can be clustered in terms of the method of calculation, which means that both dispersion corrected and non-dispersion corrected functionals give quantitatively the same image, provided the geometry is kept fixed. This highlights the relevance of the geometry in NCI calculations, better than the method used to obtain the electron density. A comparison has also been carried out for promolecular densities. In these cases, noticeable quantitative -but not qualitative-changes appear which reflect the effect of the relaxation on the densities: mainly the steric regions are compacted and reduced densities localize around ICPs.

The fact that NCI quantities are not absolute also stands out from these studies. They can be used to draw conclusions within a given family, but electron densities and NCI volumes cannot be transferred from a system to another. Thus, it is very important to have this limitation in mind when analyzing results. It had already been observed that the electron density at the AIM critical point was dependent on the system, and only comparison among related families were justifiable. 45 NCI volumes, as a generalization of the critical point concept, suffer from the same limitation.

In order to dissect the influence of the geometry on the analysis while staying in the same family, we have carried out two conformational studies where including dispersion has different effects on the relative energies. We have analyzed a conformational case where the changes in the geometry do not lead to changes in the relative energies. NCI volumes were found to correlate within accuracy with the energies for both dispersion corrected and non-dispersion corrected functionals. In other words, the changes in the geometry were followed by electron density changes that mirrored the energy optimization. In this case, we found again that volume values could only be used as relative measures, so that different functionals could not be simultaneously quantitatively compared. In order to introduce even more variability in our tests, we have also analyzed a system where including dispersion changed the geometry and the relative energy of diastereoisomers. Even in this sublte case, NCI volumes and ICPs were able to dissect which stabilizing mechanisms were favored by each functional and the final relative stabilities.

All in all, NCI is very stable with respect to the method of calculation, but it mainly depends on the geometry. From a visual point of view, results are very stable with respect to the geometry, but from the quantitative point of view, NCI is able to reveal changes in the geometry and electron distribution when very subtle effects are at play. As the geometries are allowed to relax on different grounds (e.g. different functionals), subtle effects arise that (i) explain relative stabilities, and (ii) highlight the abilities/deficiencies of constructing a proper electron density of a given method. We suggest that such procedures should be routinely included in protocols for mechanistic exploration.

This geometry/density dependence clarifies how NCI is able to enable the visualization of weak interactions such as van der Waals. These interactions can usually be directly retrieved from the molecular geometry (e.g. the Lennard potential). Since NCI only depends on the electron density, which is known to be fixed to a good extent by the geometry, it directly reflects this fact. The calculation of the energy from these densities (e.g. with functionals) is a step built on the electron density. Indeed, it is common use in functional development to calculate energies on-top of fixed densities calculated at a lower level. This highlights the high accuracy needed for calculating correlation, as well as the ability of electron densities to contain this information, even when they have been obtained from lower levels of calculation.

This suggests that an inspection of the computed NCI surfaces can give rapid visual clues for understanding non-covalent chemistry. Moreover, the visual approach opens the route for inverse design.

FIG. 1 :

 1 FIG. 1: Adenine-thymine complex (geometry taken from Ref. 5) with NCI surfaces (s=0.3, color range: -0.03 to 0.03 a.u.). AIM critical points have been highlighted for clarity: (3,-1) in blue, intermolecular (3,+1) in green and intramolecular (3,+1) in red.

FIG. 2 :

 2 FIG. 2: Behavior of s for a model density ρ = e -αr with α = 0.5288 (hydrogen 5 ). a) s(ρ) atomic b-c) diatomic considering R=1 Å, b) s as a funtion of the H-H distance and c) s(r) for the H 2 molecule.

FIG. 3

 3 FIG.3: a) 2D graph for the monomer and the dimer of formic acid and b) the interacting volume that results from the peak (0.80 a.u. at the CCSD/6-311G level).

Figure 4

 4 Figure4collects the 2D and 3D NCI diagrams of the selected molecules. The hydrogen bond in water dimer shows a negative value of λ 2 at the critical point (Figs.4a,f), whereas λ 2 ≃ 0 (either positive or negative) for van der Waals interactions in methane dimer (Figs.

  FIG. 4: 2D (top) and 3D (bottom) NCI plots of test molecules (computed at the B3LYP/6-311G level of theory). From left to right: water dimer, methane dimer, bicyclooctene, formic acid dimer and ethanediol. The gradient isosurfaces (s=0.5 a.u.) are colored on a BGR scale according to the sign(λ 2 )ρ over the range -0.03 to 0.03 a.u.

aaFIG. 5 :

 5 FIG. 5: Comparison of NCI CCSD (solid) and promolecular (transparent) integrated volumes for test molecules at the 6-311G level: water dimer, methane dimer, bicyclooctene, formic acid dimer and ethanediol.

FIG. 6 :

 6 FIG. 6: 2D NCI graph for a methane hydrate. PP in green, PAW in red, PP+core density in violet and PAW+core density in blue (acronyms like in the text).

reactions 6 )

 6 . In those cases, having a good description of the core-valence separation can become important. The bare treatment with pseudopotentials and PAW leads to the occurrence of false peaks in the s(ρ) diagram (see false peaks in the range ρ=0.2-0.3 a.u.

  FIG. 7: a) numbering of the conformationally relevant dihedral angles ϕ=1 to 6 in sexiphenyl; b-c)pair of structures differing by orientation of the ϕ5 dihedral angle (designated by blue arrow): b) structure sexiphenyl-1, where ring "V" forms a T-shaped contact with ring "Y" and c)structure sexiphenyl-2, where the T shaped contact is absent. d-e)pair of structures differing by orientation of the ϕ3 dihedral angle (designated by red arrow): d) structure sexiphenyl-3, where ring "O" forms a T-shaped contact (red dashed line) with ring "C" and e) structure sexiphenyl-4, where the T shaped contact is absent. The optimized structures are taken from the M06-2X calculations.

FIG. 8 :

 8 FIG. 8: 3D NCI image of the (S,R) and (S,S) diastereoisomers from the Houk-List transition states in the acetone-Pyrrole-2-Carboxaldehyde-Proline model optimized at the B3LYP (top) and B3LYP+D3 (bottom) levels.27 The gradient isosurfaces (s=0.5 a.u.) are colored over the range -0.03 to 0.03 a.u. in terms of sign(λ 2 )ρ.

aa

  Grid increments of 0.09 bohr along each axis and a threshold of 1.0 a.u. and 0.5. in the density and s, respectively, were used in the integrations. b Grimme D3 dispersion correction (Ref.28)

TABLE I :

 I Method dependency of the ICP electron density for test molecules: water dimer, methane dimer, bicyclooctene, formic acid dimer (HB in first column, vdW in the second) and ethanediol.

	Method	H 2 O dimer	CH 4 dimer	Bicyclooctene	Formic acid dimer	Ethanediol b
	CCSD	0.027328	0.0025130	0.023781	0.049234	0.0081382	0.014290 †
	HF	0.025641	0.0023633	0.023113	0.046651	0.0072130	0.013501 †
	MP2	0.027456	0.0025443	0.024157	0.049405	0.0084400	0.014479 †
	CISD	0.027025	0.0024832	0.023535	0.048359	0.0077967	0.014072 †
	B3LYP	0.027572	0.0025944	0.023394	0.049895	0.0082705	0.014206 †
	B97D	0.027872	0.0025270	0.023435	0.050661	0.0083902	0.014206 †
				a

a All calculations with 6-311G basis set. b Due to their similarity (Ref. 8) only attractive point is shown in ethanediol. † Approximated value obtained from 2D plot due to the absence of AIM critical point.

TABLE II :

 II Method dependency of the NCI volume (pseudomolecules as reference) for test molecules at the 6-311G level: water dimer, methane dimer, bicyclooctene, formic acid dimer and ethanediol.

	Method	H 2 O dimer	CH 4 dimer	Bicyclooctene	Formic acid dimer	Ethanediol
	CCSD	0.24	0.86	1.55	0.80	0.36
	HF	0.21	0.79	1.33	0.69	0.28
	MP2	0.25	0.88	1.64	0.84	0.38
	CISD	0.24	0.85	1.48	0.76	0.33
	B3LYP	0.25	0.90	1.50	0.82	0.37
	B97D	0.25	0.88	1.48	0.81	0.37
				a

a Grid increments of 0.09 bohr along each axis and a threshold of 0.05 a.u. and 0.5. in the density and s, respectively, were used in the integrations.

TABLE III :

 III Variation from promolecular to CCSD in the ICP electron density for test molecules: water dimer, methane dimer, bicyclooctene, formic acid dimer (HB in first column, vdW in the second) and ethanediol.

	Method	H 2 O dimer CH 4 dimer Bicyclooctene	Formic acid dimer	Ethanediol b
	CCSD(6-31G)	0.027328	0.0025130	0.023781	0.049234	0.0081382	0.014290
	Promolecular	0.029861 †	0.0063881 †	0.053456 †	0.051855 † 0.0013996 †	0.018327 †

TABLE IV :

 IV Variation of the NCI volume between CCSD (pseudo-densities of monomers are used as reference) and the promolecular approach for test molecules at the 6-311G level: water dimer, methane dimer, bicyclooctene, formic acid dimer and ethanediol. Two different thresholds for s Grid increments of 0.09 bohr, ρ threshold <0.05 a.u. b A threshold of 0.5 in s was used in the integrations for the CCSD volumes. c s <0.6 in bicyclooctene due to the big shift in the peaks. .

	were used in the integrations of the volume for promolecular densities: 0.5 and 0.4.
	Method	H 2 O dimer CH 4 dimer Bicyclooctene c Formic acid dimer Ethanediol
	CCSD b	0.24	0.86	1.55	0.80	0.36
	Promolecular s =0.5	0.63	17.28	9.17	2.22	1.53
	Promolecular s =0.4	0.30	8.05	3.28	0.98	0.22
			a
	a

TABLE V :

 V Basis set dependency of the ICP electron density for test molecules: water dimer, methane dimer, bicyclooctene, formic acid dimer (HB in first column, vdW in the second) and All calculations at the CCSD level b Due to their similarity, only attractive point is shown in ethanediol. † The calculation did not converge due to nearly linear dependency ‡ Approximated value obtained from 2D plot due to the absence of AIM critical point.

	ethanediol.
	Method	H 2 O dimer CH 4 dimer Bicyclooctene	Formic acid dimer	Ethanediol b
	6-31G	0.028449	0.0025736	0.023113	0.049204 0.0075955	0.014496 ‡
	6-311G	0.027328	0.0025130	0.023781	0.049234 0.0081382	0.014290 ‡
	6-311G**	0.023977	0.0023279	†	0.046353 0.0074230	0.014234 ‡
	6-311G++	0.026571	0.0025064	0.024020	0.049072 0.0083187	0.014035 ‡
				a
	a

TABLE VI :

 VI Basis set dependency of the NCI volume (pseudo-densities of monomers are used as reference) for test molecules at the CCSD level: water dimer, methane dimer, bicyclooctene, formic acid dimer and ethanediol. All calculations at the CCSD level. Grid increments of 0.09 bohr along each axis and a threshold of 0.05 a.u. and 0.5. in the density and s, respectively, were used in the integrations.

	Method	H 2 O dimer	CH 4 dimer	Bicyclooctene	Formic acid dimer	Ethanediol
	31G	0.28	0.83	1.28	0.77	0.32
	311G	0.25	0.90	1.50	0.82	0.36
	311G**	0.24	0.97	1.43	†	0.32
	311G++	0.24	0.85	1.57	0.84	0.38
				a
	a

†

The calculation did not converge due to nearly linear dependency.

TABLE VII :

 VII Interaction energies (kJ/mol) and NCI volumes (a.u.) for the four rotamers (1 to 4) of sexiphenyl at the M06-2X and B3LYP levels of theory with the 6-31G(d,p) basis set.Structures 1 and 3 present a T-shape interaction, whereas it has been respectively released in 2 and 4 (TablesS2 and S3in S.I.). Table VII presents the relative energies and NCI volumes for M06-2X and B3LYP.

	Conformer		1	2	3	4
	T-shape		yes	no	yes	no
	M06-2X	E int	0.00	9.46	1.29	8.12
		V N CI	110.54	86.60	100.87	87.74
	B3LYP	E int V N CI	0.00 99.26	2.22 87.33	0.12 98.22	2.49 84.74
	interactions. Rotation of ϕ5 gives rise to structures 1-2 (Figs. 7b and 7c, respectively);
	whereas structures 3-4 (Figs. 7d and 7e, respectively) result from the rotation of ϕ3.

TABLE VIII :

 VIII Relative energy of conformers (S,R) and (S,S) optimized at the B3LYP/TZVP or at the B3LYP-D3/TZVP level (kcal mol -1 ) along with their D3 dispersion correction (kcal mol -1 ) and their NCI volumes (a.u.). Most stable diastereoisomer highlighted in bold.

	Method	B3LYP		B3LYP+D3
	Diastereoisomer	(S,R)	(S,S)	(S,R)	(S,S)
	∆E diast	0.0	2.35	0.10	0.0
	-D3 b	46.910	49.377	47.646	50.303
	V N CI	8.41	10.83	9.06	11.65

Giannozzi P, Baroni S, Bonini N, Calandra M, Car R, Cavazzoni C, Ceresoli D, Chiarotti

Acknowledgments

We want to thank Prof. Rzepa for useful discussions and his always present metadata which make research so much easier. This work was supported partially by the framework of CALSIMLAB under the public grant ANR-11-LABX-0037-01 overseen by the French National Research Agency (ANR) as part of the "Investissements d'Avenir" program (reference: ANR-11-IDEX-0004-02). FIR thanks to FPU program from MECD for a PhD grant and financial support from Spanish Ministerio de Economia y Competitividad and FEDER programs under Projects No. CTQ2012-35899-C02 and No. CTQ2015-67755-C2. M. A. thanks the Fund for Scientific Research-Flanders (FWO-12F4416N) for a postdoctoral fellowship and the Free University of Brussels (VUB) for financial support. The Coimbra Chemistry Centre (CQC) is supported by the Portuguese Fundação para a Ciência e a Tecnologia (FCT), through the project UI0313/QUI/2013, co-funded by COMPETE-UE. I. R. acknowledges FCT for the Investigador FCT grant.