A Topological Data Analysis Perspective on Non-Covalent Interactions in Relativistic Calculations

TL;DR

This paper applies Topological Data Analysis to electron densities in relativistic calculations, revealing that previous claims of non-covalent interactions due to relativity are not supported by topological persistence measures.

Contribution

It demonstrates the use of TDA in analyzing topological features of electron densities and challenges prior interpretations of relativistic effects on non-covalent interactions in a specific gold-containing molecule.

Findings

TDA shows no significant topological features indicating non-covalent interactions due to relativity.

Persistence values decrease with larger basis sets, suggesting features are not meaningful.

Reduced density gradient analysis supports the absence of relativity-induced non-covalent interactions.

Abstract

Topological Data Analysis (TDA) is a powerful mathematical theory, largely unexplored in theoretical chemistry. In this work we demonstrate how TDA provides new insights into topological features of electron densities and reduced density gradients, by investigating the effects of relativity on the bonding of the Au4-S-C6H4-S'-Au'4 molecule. Whereas recent analyses of this species carried out with the Quantum Theory of Atoms-In-Molecules [Anderson et al., Chem. Eur. J. 25, 2538, 2019] concluded, from the emergence of new topological features in the electron density, that relativistic effects yielded non-covalent interactions between gold and hydrogen atoms, we show from their low persistence values (which decrease with increased basis set size) these features are not significant. Further analysis of the reduced density gradient confirms no relativity-induced non-covalent interactions in Au4-S-C6H4-S'-Au'4. We argue TDA should be integrated into electronic structure analysis methods, and be considered as a basis for the development of new topology-based approaches.

Equations

$$s(\vec{r})=\frac{1}{2(3{\pi }^2{)}^{1/3}}\frac{|\nabla \rho (\vec{r})|}{\rho (\vec{r}{)}^{4/3}}.$$

$$p_{BCP-RCP}$$

