The Volumetric Source Function: Looking Inside van der Waals Interactions

The study of van der Waals interactions plays a central role in the understanding of bonding across a range of biological, chemical and physical phenomena. The presence of van der Waals interactions can be identified through analysis of the reduced density gradient, a fundamental parameter at the core of Density Functional Theory. An extension of Bader’s Quantum Theory of Atoms in Molecules is developed here through combination with the analysis of the reduced density gradient. Through this development, a new quantum chemical topological tool is presented: the volumetric source function. This technique allows insight into the atomic composition of van der Waals interactions, offering the first route towards applying the highly successful source function to these disperse interactions. A new algorithm has been implemented in the open-source code, CRITIC2, and tested on acetone, adipic and maleic acids molecular crystals, each stabilized by van der Waals interactions. This novel technique for studying van der Waals interactions at an atomic level offers unprecedented opportunities in the fundamental study of intermolecular interactions and molecular design for crystal engineering, drug design and bio-macromolecular processes.

Recently, Johnson et al. [26][27][28][29] developed an approach to visualize NCIs based on topological analysis of the electron density, dubbed NCI surfaces. Within their method, the reduced density gradient (RDG) is analyzed and an interaction map is generated where → s r ( ) 0; i.e. critical points. The magnitude of ρ r ( ) is used to scale the relative strengths of these interactions, and the nature of the interaction is identified by analysis of the Hessian of the Laplacian of electron density. This analysis leads to a highly visual representation of NCIs, without need for expert knowledge in electron density analysis. The intuitive nature of this approach has made it immensely popular among both the experimental and theoretical communities. In particular, NCI surface analysis has found widespread applications in the materials, chemistry, and biology communities, offering new insight into the structure of bonding interactions [26][27][28][29] and self-assembly phenomena 10,[29][30][31] . Furthermore, this technique has allowed for the identification of non-conventional hydrogen bonds, which are typically unobserved within the frameworks of other theoretical approaches [31][32][33][34][35] .
Despite its conceptual power, analysis of NCI surfaces does not permit atomic-level interpretation. That is to say that while the NCI surface shows the entire interaction surface between molecules, no information is obtained as to which atoms contribute to this interaction. Such information would offer invaluable insights into the tunability of NCIs. A topological tool -the source function (SF) [36][37][38][39][40][41][42][43][44][45] -is already known to provide this insight, doing so within the framework of the Quantum Theory of Atoms in Molecules (QTAIM) 43 . The NCIs are typically descripted by SF at well-defined critical points. Albeit, the SF reconstruction along a line, on a surface or within a volume were previously introduced and discussed, both for the electron density and the electron spin density [46][47][48] . Here, we describe an interesting application of the SF double integration procedure through its combination with the study of NCI surfaces. This extension allows, for the first time, insight into the atomic contributions to NCIs, and hence represents a new tool in the targeted design of molecular recognition. Through this work, we develop a fundamental extension to the traditional bond critical point 43 (BCP) seen by QTAIM, dubbed the vdW volume (V vdW ). Extension of the SF to accommodate the V vdW leads to generation of a further development: the volumetric source function (VSF). Analysis of NCI surfaces by the VSF provides, for the first time, a chemically intuitive, quantitative approach to the study of non-directional interactions. This technique adds a new method to the chemists' toolbox to investigate the nature and structure of vdW interactions, with unique atomic-level insight.
To demonstrate the potential of our extension to QTAIM a brief introduction to the fundamental theories and developments is provided, which is followed by the description of the corresponding algorithmic developments. The development of theoretical considerations is tested on acetone, adipic, and maleic acids molecular crystals discussing the interpretation of vdW interactions with the new topological tool so called VSF here presented.
theoretical Background introduction to QtAiM and the source function. QTAIM is based on the separation of molecules into distinct atomic subunits, familiar to chemists. These atomic subunits (known as atomic basins) are delimited by surface so called zero flux surface. Such surface is made by an infinity of points r, for which the dot product of the gradient of electron density, ∇ρ, and the normal vector, n, is zero (zero-flux boundary condition).
This is to say that the surface of an atom is defined by the set of vectors n ( ) that align perpendicular to these points of ∇ρ = 0. All points which sit on this surface are denoted r s . Thus, there is no change (or flux) of electron density across this surface, and it is 'zero-flux' . Particular critical points (namely where ∇ ρ r ( ) 2 is negative in only one direction) are known as BCPs and define the unique bond path of covalent (or semi-covalent, e.g. HB) interactions. In order to describe bonding at BCPs intuitively, an extension to QTAIM was developed, dubbed the SF. The SF describes the contribution of each atomic basin to the BCP and hence describes the atomic contributions to individual covalent bonds 43 , or intermolecular hydrogen bonds 44,45 . Only in very specific cases the SF can be used with success to describe bonding in vdW interactions 45 , because in principle vdW interactions may be described in terms of SF reconstructions as for any other point or surface or volume, but there might arise serious numerical accuracy problems. The electron density at point r that comes from atomic basin Ω is given by the sum of two contributions: (i) the integral of ′ ∇ ρ r ( ) 2 evaluated over all points ′ r within Ω, each weighted as an inverse function of its distance from the point of interest − ′ − r r 1 ; and (ii) the flux of the electric field density, ε − r r ( ) s , across the boundary of Ω, calculated at r. This term depends on the electron density of the surface of Ω, ρ r ( ) S [36][37][38][39][40][41][42][43][44][45] . The electron density at point r within an atomic basin Ω can be written as, S S 2 For a closed system with boundary at infinity the second term of Equation 2 is reduced to zero. For a system with more than one atomic basin, ρ r ( ) becomes a sum of the contribution of each atomic basin, with each individual contribution known as the SF for the respective basin, 2 The integrand in Equation 3 is defined as Local Source (LS).

