A Study of Differential Topology on the Magnetically Induced Isotropically Averaged Lorentz Force Density of Few Simple MoleculesSupporting Information

1. Introduction

Quantum Chemical Topology (QCT), the field of research intending to study, via differential topology, scalar and vector fields which are useful to characterise molecules [1,2], has its origin in the work of Bader, who first showed that basic features of a chemical structure, $i.e.$ atoms, bonds, faces and cages, can be associated to the occurrence of critical points (CPs) of the gradient of the electron density, i.e. the points where $\nabla \rho =\mathbf{0}$ [3,4]. While the location of maxima of the density turned out unexpected (non-nuclear maxima) only in few special though interesting cases [5,6,7], the possibility to identify the most relevant atomic interactions, the "chemical bonds", via a topological analysis, i.e. by assessing whether the maxima of the electron density of two atoms are connected by a trajectory passing through the so-called Bond Critical Point (BCP), has animated an amazingly long list of contributions in literature, as in the case, to cite just one of the many discussed systems, of the H-H bonding [8,9,10,11,12,13,14,15,16], which well testifies how much the chemical community is reluctant to any shift of paradigm on bonding. In the spirit of the original Atoms In Molecules (AIM) theory [3] the topological analysis of $\nabla \rho$ is sufficient to extract a molecular graph, capturing the main elements of the molecule: the atoms and the bond paths linking a subset of atom pairs. Although the choice of density gradient as the field to unravel the chemical structure finds its physical justification in its homeomorphism with the virial field, i.e. the electronic potential energy density of an electron [17], and in the identification of density-driven bond paths as privileged exchange channels [18,19], different fields have also been used to that end [20]. We have recently studied the magnetically induced isotropically averaged Lorentz force density $〈\mathbf{F}〉$ [21], which offers a promising avenue to reveal the chemical electronic structure starting from the magnetic response.

In the present paper we are giving a thorough discussion of the topology of $〈\mathbf{F}〉$ in few simple molecules built from atoms of the first and second period, aiming at showing some of its features in the realm of QCT. The attractors of $〈\mathbf{F}〉$ occur in correspondence of both atoms and bonds, and thus their connection, as captured by trajectories passing through the saddles equivalent to the BCPs, can be expected to be more complicated than that of $\nabla \rho$, but likely similar to that obtained from the Laplacian of the density. The latter has been used to introduce $L-$graph [22], and, by analogy, we will here introduce $〈\mathbf{F}〉-$graphs. The analysis will be enriched with a discussion on the count of critical points, whose topological indices should be bound by relations of the kind of the Poincaré-Hopf theorem.

The paper is organised as follows: after some basic information on the Lorentz force field (Section 2.1), an introduction to the needed items of differential topology (Section 2.2), and some computational details (Section 3), we will first present in Section 4 a systematic report of the topological features of $〈\mathbf{F}〉$ computed in prototypes of acyclic molecules (lithium hydride, methane, ethane, ethylene, acetylene, propyne and 3-borapropyne), aromatic, antiaromatic and non-aromatic molecules (benzene, planarized cyclooctateraene and borazine, respectively) and a cage molecule (cubane), and we will then discuss and resume the main finding in Section 5.

2. Theory

2.1. The Magnetically Induced Isotropically Averaged Lorentz Force Density

At second order in the external magnetic field $\mathbf{B}$, the energy acquired by a freely tumbling isolated closed shell molecule is [23,24]

$$W=-{\displaystyle \frac{1}{2}}{\xi }_{\mathrm{i}so}{B}^2$$

(1)

where the isotropically averaged magnetizability is obtained as one third of the trace of the magnetizability tensor, ${\xi }_{\mathrm{i}so}={\displaystyle \frac{1}{3}}({\xi }_{xx}+{\xi }_{yy}+{\xi }_{zz})$, or more briefly, using here and in the following the Einstein convention of summing over repeated indices, ${\xi }_{\mathrm{i}so}={\displaystyle \frac{1}{3}}{\xi }_{\alpha \alpha }$.

Magnetizability can be computed from the magnetizability density [25,26],

$${\xi }_{\mathrm{iso}}=\int {\Xi }_{\mathrm{iso}}{d}^{3}r,$$

