Analysing Magnetically Induced Currents In Molecular Systems Using Current-Density-Functional Theory

A suite of tools for the analysis of magnetically induced currents is introduced. These are applicable to both the weak-ﬁeld regime, well described by linear response perturbation theory, and to the high-ﬁeld regime, which is inaccessible to such methods. A disc-based quadrature scheme is proposed for the analysis of magnetically induced current susceptibilities, providing quadratures that are consistently deﬁned between diﬀerent molecular systems and applicable to both planar 2D and general 3D molecular systems in a black-box manner. The applicability of the approach is demonstrated for a range of planar ring systems, the ground and excited states of the benzene molecule and the ring, bowl and cage isomers of the C 20 molecule in the presence of a weak magnetic ﬁeld. In the presence of a strong magnetic ﬁeld, the para- to dia-magnetic transition of the BH molecule is studied, demonstrating that magnetically induced currents present a visual interpretation of this phenomenon, providing insight beyond that accessible using linear-response methods.

In the present work, we utilise LAOs, as done by Jusélius, Gauss and Sundholm 21,22 to determine magnetically induced current susceptibilities in molecular systems. The physical current is routinely calculated in the application of non-perturbative current-density-functional approaches [23][24][25] for molecules in magnetic fields. The magnetically induced currents can then be easily computed by finite differences from the physical current, evaluated for a small perturbative magnetic field applied along one Cartesian axis. In particular, we use a recently constructed family of current-dependent meta-generalized-gradient approximation (mGGA) functionals to determine magnetically induced current susceptibilities. The use of such a non-perturbative approach also allows for the determination of current densities directly as a function of magnetic field strength; this aspect is explored for the BH molecule in this work.
We commence by outlining the theoretical background underpinning the non-perturbative calculations used in the present study, which employ LAOs as the basis functions in which to expand the molecular orbitals. We then briefly review a recent implementation of currentdensity-functional theory (CDFT) in a magnetic field, highlighting the ease with which the physical current density j induced by the field B and the related current density susceptibility tensor J Bτ ν (r) = ∂j ν (r)/∂B τ (where ν and τ label Cartesian components) may be extracted from CDFT calculations.
The physical current density is a rich source of chemical information; its topology reflects the chemical structure of the molecule and the interaction of the electronic structure with an externally applied magnetic field. Magnetically induced current susceptibilities have long been used in ring-current models to provide insight into nuclear-magnetic-resonance (NMR) chemical shifts 1,2,[26][27][28][29] and as a criterion for assessing the aromaticity of molecules. 30 They give insight into electron delocalisation [31][32][33][34] and their analysis has also been applied to probe hydrogen bond strengths. [35][36][37] Nonetheless, the physical induced current j is a complicated vectorial quantity associated with a particular orientation of the applied magnetic field B with respect to the molecular frame. The induced current density susceptibility J Bτ ν is a tensorial quantity, reflecting the vectorial nature of both the applied magnetic field and the induced physical current density. As such, the analysis of these quantities is less straightforward than that for simple scalar quantities.
In general, there are two main approaches that may be used for the analysis of currents induced in molecular systems by external magnetic fields. The first of these are integration techniques which, by constructing numerical quadratures over two-dimensional planes, allow the current density to be probed in specific parts of a molecule. 21,22 Secondly, topological techniques employing concepts from vector-field analysis such as separatrices and stagnation graphs are used to analyse the induced current fields. 8,[38][39][40] Both approaches can provide quantitative information on the nature of the electron delocalisation in chemical species and their interactions with external fields. In the present work, we focus on the approach of integrating current densities and propose a simple disc-based quadrature method, with which a measure of current flow along a chemical bond may be obtained. We show that this approach complements those designed to calculate ring currents and gives similar insight.
In addition, our approach can readily be applied to general three-dimensional molecular structures, owing to the finite spatial extent of the quadrature used. The advantages of this flexibility are discussed and the application of this method to non-planar molecules demonstrated.
Finally, we show how the non-perturbative nature of our implementation allows us to explore how the physical induced current changes as a function of the external magnetic field. The nature of the induced currents can be used to help to rationalise the changes in the energy of molecular systems in strong magnetic fields of strength up to one atomic unit