non-covalent interaction surfaces. The RDG
For regions far from the various nuclei of a system (i.e. where density decays exponentially to zero), the RDG adopts large positive values. In contrast, the RDG values approach zero for regions of covalent and non-covalent bonding. Hence, the magnitude of s offers a good indication of the position of NCIs. The nature of the NCI is subsequently defined by analyzing the Laplacian of electron density ∇ ρ r ( ) 2 . To characterize interactions, the Laplacian is often decomposed into its three principal axes of maximum variation, the three eigenvalues λ n of the Hessian of electron density, ∇ ρ = λ + λ + λ λ ≤ λ ≤ λ r ( ) ( ) 2 1 2 3 1 2 3 . As per convention, λ 1 , λ 2 and λ 3 are in ascending order. Then for BCPs it is clear that λ 3 is always greater than zero. This was arbitrarily assumed, in view of the fact that the low RDG isosurfaces bound volumes where in most cases RDG goes to zero and so enclose an electron density critical point (CP). When this is a BCP (bonding interaction) λ 2 is negative and when this CP is a ring CP (associated to non-bonding/repulsive interactions at the ring CP between the atoms forming the ring) λ 2 is positive and for vdW interactions λ ≈ 0 2 . Hence the nature of interactions at points of → r s( ) 0 can be defined by the corresponding value of ρ ⋅ λ r ( ) sign( ) 2 . It is generally accepted that the magnitude of ρ r ( ) corresponds to the relative strength of interaction. As a rule of thumb, it has been suggested that the limits ± .