(2)

where [21]

$${\Xi }_{\mathrm{iso}}={\displaystyle \frac{{ϵ}_{\alpha \beta \gamma }}{6}}(r_{\beta }-r_{0\beta }){\mathcal{J}}_{\gamma \alpha },$$

with ${ϵ}_{\alpha \beta \gamma }$ the Levi-Civita symbol, $r_{\beta }$ a Cartesian displacement from an origin with Cartesian coordinates $r_{0\beta }$, and $\mathcal{J}$ is the $3\times 3$ current density tensor, which allows the computation of the magnetically induced first order current density via ${\mathbf{J}}^{\mathbf{B}}=\mathcal{J}\mathbf{B}$. The magnetizability density depends on the arbitrary origin of the coordinates ${\mathbf{r}}_0$, so that it is not considered a useful tool for QCT. However, by manipulation of the magnetically induced energy [21],

$$\begin{array}{ccc}\hfill W& =& -{\displaystyle \frac{B^2}{2}}\int {\Xi }_{\mathrm{i}so}{d}^{3}r=-{\displaystyle \frac{B^2}{12}}\int {ϵ}_{\alpha \beta \gamma }(r_{\beta }-r_{0\beta }){\mathcal{J}}_{\gamma \alpha }{d}^{3}r\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& =& -\int {〈F〉}_{\beta }(r_{\beta }-r_{0\beta })d^{3}r,\hfill \end{array}$$

(4)

it can be realised that the isotropically averaged Lorentz force density

$$\begin{array}{ccc}\hfill {〈F〉}_{\beta }& =& {\displaystyle \frac{B^2}{12}}{ϵ}_{\beta \gamma \alpha }{\mathcal{J}}_{\gamma \alpha },\hfill \end{array}$$

(5)

is an origin-independent vector (and not tensor) field amenable to QCT studies, like $\nabla \rho$. A positive or negative value of the Divergence of the Isotropically Averaged Lorentz force density (DIAL) indicates whether the local contribution to the magnetizability is positive or negative. The usefulness for the study of magnetic aromaticity of $\pi$-DIAL (the contributions of $\pi$ electrons to DIAL) has been recently discussed [27].

Apart from a factor, DIAL, the divergence of 5, has been independently introduced in quantum chemistry by Barquera-Lozada, who called it $\mathrm{tp}J\left(\mathbf{r}\right)$ [28].

2.2. Topological Analysis of the Critical Points