Theory
In the present work, systems are considered in the presence of a uniform magnetic field represented by the vector potential, where R G is an arbitrary gauge-origin. The effects of the magnetic field are treated nonperturbatively, employing a basis of LAOs, 20 consisting of a standard Gaussian function χ µ (r) centred at R K multiplied by a fielddependent complex phase factor, yielding properties that are gauge-origin invariant and which are correct to first order in the magnetic field, resulting in more rapid convergence towards the basis set limit than a standard Gaussian basis.
The London code 41 was the first to provide a range of functionality for electronic structure calculations in strong magnetic fields, including Hartree-Fock (HF), 42 current-densityfunctional, 23 configuration-interaction, 43 Møller-Plesset perturbation theory, coupled-cluster 44 and complete-active-space (CAS) self-consistent field (SCF) theories. 43 This has been followed by several other electronic structure packages including BAGEL, 45,46 ChronusQ 47 and our own code QUEST, 48 which also provide functionality for calculations using LAOs. In the present work, we use the implementations of HF and CDFT in the QUEST program to determine magnetically induced physical currents and current susceptibilities.

Current-density-functional theory
In the Vignale-Rasolt formulation of CDFT, 49-51 the Kohn-Sham (KS) equations take the form 1 2 where p is the canonical momentum operator and s is the spin operator. In CDFT a noninteracting system is introduced to reproduce both the charge density, where i labels occupied orbitals with spin σ, and paramagnetic current density, of the physical system. The KS potentials (u s , A s ) are where (v ext , A ext ) are the physical external potentials, v J , is the Coulomb potential and the exchange-correlation scalar and vector potentials are A central challenge for CDFT calculations is to define an exchange-correlation functional This leads to a well-defined and properly bounded iso-orbitial indicator when applied to the Tao-Perdew-Staroverov-Scuseria (TPSS) functional (see, for example, Ref. 58 for comparisons) and the resulting form has been called cTPSS. In the present work, we use the same modification for non-perturbative calculations in the presence of external magnetic fields.

Results
Determining Magnetically Induced Currents and Current Susceptibilities at the CDFT level The GIMIC program 21,22 determines magnetically induced current susceptibilities by interfacing to programs that can perform calculations for magnetic response properties-for example, the determination of NMR shielding constants. This approach has the advantage that a wide variety of programs can be used, along with a wide range of electronic structure methods. However, for CDFT functionals, few linear response implementations exist. In the present work, the need for linear response calculations is avoided since we can access these quantities as a direct by-product of our non-perturbative calculations.

The physical magnetically induced current
Our implementation in the QUEST program 48 allows for the calculation of the LAO oneparticle density matrix (1-RDM) at the HF and CDFT levels in the presence of static, uniform magnetic fields of arbitrary strength. Once the density matrix has been determined, the charge density at a grid point can be calculated as and the paramagnetic current density as The physical current density is then constructed as The induced-current susceptibility tensor can be evaluated from three SCF calculations, in which the external magnetic field is aligned along each of the three orthogonal Cartesian axes, respectively. With a weak field along each Cartesian direction τ , the partial derivative can be evaluated numerically, yielding the current density susceptibility tensor.

The Biot-Savart Law
Once the induced-current susceptibility tensor is determined, it is possible to determine the NMR shielding tensor for a nucleus K with associated magnetic moment M K . To lowest order, the energy of a nuclear magnetic moment M K with vector potential A K (r) = µ 0 4π M K ×r K r 3 K in an electronic system with current density j(r) is given by A K (r) · j(r)dr, yielding the following expression for the shielding constant: 59 Here the Einstein summation convention is used and ε αβγ is the Levi-Civita tensor. The isotropic shielding constant is then given as σ K iso = 1 3 Trσ. We note that the Biot-Savart law has long been used in the context of ring-current calculations. [60][61][62][63] The evaluation of NMR shielding tensors has been implemented in the QUEST program. The integrals required in Eq. (14) have been implemented both numerically, using standard DFT quadrature, and analytically for LAOs. 64 In both cases, NMR shielding tensors can be calculated either for specified nuclei or for all nuclei in a molecule. It may be possible to further optimize the numerical evaluation of the integrals in Eq. (14) by exploiting their locality to each nucleus K, using tailored quadratures. However, we do not explore that possibility here and all results in the present work use the analytic implementation.