Development of volumetric source function
The above discussion of SF relates to the decomposition of atomic basin contributions to a well-defined point. This is generally taken to be a BCP, and hence SF has traditionally been used to interpret covalent and well-defined NCIs such as HBs. Unfortunately, due to their non-directionality, vdW interactions are not associated with such well-defined critical points. Hence, these interactions have remained beyond the scope of QTAIM analysis. In order to extend QTAIM for this purpose, we define a V vdW via the infinite points within the volume defined by NCI analysis, ∈ V r vdW v dW . The number of electrons within V vdW for a closed system with boundary at infinity is given by the sum on all M atomic basins of VSF. Such equation is defined as the integral of SF on the V vdW giving the contribution of an atomic basin to such volume, The general workflow developed for calculation of the VSF is summarized in Scheme 1. In order to investigate the nature of vdW interactions within the solid state, the input structures were first relaxed in the solid-state using plane-wave Density Functional Theory (PW-DFT). The electron density of the relaxed structure was subsequently obtained with a plane wave kinetic energy cut-off of 100 Ry. This value has been previously shown to facilitate the reliable calculation of the SF 48 . The resulting electron density was analyzed for NCI surfaces using the NCIplot code, as implemented in CRITIC2 49,50 . The NCI surfaces were generated for ρ − .
≤ ⋅ λ ≤ . 0 02 s ign( ) 0 02 2 a.u. The NCI surfaces corresponding to vdW dimers that compose the reticular motif of acetone, adipic and maleic acids crystal structures were subsequently isolated and taken as definition of the ... To assess the suitability of PW basis sets, the electron density in the dimers was recalculated with localized basis set (LBS) at the MP2 51-55 /aug-cc-pVTZ 56 level. This density was subsequently projected onto the NCI surface obtained above (i.e. onto the V vdW ). The SF for each atomic basin on each sampled point within V vdW was calculated using YT integration scheme 57    The error was found to achieve the asymptote for 1000 points. Thus, to guarantee the most feasible NCI description V vdW were sampled with more than 1000 points neglecting the point with ρ less than 10 -6 e/bohr3. Scheme 1. Workflow diagram describing the methodology used for correct reproduction of the VSF using localised basis sets (LBS) within the MP2/aug-cc-pVTZ and plane waves (PW) electron densities. The function f1 is the Gatti's reliability parameter 42 , which describes the percentage error in the reconstruction of the electronic density within the vdW volume (V vdW ).  63 . The PW86PBE 64,65 exchange-correlation functional was used in combination with the exchange-hole dipole moment (XDM) dispersion correction with damping factors a 1 =0.6836 and a 2 =1.5045 [66][67][68] . XDM uses the interaction of induced dipoles (and higher-order multipoles) to model dispersion: the source of the instantaneous dipole moments is taken to be the dipole moment of the exchange hole [66][67][68] . The wavefunction was expanded in a plane wave basis set to 100 Ry, and Brillouin zone integration was performed on 6×6×6 Monkhorst-Pack k-point grid 69 . The SCF convergence was accepted when less than 10 -8 Ry, and geometry relaxation accepted once residual atomic forces fell below 10 -8 Ry/bohr. Comparison of the experimental and relaxed structures are given in Table  E.S.I. 1. The NCI surfaces were generated from the PW density, using the NCI implementation within CRITIC2. Values of . ≤ ≤ . s 0 45 055 for computational generated densities as suggested in the original NCI work 26-29 were chosen. The VSF was calculated using an in-house code, according to the work scheme shown in Scheme 1. The LBS electron density used in the final integration step was obtained for the isolated dimers using Gaussian v.16 70 at the MP2 51-55 level with aug-cc-pVTZ basis set 56 for each atom. The dimers extracted for electron density analysis were used as input for SAPT 25 analysis of NCI energies. The scaled SAPT zero energies were calculated within the PSI4 software using the jun-cc-pVDZ basis set 71 , the bronze level SAPT method identified by Parker et al. with overall error of ± 2.05 kJ/mol 72 .
To offer a more direct numerical comparison the VSF values in this article they are provided in term of VSF%,
As analogue to the role of ρ bcp in HB systems 44,45 , one would expect that the value of n V e vdW should correlate to the total interaction energy across the dimer. As previously explored by Saleh et al. 73 (2015) the correlation between the value of electron density and the interaction energies calculated with SAPT 25 interaction energies, Table 1 and once again we do not find any notable correlation. This suggests that -unlike ρ bcp -the integral of electron density has different character and is not directly indicative of the strength of this category of intermolecular interaction. This is presumably due to the fact that a V vdW composed of points with λ ≈ 0 2 contains loci of both attractive and repulsive character. www.nature.com/scientificreports www.nature.com/scientificreports/ VSf from localized basis set (LBS) electron density. In order to consider the validity of the VSF, we adopt the approach suggested by Gavezzotti for electron density decomposition 24 . Here, the nuclear positions are extracted from the relaxed crystal structure and the electron density is recalculated using a LBS with the post Hartree-Fock correlation scheme, MP2. This approximation is supported by the fact that the Laplacian of electron density exponentially decreases far away from the nucleus as − e x 2 , where x represents the relative distance from nucleus at the point R 43 . The dimers can hence be extracted from crystals given that the primary neighbours will negligibly affect the intermolecular interactions between dimer pairs. To best correlate our values to crystalline structures, the electron density obtained using the MP2/LBS approach were projected onto the V vdW s calculated from periodic calculations. Decomposition of LBS-based ρ V vdW via  Table 1. Number of electrons, n V e vdW , within vdW volumes (V vdW ) coming from local basis set electron density for each of the three vdW dimers: acetone (AC), adipic acid (AA) and maleic acid (MA). Note that symmetryadapted perturbation theory (SAPT) energies are given in kJ/mol. The SAPT energy is decomposed into electrostatic (ele.), exchange (exch.), inductive (ind.) and dispersion (disp.) contributions. www.nature.com/scientificreports www.nature.com/scientificreports/ we obtain values f1 of 0.14%, 1.19% and 1.32% for acetone, adipic and maleic acids, respectively. Considering the VSF more closely, we simply decomposed n V e vdW into its atomic contributions for each dimer starting from the electron density contribution of each atom inside the V vdW , Fig. 2. In each structure there is a single oxygen atom which clearly dominates in its contribution to n V e vdW . Hence, there is a single oxygen atom in each case which can be suggested as being largely responsible for the vdW interaction. There is no corresponding density sink visible in any of these three systems. Instead, it appears that the electron density donated by the oxygen atom is spread largely uniformly across the remaining atoms in the system. This is consistent with prevailing theories of vdW interactions, suggesting that the interaction is indeed dependent on the total number of atoms into which the donated density can distribute.
VSf from plane waves (pW) electron density. Recent work has, however, demonstrated that the SF could be reliably reproduced from PW densities, provided sufficiently high kinetic energy cut-offs 48 . It is not immediately evident whether the same holds for reconstruction of regions with very low electron density, as V vdW , is possible and if from SF evaluated in such regions is possible to calculate VSF. We therefore extend the above discussion to consideration of VSF as obtained by a PW basis, with contributions to the V vdW taken from atoms throughout the entire structure, i.e. throughout the unit cell. Hence, the n V e vdW will be different with respect to that calculated from LBS electron density of extracted dimers. As the case for LBS electron densities, the f1 values for PW electron densities were obtained for acetone, adipic and maleic acids (i.e., 22.18%, 26.33% and 30.22% respectively). As expected from the above discussion, f1 values for PW electron densities are highest than for LBS electron densities.
As compared with the atomistic VSF contributions obtained by LBS above, the qualitative picture obtained for acetone is well reproduced, Fig. 3A. Again, only a single oxygen atom contributes to the NCI surface, while the remaining atoms exhibit approximately equal acceptor contributions. Unfortunately, the systems which exhibit much lower values of electron density and consequent n V e vdW is poorly reproduced when the VSF is generated by PW densities, Fig. 3B,C. Hence, the attempt to use PW densities to derive SF values at the V vdW was seen to not be so feasible like the one with LBS densities. Albeit in previously works SF values derived by PW densities were seen comparable with those obtained with LBS densities at the BCP 43 . This can be due to the fact that the reconstruction of electron density into a point is amenable of the value that it itself assumes at that point. Thus, for electron densities in the order of 10 -3 as seen at the BCPs PW densities work well, but for electron densities less than 10 -3 e/bohr 3 PW densities are not enough even if they are taken in large numbers.
Furthermore, a general note it is worth mentioning that the atomistic VSF contributions, Tables E.S.I. 5-7, suggests that atoms which exists within unit cell but outside of the dimer pairs contribute on average no more www.nature.com/scientificreports www.nature.com/scientificreports/ than 1% to n V e vdW . Thus, it is clear that all atoms within the periodic structure do contribute to the structure of the V vdW , but the atomic contribution for those are not directly involved is reasonable neglectable confirming ulterior the work of Gavezzotti with PIXEL 24 .