The critical points of a 3-dimensional vector field v can be classified according to their rank r and signature s as $(r,s)$[3], where r is the number of non-null eigenvalues of the Jacobian matrix at the critical point, and $s=p-n$ is the difference between the number p of positive and the number n of negative real components of the eigenvalues of the Jacobian. The sign of the determinant of the Jacobian matrix at the critical point $x^{*}$ is known as its index: ${\mathrm{ind}}_{{\mathrm{x}}^{*}}\left(v\right)=\mathrm{sign}\det J[v\left(x^{*}\right]={(-1)}^n$[29]. In appropriate cases the sum of the indices of the critical points of a vector field must follow the Poincaré-Hopf theorem, which has been considered a check of consistency of the topological analysis [30].

Before discussing this constraint for our field $〈\mathbf{F}〉$, we find it useful to introduce the Euler characteristic $\chi \left(K\right)$ of a CW-complex K, which is a collection of cells[31]. In our case it will be sufficient to consider 4 kinds of cells: 0-cells (isolated points in space), 1-cells (deformable lines connecting two 0-cells), 2-cells (the deformable 2-dimensional space contained within a closed loop of 1-cells) and 3-cells (the deformable 3-dimensional space contained between 2 or more 2-cells). The Euler characteristic of the CW-complex $\chi \left(K\right)$ can be obtained as the difference between even-dimensional cells and odd-dimensional cells:

$$\chi \left(K\right)=\sum _k{(-1)}^k{c}_k=c_0-c_1+c_2-c_3,$$

(6)

where $c_k$ is the number of k-cells. If the CW-complex is compact, it can be contracted to a single point and $\chi \left(K\right)=1$. Contraction can be equally performed on a subset of cells of the CW-complex, while leaving the Euler characteristic unchanged. A sketch of contraction of 2 CW-complex, each of them compact, is given in Figure 1.

Eq. 6, which, for a single compact CW-complex, can be read as $V-E+F-C=1$, where V, E, F and C denote the number of vertices, edges, faces and cages, is consistent with the classical Euler equation $V-E+F=2$, proposed long ago for a convex polyhedron.

Coming back to the critical points of a 3-dimensional vector field, we will be interested only in full-rank CPs (those with $r=3$, and thus no null eigenvalue of the Jacobian) and Table 1 resumes their indices.

The four possible full-rank CPs, when ordered for decreasing number of negative eigenvalues, are $(3,-3)$, $(3,-1)$, $(3,1)$ and $(3,3)$, and will occur in number $n_3$, $n_2$, $n_1$, $n_0$, respectively. In that order, the indices have alternating sign, so that, if the one-to-one correspondence given in column 5 of Table 1 occurs between the critical points of a given kind and the cells of a compact CW-complex ($n_k=c_k$ for each k), Eq.6 can be rewritten as

$$\sum _{k=0}^3{\mathrm{ind}}_k\left(v\right)n_k=1,$$

(7)

which matches the 3-dimensional Poincaré-Hopf theorem,

$$\sum _{k=0}^3{\mathrm{ind}}_k\left(v\right)n_k=\chi \left(M\right),$$

(8)

where $\chi \left(M\right)=1$ is the Euler characteristic of the compact 3-dimensional manifold where the vector field is computed. The Euler characteristic $\chi \left(M\right)$ is an invariant in differential topology, so that it only depends on the manifold M and not on the vector field. To understand how this statement can be compatible with a non-zero value of $\chi \left(M\right)$, it is important to recall that, if the manifold M has a boundary $\partial M$, the theorem holds provided that the vector field is everywhere outward oriented on the boundary [32]. For a vector field v that is everywhere inward-oriented on the boundary, the outward orientation can be recovered simply changing sign to the field, $v^{\prime }=-v$, which implies a change of sign of the signature, of all indices, and therefore, for the one-to-one correspondence given in the last column of Table 1 ($n_{3-k}=c_k$ for each k),

$$\sum _{j=0}^3{\mathrm{ind}}_j\left(v\right)n_j=\sum _{j=0}^3{\mathrm{ind}}_{3-j}\left(v\right)n_{3-j}=-\sum _{k=0}^3{\mathrm{ind}}_k\left(v^{\prime }\right)c_k=-1,$$

(9)

which, considering that the Euler characteristic of a spherical surface is $\chi (\partial M)=2$, can be recognised as a special case of the equation for an inward-oriented field [33]

$$\sum _{k=0}^3{\mathrm{ind}}_k\left(v\right)n_k=\chi \left(M\right)-\chi (\partial M).$$

(10)

For a 3-manifold, $\chi \left(M\right)=1$. However, if one intends considering the real coordinate space including points at infinity, one must consider the real projective space $\mathbb{R}{\mathbb{P}}^3$, which can be thougth as a closed 3-dimensional ball with antipodal points identified, and has $\chi \left(M\right)=0$. To get a zero in Eq. 8, it is necessary to sum also the asymptotic CPs: the $n_{3a}$ asymptotic attractors (with index ${\mathrm{ind}}_3=-1$, which are maxima if the vector field is a gradient) and the $n_{0a}$ asymptotic repellers (with index ${\mathrm{ind}}_0=1$, which are minima if the vector field is a gradient). Consistency with the alternative description of a 3-manifold with a boundary $\partial M$, shows that if the field is everywhere inward/outward oriented on the boundary there is a single asymptotic repeller/attractor. In the general case, when the orientation of the vector field is changing over a boundary enclosing all isolated CPs, the sum of indices of the isolated CPs, $\Sigma {\mathrm{ind}}_k{n}_k$ can differ from $\chi \left(M\right)$. It has been proposed, in a study of the molecular electrostatic potential (MEP), to deduce the number of asymptotic CPs by plotting the MEP on a spherical surface large enough to contain all isolated critical points, and counting the numbers of closed island-like regions of negative and positive values [30]. In our case, the vector field generally has $(3,-3)$ CPs on the outer part of the molecule (an inward-directed field) and thus will have $\Sigma {\mathrm{ind}}_k{n}_k=-1$. This should not hold when $(3,+3)$ CPs occur on the outer part of the molecule.

Eq. 8, the celebrated Poincaré-Hopf theorem assumes that CPs are isolated, although generalization are possible [33]. In acetylene, we will find loops of degenerate $(2,s)$ CPs. It has been proposed that similar cases should be dealt considering a symmetry-breaking and then considering the resulting isolated CPs [30]. Indeed, one feature of the indices is their stability, which means that small perturbation do not change the sum of the indices of close critical points. A loop of $(2,s)$ CPs can always be thought as a loop of alternating $(2,s-1)$ and $(2,s+1)$ CPs having a null sum of indices.

3. Computational Methods

Optimised geometries and magnetically perturbed wavefunctions have been obtained using Gaussian 16 [34] and the BHandHLYP functional, recently assessed to perform well in calculation of magnetic properties [35], combined with the pcSseg-4 basis set [36]. Essential parameters of the optimised geometries are given in Table S1.

Critical points have been determined by the Newton-Raphson algorithm as implemented in SYSMOIC [37], starting from nuclei and then from points sampled out from a parallelepiped containing the molecule, with a step of 0.1 au, reduced to 0.05 au, for molecules with triple bonds, and checking that further reduction did did not increase thje number of CPs.

The CPs found are sometimes in disagreement with a previous calculation carried out with the smaller 6-31G(d) basis set [21]. We have checked that number and kind of critical points is stable upon basis set enlargement, by performing a calculation at the aug-pcSseg4 level for benzene, as well as adding additional basis functions on points that were missing in lower level calculations (Table S2). Calculation with additional basis set functions on the $(3,+3)$ CPs above faces in benzene, borazine, cyclooctatetraene and cubane gave the same unaltered set of CPs. Considering the small size of the molecule studied, all calculations (both geometry optimisation and magnetic perturbation) have been carried out at the BHandHLYP/pcSseg-4 level.

The critical points are displayed according to the following convention: $(3,-3)$ CPs are small red spheres, (3,+3) CPs are small blue spheres, $(3,\pm 1)$ CPs (with distinguished eigenvalues) are displayed as three crossing segments going along the eigenvectors of the Jacobian, blue for positive eigenvalues, red for negative eigenvalues. In case of degenerate eigenvalues a circle is shown in the plane of the eigenvectors corresponding to the degenerate eigenvalues, with the same colour code convention: a blue circle for positive degenerate eigenvalues, a red circle for negative degenerate eigenvalues. Once the CPs have been found, we have computed trajectories moving out of the $(3,-1)$ CPs along the single positive eigenvalue, to define what we call the $〈\mathbf{F}〉$-graph. $〈\mathbf{F}〉$-graphs are fairly rich 3-dimensional objects; they can be best appreciated moving them interactively, which is possible operating on the 3d files generated for all of them, using the the v3d code, freely distributed on the SYSMOIC website.

4. Results

We are going now to report the location and number of the CPs found for the molecules chosen. A summary of the results is given in Table 2. Lists of location and eigenvalues of the Jacobian at the CPs are given in the Supporting Information (Table S3).

4.1. Lithium Hydride

As a first example we show the $〈\mathbf{F}〉$ field for LiH. In this simple molecule the field has two attractors separated by a $(3,-1)$ CP. The field is everywhere inward-oriented (on a boundary containing all isolated CPs) and the sum of indices is consistently equal to $-1$.

4.2. Methane

Figure 3 shows the 29 critical points of $〈\mathbf{F}〉$ for methane (their list is given in the Supporting Information). There are 5 $(3,-3)$ attractors, 1 on the C atom and 4 on the C–H $\sigma$ bonds. If connected by trajectories passing through the $(3,-1)$ CPs, they form a tetrahedral shaped complete graph $K_5$: every attractor is bound to any other one, and indeed the number of $(3,-1)$ CPs is $n_2=\left({\displaystyle \genfrac{}{}{0pt}{}{5}{2}}\right)=10$. The number of topologically triangular faces is also $n_1=\left({\displaystyle \genfrac{}{}{0pt}{}{5}{3}}\right)=10$ and a $(3,+1)$ face critical point can be associated to each of them. Eventually, the $n_0=4$$(3,+3)$ CPs can be associated to tetrahedral cages formed by three H atoms and the central C atom. The sum of indices is $-1$ as expected for an inward-oriented vector field, with a single repeller at infinity.

The complete graph $K_5$ is one of the two basic $gt$-nonplanar graphs, those that cannot be drawn on a sheet of paper without crossing of edges. $gt$-nonplanar graphs are unusual topologies in organic chemistry [38,39], when one considers standard bonding schemes, often compatible with the molecular graphs obtained from $\nabla \rho$. However, we see that even the simplest organic molecule, namely methane, shows a $gt$-nonplanar graph, when studied with the $〈\mathbf{F}〉$ field.

4.3. Ethane

Figure 3 shows the 57 critical points of $〈\mathbf{F}〉$ and the trajectories defining the $〈\mathbf{F}〉-$graph for staggered ethane. All topological features are preserved in eclipsed ethane (see Supporting Information, Figure S1). The number of CPs is just one less than the sum of the CPs of two methane molecules: indeed one of the $(3,-3)$ CPs is shared between the two methyl moieties. The methyl unit has thus a conservative pattern of CPs, much as reported for the Laplacian of the density [22]. The sum of indices equals -1, as expected for an inward-oriented vector field.

4.4. Ethylene

The $〈\mathbf{F}〉$ field for ethylene has 31 isolated CPs. 8 CPs are attractors of $〈\mathbf{F}〉$: 2 on the C atoms, 4 on the C–H $\sigma$ bonds and a pair above and below the molecular plane in correspondence with the double bond. As a consequence, in the middle of the C–C bond, rather than finding an attractor, we find a $(3,-1)$ CP. In addition to these CPs, on each methylene moiety we find 5 $(3,-1)$, 4 $(3,+1)$ and 2 $(3,+3)$ CPs. The $(3,-1)$ CPs are half of those found in the methyl moiety in ethane: The CH attractors are linked to the $\pi$-bond attractors from above and below (4 paths for each methylene); the $(3,-1)$ CP close to the C attractor, in the HCH region, has paths going on either of the close CH attractors, separated by a single path that goes directly in the $(3,+1)$ CP. The (3,+3) CPs appear between two topological triangles formed by the two $\pi -$bond attractors and either the CH-bond attractor or a $(3,+1)$ CP located close to the C atom. The sum of indices is $-1$ as expected for an inward-oriented vector field.

4.5. Acetylene

The $〈\mathbf{F}〉$ field of acetylene is not formed by isolated CPs only, but three circles of critical points appear: a circle of $(2,-2)$ CPs in the middle of the molecule and two circles of (2,0) CPs close to the location of C atoms. Along the molecular axis 9 CPs are found: 4 $(3,-3)$ CPs located on C atoms and on C–H bonds, a $(3,+1)$ CP in the middle of the C–C bond, two (3,+3) CPs between the middle of the C–C bond and the C atoms, and two $(3,-1)$ CPs between the C atoms and the attractors on C–H bonds. The latter $(3,-1)$ CPs are only 0.097 au far apart from the $(3,-3)$ CPs located on C atoms. Disregarding the continuous circles of CPs, the sum of indices is $-1$ as in all previous molecules. It is to be noted that, as far as the change of the vector field is considered only along the molecular axis, the sequence of the 3 CPs $(3,+3)-(3,+1)-(3,+3)$ is equivalent to a sequence $(3,-1)-(3,-3)-(3,-1)$ which could have been expected for a $\sigma$ bond.

4.6. Propyne

The $〈\mathbf{F}〉$ of propyne has 55 CPs. The insertion of the methyl has perturbed each of the three circles of rank-2 CPs of acetylene, producing two sets of three alternating CPs. In the case of the central circle of $(2,-2)$ CPs, 3 $(3,-3)$ points are generated alternating with $(3,-1)$ CPs. In the case of the C-centred circles of $(2,0)$ CPs, 3 $(3,-1)$ CPs are formed alternating with 3 $(3,+1)$ CPs. The sum of indices along each circle is zero. Excluding these $3\times 6=18$ CPs, we have 37 CPs, just one less than the sum of the isolated CPs of acetylene and methyl: 9+29=38. Indeed, location and kind of CPs are very similar to those of the smaller molecules, and the $(3,-3)$ on the C-Me bond is shared by the two moieties.

4.7. 3-Borapropyne

The $〈\mathbf{F}〉$ field of 3-borapropyne has 43 isolated CPs. 7 out of the 9 $(3,-3)$ CPs occur close to atoms or to the centres of $\sigma$ bonds. Two further attractors occur in the plane bisecting the HBH angle, in place of the circle of $(2,-2)$ CPs found in acetylene. The $(3,+3)$ repellers occur 3 around the B atom, 2 around the central C atom, and 1 close to the peripheral C atoms. Of the 15 $(3,-1)$ CPs, 12 are out-of-plane (the 4 connecting the attractors on BH bonds and the B atom have 2 complex eigenvalues with negative real component) and the remaining 3 occur on the symmetry axis. No trajectory passing through a $(3,-1)$ CP links directly the BH attractors. Of the 13 $(3,+1)$ CPs, 7 can be associated to loops of alternating $(3,-3)$ and $(3,-1)$ CPs, 3 alternate with repellers close to the B atom, 2 alternate with repellers close to the central C atom, and 1 acts as a repeller of trajectories (not shown) originating from the $(3,+3)$ CP close to the B atom along the bisector of the HBH angle.

4.8. Benzene

The $〈\mathbf{F}〉$ field of benzene has 99 CPs. The 18 $(3,-3)$ CPs are located on C atoms, CH and CC bonds. Of the 14 $(3,+3)$ CPs, 12 are close to the C atoms, located towards the inside and the outside of the ring, and 2 are above and below the ring. The 36 $(3,-1)$ CPs come from 6 sextets around the C atom. Each sextet, in turn, is composed by two triplets of CPs placed symmetrically above and below the ring. The 3 symmetry-unique $(3,-1)$ CPs all occur in the plane of symmetry containing two opposing CH bonds, one very close to the C atom, one at some distance from the centre of the CH bond and one inside the ring, at an height of 0.952 bohr, close to the height traditionally chosen to display maps of the $\pi$ current density [37,40]. Notably, along the CH bond, we do not find a $(3,-1)$ CP, but the two attractors are separated by a $(3,+3)$ CP. The 4 symmetry-unique CPs out of the 31 $(3,+1)$ CPs occur one in the middle of the ring, one in the CH bond direction inside the ring, and two flanking the CH bond on either side. The sum of indices is +1. This can be understood in terms of presence of a repeller and two attractors at infinity, according to Leboeuf et al. Indeed the map of the flux of $〈\mathbf{F}〉$ out of a spherical surface shows two zones of positive values, compatible with the presence of two asymptotic attractors (see Supporting Information, Figure S2).

4.9. Cyclooctatetraene

The $〈\mathbf{F}〉$ field of planarised cyclooctatetraene has 123 CPs, with 28 $(3,-3)$ CPs, occurring on atoms on CH and long CC bonds, and in pairs above and below the middle of short CC bonds. The trajectories following the eigenvector corresponding to the positive eigenvalue of the 44 $(3,-1)$ CPs connect the C atom attractor with the CC long bond attractor, the two $\pi$-bond attractors, the CH-bond attractor with the adjacent CC bond and $\pi$-bond attractors. The trajectories passing through the CH-bond attractor form topological triangles or squares each associated with a $(3,+1)$ CP. Of the remaining 17 $(3,+1)$ CPs, 16 occur in pairs close to the C atom, alternating with an equal number of $(3,+3)$ CPs (much alike the case of benzene, although with less symmetry), and 1 in the ring centre separating two $(3,+3)$ CPs located above and below the ring. The sum of indices is +1. This can be understood in terms of presence of a repeller and two attractors at infinity, according to Leboeuf et al. Indeed the map of the flux of $〈\mathbf{F}〉$ shows two zones of positive values, compatible with the presence of two asymptotic attractors (see Supporting Information, Figure S2).

4.10. Borazine

The $〈\mathbf{F}〉$ field of borazine has 129 isolated CPs. The 24 attractors are located on N and B atoms, on BH and NH bonds (closer to the more electronegative atom), and in pairs, above and below the BN bonds, shifted towards the more electronegative N atom. Around B and N atoms there are 2 and 3 $(3,+3)$ repellers, respectively, and other 2 repellers occur above and below the ring as in benzene and cyclooctatetraene. The 48 $(3,-1)$ CPs are 1 for each of the 6 close pairs of $\pi$ CN bonds attractors, 8 around each N atom and 6 around each B atom. Both atoms have 2 $(3,-1)$ CPs above and below them, shifted towards the inside of the ring in case of B atoms; through any of these CPs a trajectory passes linking two up or down pairs of $\pi$ attractors. The pair of $\pi$-linking trajectories are slightly on the outside of the ring for N atoms, and in the inside for B atoms. The inside location of the B trajectories is reminiscent of those of benzene. A second difference between B and N atoms is that, for the B atoms, 4 CPs occur on the sides of the BH bond and hold 4 trajectories linking the BH attractor to the $\pi$ attractors, while for the N atoms 2 CPs occur in the vertical plane containing the NH bond; from these CPs trajectories are also found towards the $\pi$ maxima, but they do so after being scattered away by the $(3,-1)$ CP located above the N atom. The latter $(3,-1)$ CP works as a scatterer towards the $\pi$ attractors also for the trajectories starting from the C and N nuclear attractors. The 40 $(3,+1)$ CPs occur 1 at the ring centre, 2 and 3 alternating with repellers close to N and B atoms, respectively, and then 3 and 5 around the trajectories of the $〈\mathbf{F}〉$ field clustering around B and N atoms and forming topologically square pyramids, one squeezed (for N atoms).

4.11. Cubane

The $〈\mathbf{F}〉$ field of cubane has 137 CPs. The 28 $(3,-3)$ attractors correspond to the 8 C nuclei and the 20 $\sigma$ bonds. Notably the $(3,-3)$ CPs of C–C bonds are shifted towards the outer part of the cube, much as happens for the bond paths obtained from $\nabla \rho$[41]. 23 $(3,+3)$ CPs are found: 6 above the faces, 1 in the middle of the cube and, 2 for each C atom, along the diagonal of the cube, one towards the exterior and one towards the interior. The 48 $(3,-1)$ CPs are 8 symmetry-equivalent sextets located around C atoms. The distortion of the valence angles from the tetrahedral values is associated with reduction of the number of CPs, if one compares them with the 10 $(3,-1)$ CPs of methane. The 5 attractors close to a C atom, if considered as graph vertices, do not form any more a complete graph using $(3,-1)$ CPs as edges: the central attractor on the C atom is connected via $(3,-1)$ CPs to only the three attractors on the C-C bonds but it is connected to the $(3,-1)$ CH-attractor via a $(3,+3)$ CP. The CH-attractor is connected to the CC-attractors, but there are no $(3,+1)$ CPs on the topological triangle formed by two CC-attractors and the CH-attractor. Each of the $4=\left({\displaystyle \genfrac{}{}{0pt}{}{4}{3}}\right)$ topological triangles formed of by the C atom and the 3 CC-attractors is associated with a $(3,+1)$ CP. 6 More $(3,+1)$ CPs occur on the faces of the cube; they can be thought as the result of the disappearance of the $(3,+1)$ CPs on the topological triangles CH-CC-CC upon joining the 8 CH units forming the molecule. The sum of indices is +5, which can be understood, in terms of the algorithm proposed by Leboeuf in terms of 6 positive zones inside an overall negative zone on the sphere (Figure 9).

5. Discussion

The molecules studied in this paper give some insight on the capability of the isotropically averaged Lorentz force density $〈\mathbf{F}〉$ to contribute in the study of the molecular electronic structure.

As a first point, we note that the attractors of $〈\mathbf{F}〉$ tend to be greater in number than the attractors of $\nabla \rho$, where their number corresponds to the number of nuclei (apart from the seldom observed non-nuclear maxima). In the case of $〈\mathbf{F}〉$, the attractors correspond to both nuclei and bonds. In effect, attractors are found very close to H nuclei, only for the ionic Li–H bond, while for B–H, C–H and N–H bonds the attractors are somewhere along the bonds, closer to the more electronegative atom. In all cases where the electronic structure can be described by core $1s$ electrons pairs and $\sigma$ bonds, the total number of attractors $n_3$ equals the number of electron pairs (see Table 1). This holds true also for ethylene and cyclooctatetraene, but not for benzene and borazine, showing less (18 vs 21) and more (24 vs 21) CPs than the stoichiometric electron pairs. The case of the triple bond also deserves attention. In acetylene the number of attractors $n_3$ is 3 less than the number of stoichiometric electron pairs, which hints at the missing triple bond. Then, in propyne, $n_3$ matches the number of stoichiometric electron pairs, but in 3-borapropyne $n_3$ is 1 less than the number of stoichiometric electron pairs. To summarise, although a one-to-one correspondence between electron pairs and attractors holds for saturated molecules, it fails in most cases where bonding can be expected to develop far off of the lines joining the nuclei.

We find most interesting that C–C bonds of different order are not characterised by the same $(3,-3)$ CP as for $\nabla \rho$: only for $\sigma$ bonds, we find a single $(3,-3)$ attractor, while for the double bond of ethylene we find a $(3,-1)$ saddle flanked by two close attractors, for the triple bond of acetylene we find a $(3,+1)$ saddle surrounded by a circle of $(2,-2)$ CPs. Notably, in propyne, the circle of $(2,-2)$ CPs splits in 3 attractors and 3 saddles, while the $(3,+1)$ CP is preserved. The appearance of three $(3,-3)$ CPs is a feature of the symmetry of the molecule, as exemplified by the $adhoc$ studied $C_{2\mathrm{v}}$BH₂−CCH molecule, where only two $(3,-3)$ CPs are found. However, the $(3,+1)$ CPs in the middle of the bond remains, indicating that it is a characteristic feature of the triple bond. The correspondence between bond orders and number and kind of attractors seems preserved to some extent. In cyclooctatetraene, long and short bonds are associated with 1 and 2 attractors, respectively. For the partial bond orders of benzene we found a single attractor. We note that the possibility to reveal different orders of chemical bonds has been very early noted in QCT, in studies of the electron localisation function [42] and the Laplacian of the density [22]. To our knowledge, this is the first time that this possibility is highlighted in terms of the magnetic response.

As a second point, when one considers the critical points close to an atom, a certain degree of transferability can be observed, as previously observed for the Laplacian [22]. Indeed, location and kinds of CPs in ethane and the methyl group of propyne are almost coincident with those of methane. Some transferability is also observed among the critical point found close to the CH₂ group of ethylene and the CH groups of benzene and cyclooctatetraene. Differences are observed instead in the trajectories connecting the attractors on the ring bonds: in benzene these trajectories are shifted towards the inside of the ring, while in cyclooctatetraene they are shifted towards the outside. Notably, in borazine these trajectories are shifted alternatively inside and outside. This hybrid behaviour is well compatible with the main characterisation of borazine as a non-aromatic molecule [43], although the different dimensions of the domains of N-centred and B-centred trajectories could well be compatible with some weak aromaticity, which, using ring current strength was proposed several years ago [44] and recently rediscovered, and enriched with considerations based on excited state currents [45].

As a last point, we notice that for the 4 molecules studied, endowed with a ring structure, the $〈\mathbf{F}〉$ field turns out not everywhere inward-oriented like $\nabla \rho$, but lines of positive (outward) flux of the force have been found above and below all rings. In this case, the straightforward expectation of a sum of topological indices equal to -1 is not valid. An adaptation of the algorithm proposed by Leboeuf [30] of counting islands of positive and negative flux of the vector over a sphere surrounding the molecule allows in these 4 cases to justify the sum of indices. We find it noteworthy that closed-shell diamagnetic molecules can have at large distance narrow spherical angles where the local magnetic response is paramagnetic.

We are confident that further studies on the $〈\mathbf{F}〉$ field will be useful in the study of molecular electronic structure from the QCT perspective.