Functional Dependence
The Biot-Savart Law in Eq. (14) provides a direct way to assess the quality of magnetically induced currents in the vicinity of nuclei. The calculation of NMR shielding constants can be challenging for density-functional methods, particularly for nuclei such as 15 N, 17 O, 19 F. It was shown in Ref. 24 that cTPSS gives modest improvements over GGA level functionals such as the Perdew-Burke-Ernzerhof (PBE) functional for the evaluation of NMR shielding constants. Here, we explore the use of hybrid (cTPSSh) and range-separated hybrid (cTPSSrsh) forms of the cTPSS functional for the evaluation of this quantity. For the cTPSSh functional, we add 10% orbital-dependent exchange to 90% of the cTPSS exchange functional.
The construction of cTPSSrsh follows that of Ref. 65. In particular, we range separate the exchange component so that the exchange energy per particle becomes yielding the cTPSS contribution to the exchange energy, where cTPSS ensuring that, for large inter-electronic separations r 12 , the exchange integrals, which are evaluated with the [erf (µr 12 )] /r 12 operator, approach 100% orbital-dependent exchange. MAEs of 24.6 and 25.3 ppm respectively. In the remainder of this paper, we will examine current densities associated with chemical bonds.
At bond mid-points, these currents are sensitive to the delocalisation of electronic charge.
Delocalisation errors associated with GGA-type functionals 70 can lead to less accurate current densities, which to some extent can be corrected for by utilising hybrid functionals.
Combining this observation with the results obtained from the evaluation of NMR shielding constants, we use the cTPSSh functional in the remainder of this study.

Quadrature schemes for analysis of magnetically induced currents
The physical current is a complicated vector field. Whilst its features potentially reveal a large amount of chemical information, tools are required to aid in its analysis and interpretation. The GIMIC program 21,22 provides flexible quadratures to allow for integration of the current passing user-defined planes, giving measures of ring and bond currents in molecular system. 2, 21 We have added similar functionality to the QUEST program. 48 As an example, consider the naphthalene molecule, shown in Figure 1.
On the left, the current is plotted for a plane 1a 0 above the molecular plane. To determine this current, a weak field of 0.001 B 0 was applied in the direction normal to the molecular plane. There is a strong current around the perimeter of the molecule following the C-C bonds. A common way to measure this global ring current is to setup integration planes.
In the right-hand panel of Figure 1, the red line indicates the position of one such plane, bisecting 3 C-C bonds and extending 10a 0 above and below the molecular plane. Integration over the whole plane gives zero current by symmetry, as seen from the bond-current profile, Whilst the naphthalene molecule is a simple planar system, it highlights two common issues encountered when attempting to analyse molecular currents by setting up local quadratures. Firstly, the integrated current susceptibility calculated is dependent on the area of the integration plane-in this case, a large plane is used to capture the whole ring current value. Secondly, it may not always be possible to use such large integration planes without intersecting another bond vector-for example, the plane used for the C9-C10 bond from ring centre to ring centre will have a different spatial extent to that starting at a ring centre, bisecting the C12-C13 bond and continuing to a distance far from the molecule. This issue is commonly encountered for more complex structures, particularly if they are non-planar. We now consider how to set up bond-centred quadratures that have consistent areas in general molecular systems.