conclusions
Here we outline a different approach to deconvolve vdW interactions in terms of their atomic composition, which description has been a long-standing challenge for theorists. This approach is rooted in the QTAIM 56 and represents a vdW analogue to the previously developed SF analysis. This development includes the description of a V vdW and the corresponding VSF. Based on the approach for electron density decomposition proposed by Gavezotti 24 , it was found that electron densities based on the localized basis set allowed reliable calculation of VSF. Analysis in this way suggests that vdW surfaces are generated through the charge donation of a single heteroatom, here oxygen. The reciprocating charge accepting behaviour is shared across the remaining atoms, consistent with prevailing theories of vdW interactions. This suggests that the VSF method indeed presents a consistent picture of these weak interactions. Despite correctly reproducing the number of electrons, n V e vdW , inside the V vdW , plane wave basis sets were found less capable of providing reliable calculation of VSF. This was particularly true for systems with low values of electron density within the V vdW , and is presumably due to the inaccuracies associated with reconstructing points of low electron density. However, this work reports on a new approach for the quantum chemical topological investigation of weak vdW interactions. The methodology is based on empirically definition of vdW interaction using a fundamental dimensionless quantity in DFT that describes the deviation from a homogeneous electron distribution so called RDG [73][74][75] . This leads to atomistic-level detail of the interacting surface and hence offers the first approach for rational design of these ubiquitous interactions.