$$\{\begin{array}{ll}|log({\rho }_{BCP})-log({\rho }_{RCP})|& \exists \text{leavevmode}\text{BCP-RCP}\\ 0& \nexists \text{leavevmode}\text{BCP-RCP}\end{array}$$

Full text

Keywords — noncovalent interactions, relativistic effects, quantum chemistry, topological data analysis.

1 Introduction

A well-exploited consequence of the Hohenberg-Kohn theorem[43] is that the exact, non-degenerate, ground-state electron density (ED) contains complete molecular information and therefore encodes all interactions in molecular systems. In order to retrieve information on interatomic interactions from the ED, this property is subjected to various analysis schemes, which often involve visual exploration. This is, however, both a technical challenge and a test of chemical concepts. The former relates to a nontrivial choice of analysis techniques of scientific data. For instance, direct visualization schemes (an example of which is representing 3D scalar fields as isosurfaces), although intuitive and simple in realization, often lead to images that are cluttered and difficult to grasp, especially when aimed to illustrate complex scientific questions. On the other hand, the so-called feature-based visualization methods are designed to recapture only the meaningful features of data, however are more complex to comprehend and therefore not as widely used.[48, 87]

Topological Data Analysis is a recent research field, at the interface between mathematics and computer science, which suggests to investigate data based on its structure.[41] Rooted in sound theoretical settings such as Morse theory,[52] TDA provides to users a set of tools for the robust extraction of topological features in their data. This includes for instance critical points, integral lines, separating surfaces, voids, etc. Moreover, Persistent Homology (PH) [31] brings an appealing framework for measuring the salience of topological features in the data, which is well established both at a theoretical and practical level.[55] In practice, it allows users to discriminate important features from non-significant configurations and provides the necessary basis for multi-scale data analysis.[85] Recent years have witnessed some successful applications of TDA to chemistry, where it helped to understand hydrogen-bonding networks in ion and molecule aggregates,[100, 101, 81] aqueous solubility of molecules,[99, 61] stability of fullerenes,[102] molecular transition pathways,[14] conformational spaces of molecules[50] or the bonding patterns in molecular systems.[95, 38, 42, 15]

In this manuscript we focus on the use of TDA to investigate non-covalent interactions (NCIs) for systems containing heavy elements. The description of NCIs is by itself a challenge for existing chemistry concepts, and has been a driving force for the development of new analysis techniques.[56, 22] For systems containing heavy elements, there is the additional challenge of including the effects of relativity in the electronic structure as these can qualitatively change the properties of a molecular system.[69, 78, 79, 27, 70] The bonding analysis in such systems is especially demanding if the spin-orbit coupling (SOC) is non-negligible.[26] There is currently intense activity in the development and practical use of new analysis tools and methods that can reliably characterize chemical bonding irrespective of the strength of SOC, and as we hope to show, TDA is a promising approach.

In order to illustrate and further extend the idea of using TDA for interatomic interactions, we have chosen a molecular system studied recently in ref.[7], Au4-S-C6H4-S’-Au’4 (1), which is interesting in the context of gold-sulfur[39, 58, 103, 72, 73] and gold-hydrogen[71, 83, 76, 77, 93] interactions due to their ubiquity in metalorganic complexes and on interface - a proof of many faces of the fascinating chemistry of gold.[18, 66, 67, 68] A remarkable finding of the study of Anderson and co-workers[7] is that it characterized the emergence of NCIs between gold and hydrogen centers, using the Quantum Theory of Atoms-In-Molecules (QTAIM).[10, 11]

QTAIM is a very popular analysis tool,[64] which describes the stationary points of the ED gradient and interprets these points as chemical objects - maxima correspond to nuclei (although some molecules exhibit non-nuclear attractors[36]), while saddle points are identified as bond critical points (BCP), ring critical points (RCP) or cage critical points (CCP). ED gradient paths emanated from a saddle point, terminating at two maxima and traced along maximum ED, originally named bond paths,[9] are interpreted as interaction lines between atoms and the areas of space delimited by zero-flux surfaces of ED gradient define atomic basins. However, recent ambiguous results raised a discussion on the limitations of this method,[21, 63, 49, 46, 35, 96, 44, 62, 19, 97] an important point in this debate being precisely the failure of QTAIM to study NCIs.

While it can be particularly appealing to use QTAIM in this context, especially due to its availability in many standard quantum chemistry codes and supposedly simple interpretation of topological features of ED, it is known that the latter is not a sensitive descriptor of subtle interatomic interactions, due to its exponentially decaying behavior away from the nuclei. This observation has led to the development of other descriptors, better suited for studying non-covalent interactions, such as the NCI index[45, 24, 54, 53] based on the reduced density gradient (RDG), $s(\vec{r})$:

[TABLE]

At critical points of ED where $\nabla\rho(\vec{r})=0$, RDG reaches its minimum value, $s(\vec{r})=0$, therefore tracing the minima of $s(\vec{r})$ can be used to locate the interaction sites in molecular systems. The minima of $s(\vec{r})$ may also appear in areas not associated with QTAIM critical points - referred to as "non-AIM critical points" or "interaction critical points (ICPs)"[23, 17] - in which case they indicate very weak interactions.

In practical applications, both ED and RDG are approximate and their quality depends on the quantum mechanical (QM) model adopted in calculations. The approximations introduced with this model,[40, 65, 74] involving the choice of the Hamiltonian (determining whether and which relativistic effects are included), the method (responsible for the description of electron correlation) and the basis set (establishing the quality of one-electron functions), as well as other numerical parameters (such as the numerical thresholds, grids), can be assessed in a systematic way, yet there are relatively few studies illustrating how all of them affect the aforementioned fields and their topological features.[32, 60, 20, 6]

Inaccuracies of QTAIM results also stem from theoretical and semantic reasons. The former is related to the fact that QTAIM uses concepts from the graph theory, hence inherits its limitations, in particular it can give the "static" or "binary" information about topological features of ED, but not on their evolution, stability or importance (contrary to TDA). The latter reason is a source of continuous debate, related to the misinterpretation of QTAIM concepts (an example of which is a confusion between "bond paths" and "chemical bonds"[12, 80, 62, 25]), and to the overuse of the adjective "topological" in chemistry which should not be confused with mathematical topology.[8, 51]

The TDA approach to the interatomic interactions is more general. TDA provides information on the topological features of a scalar field, including their properties, configuration, evolution and significance, while staying oblivious to what the scalar field represents. This detachment of a tool and an interpretation of its output is very helpful to gain a new perspective on the research problem or to redefine existing concepts in a more robust way. In particular, TDA applied to ED, RDG and a signed electron density ($\tilde{\rho}(\vec{r})=sign(\lambda_{2}(\vec{r}))\rho(\vec{r})$ with $\lambda_{2}$ being the second eigenvalue of the ED Hessian)[24] enabled to automatically extract and classify interatomic interactions merely based on the topological persistence (vide infra) of critical point pairs of these fields.[38]

2 Computational details

Following ref.[7], the geometry of the Au4-S-C6H4-S’-Au’4 molecule was optimized in the ADF software[84, 13, 34] using three different models: (i) scalar-relativistic Zeroth-Order Regular Approximation (sr-ZORA),[90, 89, 91, 92] PBE functional[59] and TZP basis set[88] at the singlet spin-restricted level (denoted throughout the text as g1), (ii) sr-ZORA, PBE, TZP at the triplet spin-unrestricted level (g2) and (iii) sr-ZORA, PBE, QZ4P[88] at the triplet spin-unrestricted level with tighter SCF convergence criteria (g3). The study of a triplet state of this molecule, motivated by its proximity to the singlet ground state, was conducted to explore how the spin state affects the ED topology.

We considered simpler model systems of that molecule to investigate the importance of a local environment to the topological features of interest, Au-S-CH-CH2 (1). Its geometry was obtained (without reoptimization) by removing atoms from the g1 model of Au4-S-C6H4-S’-Au’4 molecule, so that it contains only one gold atom linked by sulfur to the ethene group that replaced the benzene ring. As the positions of all nuclei in Au-S-CH-CH2 were preserved as in the original molecule and due to the lack of symmetry in Au4-S-C6H4-S’-Au’4, this simplification resulted in two geometries of the small model (atom numbers in parenthesis, after ref.[7]) - Au(4)-S(1)-C(1)H-C(2)H(1)H ("a") and Au(8)-S(2)-C(4)H-C(3)H(2)H ("b"), hence the nomenclature m-a-g1 and m-b-g1, respectively.

For all studied molecules, ED has been evaluated in single-point calculations with the ADF software. These employed the non-relativistic (NR) Hamiltonian as well as scalar and spin-orbit ZORA Hamiltonians (sr-ZORA and so-ZORA, respectively), DFT with different exchange-correlation (XC) functionals (PBE, B3LYP,[82] SAOP,[75, 37] M06-2X[105, 104]) and three basis sets (DZ,[88] TZP and QZ4P). Additionally, for the small Au model we performed calculations in the locally-modified version of the DIRAC software[3] with the Lévy-Leblond (LL),[47] spinfree Dirac-Coulomb (SFDC)[28] and Dirac-Coulomb (DC)[94] Hamiltonians, two methods involving different levels of electron correlation (Hartree-Fock (HF) and DFT with the PBE functional) and a basis set of triple-zeta quality.[29, 30, 4, 98]

The TDA of the logarithms of ED and RDG ($log(\rho)$ and $log(s)$, respectively[1]) was performed with the TTK software package.[86] All details can be found in the electronic supplementary information (ESI).

3 Results and discussion

Even though the optimized geometries of the Au complex are slightly different than in ref.[7], our analysis confirms the main conclusions of the cited work: additional pairs of critical points of ED (1-saddles and 2-saddles, interpreted as "RCPs" and "BCPs", respectively) between Au and H nuclei appear for certain combinations of Hamiltonians, methods, basis sets, molecular geometries and spin states, but they never appear in the non-relativistic context.

The comparison of topological features of relativistic (so-ZORA) and non-relativistic EDs of g1 molecule obtained with the PBE/QZ4P QM model is illustrated in 2, where, in particular, the topological skeleton of non-relativistic ED (without these additional saddle-saddle pairs) and relativistic ED (with these pairs) are shown. The results from TDA of these densities are described in detail in the caption to 2.

The TDA pattern of critical points and 1-separatrices of ED is principally the same as the one obtained from QTAIM analysis.[38] However, as we shall show, TDA of ED calculated with various QM models paints a more nuanced picture.

Before continuing on that thought, we shall clarify that the presence of 1-saddles accompanying 2-saddles between Au and H atoms is a result of Morse inequalities[52] 111Equivalently it can be deduced from the Poincaré-Hopf theorem, which applies to generic vector fields. and is in accordance with chemical intuition due to the ring-shaped geometry of the Au-S-C-C-H chains. In Persistent Homology, the critical points of a scalar function are paired in a non-ambiguous way (called the Elder rule[31]) and in this particular case, the additional 1-saddle and 2-saddle points, if present, form a pair on each side of the aromatic ring. The lifespan of this pair is measured by persistence, calculated as a difference of scalar field values in corresponding critical points. Accordingly, we define the persistence for the additional "BCP-RCP" pair as the difference of the (logarithm of) densities in the position of relativistic "BCPs" and "RCPs" as: 222 The mathematical background of persistence can be found elsewhere.[31]

[TABLE]

and use this definition to calculate the persistence values of extra saddle-saddle pairs of $log(\rho)$ determined with various QM models. In all cases these values are very small, constituting typically less than 0.1% of the full persistence range (represented by red bars on the persistence diagram plotted for relativistic ED in 2). Furthermore, they decrease with an increasing accuracy of the QM models, notably in molecules g1, g2 and g3, these persistence values systematically reduce along the DZ - TZP - QZ4P basis sets order for all XC functionals (3), a trend which is is observed for both spin states and all geometries. The most remarkable change is seen for the combination of QZ4P basis and the meta-GGA M06-2X functional, for which these saddle-saddle pairs disappear (as signified by the zero-persistence values on the plots for all but the Au4-H1 pair in the g2 system. This functional is also the one for which persistence values are generally lower. Taking at face value the claim of its superior performance for describing non-covalent interactions compared to the other functionals used in this work,[105, 16] our findings that this functional gives zero-persistence values would therefore imply the absence of these topological features.

At this point it is also interesting to compare results for g2 and g3, since both represent the triplet state of the Au4-S-C6H4-S’-Au’4 molecule, only with a larger basis set and tighter SCF criteria used in the geometry optimization and ED calculations of the latter. It turns out that these improvements in the calculations have a tremendous effect on the persistence values of extra saddle-saddle pairs on both sides of the aromatic ring, namely these pairs are absent in g3 in the most accurate setups (QZ4P basis set) for all selected XC functionals. This confirms the conclusion from ref.[7] that the topology of ED significantly changes depending on how close the structure is to a minimum of the Born–Oppenheimer (BO) energy surface, but also underscores that this pair is not as significant as the conventional QTAIM analysis would make it to be.

Similar analysis of m-a-g1 and m-b-g1 (3), this time focusing on the comparison of different Hamiltonians and methods, reveals the strong influence of the electron correlation on the persistence values of the additional saddle-saddle pair, seemingly counteracting the relativistic effects. This appears from the DIRAC calculations involving the fully-relativistic 4-component DC Hamiltonian with HF and DFT/PBE methods: the ED calculated at HF level does not exhibit the extra saddle-saddle pair, while this pair is present in the DFT calculations and its persistence is (slightly) larger than when an approximate relativistic Hamiltonian (SFDC) is used. This observation, together with the negligible role of the spin-orbit coupling on the persistence values (as evaluated on scalar-relativistic geometries), suggests further studies on m-a-g1 and m-b-g1 with wavefunction-based methods at the scalar-relativistic level of theory. Apart from a change in electronic structure method, it may well be that changes in equilibrium structures take place with the consideration of spin-orbit coupling in the geometry optimization procedure, as evidenced by a recent study on SOC-triggered bond contraction in gold clusters.[33] We have nevertheless opted in this work to focus on selected geometries such as those previously considered by Anderson et al.[7] This was motivated, first, by the comparison of the results for these model systems with the original molecular system, which indicates that the chemical surrounding of interacting Au and H atoms does not necessarily affect the topology of ED in these Au-H interatomic areas. And, second, that we considered that this would divert attention to what we believe is the central point of this work, which is to underscore the differences and similarities between TDA and QTAIM. We note, however, that further studies on these gold systems should indeed consider the effect of structural chances brought about by SOC.

In TDA the consensus is that the persistence of a topological feature is associated with its importance (in particular in the presence of additive noise), therefore it is frequently used in the topological simplification of complex scalar fields. Features manifested by so small persistence values would typically be deleted, which raises the question of whether the extra "BCPs" and "RCPs" of ED are an important feature of this scalar field. The answer to this question relies on the domain knowledge to provide interpretation of such feature. In the light of late research which cast "serious doubts on the interpretation of BCPs in QTAIM analyses as indicating bonding interactions, however weak",[96] confronted with the recently proposed translation of these additional saddle-saddle pairs as indicators of "favorable interaction",[7] we are inclined to answer "no".

There are however more arguments standing for this answer, which do not involve the chemical interpretation of these additional critical points. For instance, the line integral convolution plots of the ED gradient demonstrated in 2 have a very similar pattern in relativistic and non-relativistic contexts, what is understandable considering that there need to be a border separating the atomic basins of Au and H (regardless of the QM model). In addition, as also exposed in 2, the Euclidean distance between the extra "RCP" and the extra "BCP" in relativistic case is so small, that their presence can be considered as a negligible (unimportant) perturbation to the topology of ED in this interatomic area.

In order to test the robustness of the presented analysis with respect to grid sampling, we have analyzed the topology and geometry of ED of m-a-g1 molecule calculated with a selected QM model (so-ZORA Hamiltonian, PBE functional and three basis sets, DZ, TZP and QZ4P) on a series of regular grids of size ranging from 1003 to 3203 points. The results of this analysis are demonstrated in 4 and they confirm prior observations: the extra saddle-saddle pairs (if present) can indeed be considered as unimportant topological features of ED, if not artificial due to the expected exponential decay of the persistence values of the extra pairs of the $\rho$ scalar field with the grid size (linear dependence of the persistence values of the corresponding pairs exhibited by the $log(\rho)$ scalar field).

What clearly discerns the relativistic ED exhibiting "BCP-RCP" pairs and nonrelativistic EDs without these pairs (for the same QM setup) is the fact that in the former case $\rho^{R}(r^{R}_{BCP})>\rho^{R}(r^{R}_{RCP})$, what causes the change of genus of the ED surface when sweeping over all values of ED, manifested by the appearance of additional critical points, as demonstrated in 2. In the latter case, the value of ED probed in the location of "relativistic BCP", $r^{R}_{BCP}$, is lower than the value of ED in the position of the "relativistic RCP", $r^{R}_{RCP}$, hence $\rho^{NR}(r^{R}_{BCP})<\rho^{NR}(r^{R}_{RCP})$. Consequently, the surface does not change its genus and the critical points do not appear. Yet, the difference of density values in these two locations ($r^{R}_{RCP}$ and $r^{R}_{BCP}$) is very small (Table 13 in ESI), which calls for a more accurate analysis involving properties which are based on the higher-order derivatives of the ED, such as the ED gradient or the ED Hessian.[32] Other QM approximations, such as the ones determining the description of electron correlation or the quality of one-electron functions do not have such a strong influence on ED that could effectively change the total order of density values in these two locations, however as the final picture is a result of interplay of many effects (e.g. counteracting relativistic corrections and electron correlation), none of them should be neglected. Notably, relativistic formalism should be used if the interacting sites involve heavy elements, as already discussed elsewhere.[20, 32, 7]

The use of accurate QM setup is essential to assert a good quality of calculated ED, however it does not change the fact that one should be cautious when drawing conclusions on bonding interactions based on the QTAIM output only. We have therefore continued with the analysis of the RDG in studied systems. As illustrated in 5, there is a striking similarity between the relativistic and non-relativistic RDGs, in geometrical (the shape of the RDG isosurfaces) and topological (the same number of connected components ($\beta_{0}$[2]) and the presence and location of minima) senses. In particular, the interaction sites between Au and H atoms are present in both contexts. This is an additional argument for the resemblance of Au-H non-covalent interactions with or without relativistic effects included.

4 Conclusions

In conclusion, we have used TDA to reinterpret the QTAIM results for gold complexes reported recently,[7] in which relativistic effects are considered to be behind the emergence of new 1-saddle - 2-saddle point pairs between non-covalently bonded Au and H atoms.

The choice of TDA is motivated by its favorable performance (by using combinatorial algorithms instead of numerical grid searches) and mathematical robustness, allowing for an explicit translation between topological features and their chemical interpretation. TDA offers the possibility of visualizing the "topological skeleton" of a scalar field, including for the EDs the same features as QTAIM (critical points and 1-separatrices between 2-saddles and maxima that correspond to "bond paths") as well as others that do not have their QTAIM correspondence (1-separatrices between 1-saddles and 2-saddles, 2-separatrices, manifolds of the Morse-Smale complex and outputs from other TDA data abstractions).

Our results demonstrate that the relationship between the topology of ED and the existence of non-covalent interactions is more complex than what can be deduced from a QTAIM picture based on the presence or absence of "BCPs" and "bond paths" between atoms. While TDA confirms that the additional saddle-saddle pairs between Au and H atoms are indeed the features of the relativistic ED which are absent for a non-relativistic ED, these feature appear not to be meaningful, due to their low persistence values and, more importantly, to the fact that these values approach zero with increased sophistication of the calculations (larger basis set and denser grids).

Furthermore, our observations are another example of the by now established fact that the ED is not the best descriptor of such weak interactions, for which other quantities such as the RDG are better suited.[46, 23, 17, 57] Our geometrical and topological analysis of the RDG, however, does not show qualitative differences between the non-relativistic and relativistic calculations, a fact which further calls into question the importance of relativistic effects to non-covalent interactions in such systems.

Taken together, these conclusions point to the need for a deeper analysis of the influence of the relativistic effects on the non-covalent interactions. To this end, we have shown here that much smaller model systems retain the same topological features as the original molecules, paving the way for more accurate calculations, such as those based on correlated wavefunction methods.

Finally, we suggest that TDA, and notably the concept of persistence–which can be considered as a process of transforming discrete invariant into a continuous one[5]–should be considered in the development of real-space analysis techniques of the chemically-relevant scalar fields.

Acknowledgements

We acknowledge support from the Polish National Science Centre (NCN) (grant number 2016/23/D/ST4/03217), the Labex CaPPA (contract ANR-11-LABX-0005-01) and I-SITE ULNE project OVERSEE (contract ANR-16-IDEX-0004), CPER CLIMIBIO (European Regional Development Fund, Hauts de France council, French Ministry of Higher Education and Research), the French national supercomputing facilities (grant DARI A0050801859), and the European Commission grant H2020-FETHPC-2017 “VESTEC” (ref. 800904). We would like to thank Dr Valérie Vallet for her comments on the draft of this paper.

Conflict of interest

There are no conflicts to declare.

Supporting Information Available:

The supporting information is available on zenodo, DOI:10.5281/zenodo.3358788.