A disc-based quadrature
For the 2D square planar integration discussed so far, Gauss-Legendre quadrature is used, similar to that employed by the GIMIC program. 21,22 Here we set up bond-centred disc quadratures using the Elhay-Kautsky method, 71,72 where the integral of a function F in the xy plane is calculated as where r c is the radius of the disc, n θ is the number of angular nodes, n r is the number of radial nodes and w i are the quadrature weights. It is recommended to use n θ = 2n r for a balanced integration of angular and radial coordinates. In our implementation, we use the structural analysis routines in QUEST to identify all bonded atom pairs. For each bond, the average of the covalent radii of the atoms is calculated and used as r c . The quadrature is then constructed for each bond initially in the xy plane, before being translated to the bond centre and rotated so that the normal to the centre of the disc lays along the bond vector.
In this manner, small 2D disc quadratures can be rapidly constructed for each bond in the molecule.
The result of this procedure is illustrated for the benzene molecule in Figure 2. Here a small quadrature with n θ = 10 and n r = 5 is shown for illustration purposes. To test the utility of this quadrature to distill complex current-density vector fields into bond current susceptibilities, we have applied it to a range of previously studied planar ring structures. The results are shown in Figure 3 and Table 1, we emphasize that the   Figure 3: Bond current susceptibilities assigned by use of the disc-based quadrature (nA/T) for cyclobutadiene, benzene, pyridine, naphthalene, pentalene and bispentalene annelated naphthalene. C-H currents are shown in green, C-C currents in red and C-N currents in orange. All values calculated at the cTPSSh/6-31G* level.
A more challenging case is bispentalene annelated naphthalene (BPAN), previously studied using GIMIC by Sunholm, Berger and Fliegl. 83 Cao et al. 26 synthesised and characterised organic compounds with two pentalene units annelated with a naphthalene moiety in 2015.
As the molecule has 22π electrons, Hückel's rule would predict the molecule to be aromatic.
However, an upfield shift of pentalene hydrogen atoms was measured using nuclear magnetic resonance spectroscopy (NMR) indicating anti-aromaticity. The disc quadrature used in the present work leads to conclusions similar to those obtained via the GIMIC analysis in Ref. 83-namely, that the pentalene moieties remain strongly anti-aromatic, whilst the central naphthalene moiety is weakly aromatic. These observations have been used to rationalise the experimental observations of Cao et al. 26 Here they further establish the validity of our rents are quoted as a dimensionless value relative to the benzene value. As discussed earlier the planar quadrature approach yields ring currents for naphthalene and benzene whose ratio is very close to the HLPM value of 1.093. A complicating factor in making comparisons is that the ab initio results presented here vary according to subtle changes in bond length.
However, the cTPSSh/6-31G* bond current susceptibilities in Table 1 give the corresponding ratio as being in the range 1. We emphasize that the simple quadratures used are automatically set up, requiring no extra input from the user for positioning or determining the extent of the planes. All disc Table 1: Bond current susceptibilities for the planar ring systems cyclobutadiene, benzene, pyridine, naphthalene, pentalene and bispentalene annelated naphthalene (nA/T). The atomic numbering is shown in Figure 3. quadratures are determined directly from the structural analysis and bond identification that is a standard part of almost all quantum-chemical programs. Furthermore, since all bond types are treated with the same averaged covalent radii, the intensity of the currents may be consistently compared between systems. In Table 1, for example, we see that the C-C bond current susceptibilities in the anti-aromatic systems are generally more intense than those in the aromatic systems.

3D structures
A potential advantage of the disc-based quadrature is its utility for compact 3D structures.
To explore this, we consider the ring, bowl and cage isomers of C 20 as prototypical systems.
We use geometries optimized at the PBE0/6-31G* level using density fitting in the df-def2 auxiliary basis and the ADMMS approximation in the 3-21G ADMMS auxiliary basis, yielding structures close to D 10h , C 5v , and D 3d symmetries for the ring, bowl and cage respectively. The stability of these isomers has been studied extensively, 85-87 using a range of quantum-chemical methods, including accurate coupled-cluster methods. 87 In Figure 4, we show the disc-based quadratures for the bowl and cage isomers, as prototypical 3D cases. The bowl exhibits relatively weak curvature and the quadratures remain well separated. In the right panel of Figure 4, only the three nearest disc quadratures are shown for clarity. Reassuringly, none of the disc quadratures intersect in this relatively compact structure, confirming the applicability of this quadrature to general systems.
The bond current susceptibilities for each isomer are shown in Figure 5.  Table 2 for details.   Figure 5: Bond currents, assigned using the disc-based quadrature for the ring (left) and bowl (centre) isomers of C 20 . The atomic numbering is shown for the cage isomer and the bond current values are shown in Table 2 Table 2: Bond current susceptibilities for the bowl, cage and ring isomers of C 20 . The atomic numbering is shown in Figure 5.

Magnetically Induced Current Susceptibilities in Excited States
The disc-based and planar quadratures in QUEST provide a flexible suite of tools for the analysis of magnetically induced current densities in molecular systems. We have also implemented the analysis of spin-resolved current densities. This is particularly important if one wishes to examine the current densities of not only ground but also excited states. A prototypical example is the benzene molecule. In the ground state, the α and β spin currents are the same; these are plotted in Figure 6. In the left panel, the in-plane currents show the current pathways in the σ framework, whilst the currents 1a 0 above the plane show the π currents. Figure 6: Magnetically induced currents (a.u.) for the benzene molecule in its ground state, for a field of 0.001 B 0 perpendicular to the molecular plane. The current is plotted in the molecular plane (left) and 1a 0 above the molecular plane (right). Note that the plot range is restricted so as to exclude the intense induced currents in close vicinity of the nuclei.
Papadakis and Ottosson 88 have highlighted the 'Jekyl and Hyde' character of the benzene molecule, with its first triplet excited state exhibiting strong anti-aromaticity according to Baird's rule. 89 The spin-resolved magnetically induced currents are shown in Figure 7 for a field of 0.001 B 0 perpendicular to the molecular plane. In such a field, the first excited state with two unpaired β electrons is lowest. The spin-resolved currents of this state can be directly accessed via an SCF calculation. The α current in the left-right panel of Figure 7 is slightly more compact relative to the ring centre than the β current in the right-hand panel. The intensity of the α current is also lower, reflecting the larger population of β spin electrons.

Magnetically Induced Currents in Strong Fields
An advantage of the non-perturbative approach to calculating magnetically induced currents is the ability to study systems explicitly as a function of field strength, beyond the perturbative regime. The BH molecule is a classic example of a closed-shell system exhibiting paramagnetism. This paramagnetism has been rationalised in terms of a simple two-state model in Ref. 90, which leads to a ground-state energy that first decreases in the presence of a magnetic field perpendicular to the bond axis, before rising diamagnetically. This behaviour is shown in the insets in Figure 8.  Table 3, we present the potential energy curve data for the spherical-harmonic primitive basis sets aug-cc-pVXZ with 3 ≤ X ≤ 6. For each field strength, the basis-set error in kcal/mol relative to the primitive aug-cc-pV6Z basis is given in parentheses, estimated in a perpendicular magnetic field. Our estimates confirm that basis-set errors increase with increasing field strength, as expected for isotropic basis sets, but remain below 1 kcal/mol up to 0.15 B 0 for X = 3, up to 0.70 B 0 for X = 4 and over the entire range for X = 5, which has a maximum error estimate of 0.31 kcal/mol at 1.00 B 0 . Bearing in mind that these isotropic basis sets are not optimized at all for performance in a magnetic field, the results are reasonable. To maintain basis set errors below 0.1 kcal/mol at field strengths above 0.5 B 0 , it may be necessary to either consider optimization/augmentation of isotropic basis sets or the use of anisotropic Gaussian basis sets. We emphasize that the weak-field magnetic current susceptibilities analysed in the previous sections are insensitive to such basis-set effects.
In Figure 8,  continues at higher fields as the energy of this state continues to rise.

Conclusions
We have introduced a suite of tools to analyse the complex current vector field induced by exposing a molecule to an external magnetic field. In the perturbative regime, we have implemented the well-established tools of the GIMIC program, 21,22 using 2D planar Gauss-Legendre quadrature. In addition, we propose an alternative disc-based quadrature for analysing bond current susceptibilities. This quadrature may be easily and automatically set up based on the structural analysis carried out by any quantum chemical program. The radii of the discs are chosen as the average of the covalent radii of the atoms in each bond.
We demonstrated that, for a range of planar ring systems, the qualitative insights offered by this procedure mirror those of the 2D planar Gauss-Legendre quadrature. Furthermore, the bond current susceptibilities were shown to provide an accurate distillation of the complex features of the current density to simple chemical diagrams of the type in Figure 3. A key advantage of the proposed disc quadrature is its applicability to 3D structures. This was demonstrated for the ring, bowl and cage isomers of the C 20 molecule, in which the disc quadratures remain well defined and deliver bond current susceptibilities that match with expectations based on the symmetry of these systems.
The disc-based quadrature is intended as a complement, rather than replacement for rectangular quadratures. Whilst it has advantages for the analysis of currents assigned to bonds and the quadratures are consistent between structures, the rectangular quadratures can be appropriate for analysing the strength of ring currents, for example, which are a more global feature of molecular systems. The rectangular quadratures also offer more natural decomposition into bond current susceptibility profiles. The flexibility of our implementation for open-shell systems and excited states was demonstrated for the benzene molecule in its ground and first excited state, the former having aromatic character and the latter having strongly anti-aromatic character.
Finally, the utility of these tools for the anaylsis of magnetically induced currents was demonstrated for the BH molecule. At weak fields the currents reflect those obtained from the analysis in using response theory as in Ref 2. At higher fields, the currents reflect the transition from para-to dia-magnetism, an effect that cannot be visualised using the linear-response approach to study magnetically induced current susceptibilities. We expect that the alternative quadrature and tools for strong-field current analysis will provide useful methods for the interpretation of calculations in the presence of strong magnetic fields in the future.

Supplementary Information
Isotropic NMR shielding constants for the LDA, PBE, PBE0, cTPSS and cTPSSh functionals are tabulated in the supplementary information.