Current Density-Functional Theory using meta-Generalized Gradient Exchange--Correlation Functionals

We present the self-consistent implementation of current-dependent (hybrid) meta generalized gradient approximation (mGGA) density functionals using London atomic orbitals. A previously proposed generalized kinetic energy density is utilized to implement mGGAs in the framework of Kohn--Sham current density-functional theory (KS-CDFT). A unique feature of the non-perturbative implementation of these functionals is the ability to seamlessly explore a wide range of magnetic fields up to 1 a.u. ($\sim 235000$T) in strength. CDFT functionals based on the TPSS and B98 forms are investigated and their performance is assessed by comparison with accurate CCSD(T) data. In the weak field regime magnetic properties such as magnetizabilities and NMR shielding constants show modest but systematic improvements over GGA functionals. However, in strong field regime the mGGA based forms lead to a significantly improved description of the recently proposed perpendicular paramagnetic bonding mechanism, comparing well with CCSD(T) data. In contrast to functionals based on the vorticity these forms are found to be numerically stable and their accuracy at high field suggests the extension of mGGAs to CDFT via the generalized kinetic energy density should provide a useful starting point for further development of CDFT approximations.


I. INTRODUCTION
The foundations of current density-functional theory (CDFT) and its Kohn-Sham (KS) implementation were established in the late 1980's with the seminal works of Vignale, Rasolt and Geldart [1][2][3] , where it was recognised that the exchange-correlation functionals must depend not only on the electronic density, ρ, but also the paramagnetic current density j p in the presence of an electromagnetic field.Since these early works a large number of theoretical investigations of CDFT have been presented.The foundations of CDFT have sometimes been viewed as controversial.Most recently, Pan and Sahni [4][5][6] suggested that the physical current density j, rather than j p , aught to be the fundamental variable in a CDFT and attempted to establish a Hohenberg-Kohn like theorem for this physically appealing alternative choice of variable.Unfortunately, the derivation of this theorem has been shown to be in error 7,8 .Furthermore, the work of Tellgren et al. 7 showed how CDFT may be brought into Lieb's convex-conjugate formulation of DFT 9 , further strengthening its foundations and lending key insight into the more complex relationship between the key densities and potentials in the theory.In particular, it is highlighted that the lack of a Hohenberg-Kohn theorem is not an impediment to a viable CDFT.Recent theoretical works a) Electronic mail: pcxjf1@nottingham.ac.uk b) Electronic mail: andrew.teale@nottingham.ac.uk by Lieb and Schrader 10 and Tellgren et al. 11 have also addressed the issue of N -representability in CDFT.
Despite the theoretical progress in CDFT very few practical implementations of theory have been presented.Most practical studies have either presented calculations based on fixed densities (typically computed at the Hartree-Fock or standard Kohn-Sham level), or have attempted to include CDFT contributions in linear-response calculations.In the context of response theory, implementations have been presented for magnetic properties by Lee et al. 12 and for excitation energies at the meta generalized gradient approximation (mGGA) level by Bates and Furche 13 .Very few fully self-consistent implementations of CDFT capable of treating systems beyond the linear-response regime have been presented, in fact we are only aware of the work by Pittalis et al. 14,15 in the context of the optimized effective potential method, Zhu and Trickey 16 for atomic systems and our own implementation 17 for general atomic and molecular species.
A number of challenges arise when implementing quantum chemical methods for molecules in magnetic fields, the London program 18 has been specifically designed to address these and is utilized throughout the present work.In particular, London atomic orbitals [19][20][21] are employed to ensure gauge origin invariant results.For CDFT additional challenges arise since new forms are required for the exchange-correlation functional.Relatively few practical forms for CDFT functionals have been suggested in the literature.
In the present work we will examine the use of mGGAs and hybrid mGGAs for the exchange-correlation energy in the presence of magnetic fields.Functionals in this class depend on the orbital dependent kinetic energy density in the absence of a magnetic field.However, as has been noted in the literature 13,22,23 , this key quantity is not gauge invariant and so some modification is required for use in a magnetic field.One approach is to replace the kinetic energy density by a generalized form including the paramagnetic current density.This quantity naturally arises in the expansion of the spherically averaged exchange hole, as derived by Dobson 24 .Becke has already suggested the use of this approach to produce a current dependent generalisation of the Becke-Roussel 25 functional.Recently, Bates and Furche 13 have also explored a similar generalization of the Tao-Perdew-Staroverov-Scuseria (TPSS) functional 26 to calculate excitation energies via response theory.
We will consider current dependent extensions of the B98 27 , TPSS 26 and TPSS(h) 28 functionals.The use of a modified current-dependent kineitc energy density is denoted by a prefix c throughout the remainder of this work.The non-perturbative nature of the implementation in the London program will allow for testing of these functionals in both weak and strong field regimes.The availability of accurate ab initio methodologies in the London program provides a unique opportunity for the assessment and testing of CDFT functionals at field strengths upto 1 a.u.In the present work we make use of the recent implementation by Stopkowicz et al. 29 of coupled-cluster (CC) methods with single, double and perturbative triple excitations [CCSD(T)] for benchmarking the CDFT approximations.
In Section II we review the simple generalization of mGGA functionals to the CDFT framework due to Dobson and Becke, details specific to the functionals considered in this work are collected in the appendix.In Section III we outline some computational details of our calculations.Section IV summarizes our findings, assessing the quality of the CDFT approximations by comparison with CCSD(T) data; in Section IV A we explore the performance of mGGA functionals for calculating molecular properties in the weak field regime accessible via linear response theory; in Section IV B the high field regime is explored by considering the recently proposed perpendicular paramagnetic bonding.The interpretation of this bonding mechanism in the KS-CDFT framework is discussed in Section V. Finally, concluding remarks and directions for future work are given in Section VI.

II. THEORY
In this work we consider the calculation of energies and molecular properties in the presence of a static uniform external magnetic field, B, which may be represented in terms of the vector potential where R G is an arbitrary gauge origin.The London program makes use of London atomic orbitals [19][20][21] to ensure that computed energies and molecular properties are invariant with respect to choice of the gauge origin.These basis functions take the form (2) where χ µ (r K ) is a standard Gaussian-type orbital centred at position R K .These perturbation-dependent basis functions are used to expand the Kohn-Sham molecular orbitals.
The Kohn-Sham approach to density-functional theory can be extended to CDFT by searching for a noninteracting system of electrons with the same charge and current densities as the physical interacting system: (ρ s (r), j p,s (r)) = (ρ(r), j p (r)).The KS equations can be written as where the KS potentials (u s , A s ) are defined as and The first three terms in the scalar potential u s represent the external (ext), Coulomb (J) and exchange-correlation (xc) potentials defined as the functional derivative of the respective energies with respect to the density, as in the usual KS approach.The final contribution to the scalar potential arises from the non-interacting vector potential A s .In addition to the vector potential due to the external field A ext as defined in Eq. (1) the KS vector potential contains an exchange-correlation contribution defined as A central question that immediately arises in CDFT is how the exchange-correlation functional must be modified to include current effects.Whilst the paramagnetic current density is a valid quantity on which to base the universal density functional it can also be shown that the exchange-correlation component must be independently gauge-invariant 2 .This places a significant constraint on the manner in which this quantity may enter any approximate CDFT functional.In contrast to standard DFT, relatively few CDFT functionals have been proposed.The majority of these are based on the the vorticity with as proposed by Vignale, Rasolt and Geldart (VRG) 3 .The original VRG form for the exchange-correlation energy was parameterized using Monte Carlo simulations of the high density limit 1 .A number of re-parameterisations for this form have been suggested based on accurate calculations in the high-density regime 12,[30][31][32] .Higuchi and Higuchi 33 (HH) have also presented a vorticity dependent form, derived to obey known exact relations for the CDFT exchange-correlation functional.
Whilst the vorticity is a theoretically convenient choice to ensure the gauge invariance of the exchange-correlation energy it has been observed that in practical self-consistent calculations it can lead to significant numerical stability issues 16,17 .How severe these issues are depends on the exact parameterization of the the functional form, however, in all cases some degree of numerical regularization is required to ensure that the self-consistent field solution of the Kohn-Sham equations can be obtained.Furthermore, molecular properties computed by such calculations display an un-acceptably strong dependence on the regularization parameters -with no obvious convergence towards a single value.Clearly this raises questions as to how appropriate such forms are for use in quantum chemical calculations.
Most practical mGGAs make use of the kinetic energy density in their construction, where ϕ i are the occupied KS orbitals and σ is the electron spin index.This term is gauge dependent and as a result an unmodified meta-GGA type functional form cannot be used to describe a system with a non-zero magnetic vector potential.To resolve this issue the gauge independence of the exchange-correlation functional must be restored.A natural modification, which can be applied to any mGGA dependent on the kinetic energy density τ (r), arises in the work of Dobson 24,34 who generalized the expansion of the exchange-hole to include the case of non-zero current densities.The spherically averaged exchange hole at zero field can be modelled using a Taylor expansion 35 and is commonly considered 25,27,36 up to the quadratic term This expansion can be generalised to non-zero field and the curvature term becomes 24,25,34 where j pσ is the paramagnetic current density Comparing Eqs.(10) and (11) it is possible to identify a correction to the conventional τ (r) that is gauge invariant and may be utilized in mGGA functionals, The use of Eq. ( 13) has been put forward many times in the literature.Becke suggested its use in the Becke-Roussel model 25 to generate a current dependent analogue of this functional.He also suggested that this quantity could be used to define a current dependent inhomogeneity parameter in the more empirical B98 functional 27 .It has also been suggested for use to generalize the TPSS functional 26 by Tao 23 .Recently Bates and Furche 13 considered the application of the resulting cTPSS functional in the calculation of excitation energies via response theory.In the present work we consider the application of mGGA functionals with this modification to calculation magnetic properties in the weak and strong field regimes in a non-perturbative manner.In particular we consider three functional forms, cB98, cTPSS and cTPSS(h), where the prefix c denotes the use of the modified τσ in Eq. ( 13).The Appendix gives some details of these functional forms to show how these modifications enter.

III. COMPUTATIONAL DETAILS
Unless otherwise indicated all calculations in this work use the aug-cc-pCVTZ basis set 37,38 .All DFT calculations have been performed using the London 18 program.This code utilizes the XCFun library 39 for the evaluation of the density functionals and their derivatives.The modifications of Eq. ( 13) and the functionals cB98 and cTPSS have been added to the XCFun library.In addition we investigate the use of a hybrid form of cTPSS, denoted cTPSS(h), based on the TPSS(h) functional of Ref. 40.The quality of the CDFT functionals cB98, cTPSS and cTPSS(h) is assessed by comparison with CCSD(T) data.For comparison Hartree-Fock (HF), local density approximation (LDA), Perdew-Burke-Ernzerhof (PBE) 41 , and Keal-Tozer-3 (KT3) 42 density-functional results are also presented.The latter is of particular interest since it is specifically designed for the calculation of nuclear magnetic resonance (NMR) shielding constants.
The performance of these approximations will be considered in two regimes; the weak field regime accessible by linear response calculations and the strong field regime only accessible via non-perturbative calculations.In the weak field regime we will consider the calculation of molecular magentizabilities and NMR shielding constants for the 26 small molecules in Table I.Errors for these quantities are presented relative to the benchmark data of Ref. 43. Results are also compared with those including corrections from the Tao-Perdew parameterization 31 of the Vignale-Rasolt-Geldart functional 2 , taken from Ref. 17.In the strong field regime we consider the prediction of perpendicular paramagnetic bonding 44 in a IV.RESULTS

A. The weak field regime: magnetic properties
We commence by considering the molecular magnetizabilities and NMR shielding constants of the 26 small molecules in Table I.The magnetizability tensor elements, ξ α,β , are defined as where α and β label cartesian components of the tensor and magnetic field.The NMR shielding tensor for a given nucleus K is defined by where M K is the nuclear magnetic moment of nucleus K.These properties can be accessed non-perturbatively in the London program by explicit calculation of the energy in the presence of the perturbing fields.Details of this procedure are given in Ref. 17, here we compute values for each method considered in the same manner, facilitating a comparison with previous results.
Given that for many density-functional approximations these singlet second order magnetic response properties can be accessed by standard linear response methods in a variety of programs, this approach may seem cumbersome.However, it should be noted that the implementation of the new CDFT approaches in this framework is much more straightforward, requiring only an implementation of the functional and the derivatives required for construction of the KS matrix.More importantly, as we will see in Section IV B, this non-perturbative approach allows us to seamlessly explore the behaviour of new approximations in much stronger fields -inaccessible via linear response theory.This means that London provides a powerful testbed for new CDFT functionals.
To quantify the accuracy of the DFT approaches for the calculation of these properties we compare our results with the CCSD(T) benchmark values of Ref. 43.Specifically, we use the values at the CCSD(T)/aug-cc-pCV[TQ]Z level -which have been extrapolated to the basis set limit using the procedure of Refs.43 and 45.
The errors in the calculated magnetizabilities and NMR shielding constants are summarised in Table II and presented graphically in Figures 1 and 2 as box-whisker plots.In these plots individual points represent the errors for each system relative to the reference values, the upper and lower fences of the whiskers denote the maximum positive and negative errors respectively and the coloured boxes enclose errors between the 25% and 75% quantiles.Mean and median errors are marked in the plots by horizontal black and white lines, respectively.Grey diamonds are used to represent the confidence intervals.
For the molecular magnetizabilities it is clear from the error measures in Table II and their representation in Figure 1 that none of the functionals offers high accuracy.The GGA functionals PBE and KT3 in particular do not offer significant improvements over LDA.Whilst their minimum and maximum errors are slightly improved, the mean errors actually deteriorate.Similar observations were made in Ref. 45 for these type of functionals.The B97-2 functional gives slightly reduced errors and this is consistent with previous conclusions 45 that for magnetizabilities the inclusion of HF exchange may be beneficial.At the GGA level the underlying functionals are already gauge invariant but do not depend explicitly on the paramagnetic current density.To introduce this dependence the VRG functional may be added.In our earlier work 17 we found that this correction can be numerically problematic and that the most stable parameterization of this functional to date is that put forward by Tao and Perdew 31 , denoted VRG(TP).For comparison we include here the PBE+VRG(TP) and KT3+VRG(TP) results.It is clear that the inclusion of the VRG(TP) correction actually worsens the agreement of the results with CCSD(T).At the mGGA level the inclusion of current is mandatory to ensure the exchange-correlation evaluation is gauge independent.
Here we investigate the cB98, cTPSS and cTPSS(h) functionals, the TPSS and TPSS(h) forms are similar in performance to PBE -offering marginal improvements on some error measures, with TPSS(h) performing slightly better than TPSS.The cB98 form gives the best performance of all the functionals considered.It is noteworthy that in the mGGA class the mean errors reduce as more HF exchange is included in the functional -with cTPSS, cTPSS(h) and cB98 containing 0%, 10% and 19.85% HF exchange respectively.Suggesting that for this property that the treatment of exchange may be the dominant factor in the errors.Since the treatment of long-range exchange is not rectified in the transition from GGA to mGGA type functionals (and only partially corrected by a global admixture), then it may be that this factor far out weighs any improvements due to the inclusion of current effects.
It is worth emphasizing that in the course of our investigation we found that the implementation of the mGGA functionals including a generalized kinetic energy density was straightforward.In particular we found that no special care was required with respect to numerics compared with standard functionals and that in practical use the functionals are robust and self-consistent calculations using these functionals converge without significant difficulty.This sharply contrasts the behaviour for the VRG functionals as investigated in Refs.16 and 17.
The results for NMR shielding constants are presented in Figure 2.Here we see LDA is poor as expected.In addition KT3, which was designed for these properties, is the best performing functional -significantly improving over the standard PBE GGA functional and the B97-2 hybrid functional.In this case we see the addition of the VRG(TP) correction to PBE and KT3 has little effect, very slightly deteriorating the results.The mGGA results for cB98, cTPSS and cTPSS(h) are intermediate between PBE and KT3 -offering small systematic improvements over PBE.Again B98 produces the best results of the mGGA functionals.
On the whole the quality of the mGGA results at modest field strengths may be regarded as disappointing.The overall errors suggest that mGGAs may offer modest improvements over conventional GGA functionals such as PBE -though they cannot compete with GGAs tailored to specific properties.This is broadly consistent with findings by Bates and Furche for the calculation of excitation energies using cTPSS 13 .Since current corrections are known to be relatively small it is important that the underlying functional should be relatively accurate.For the mGGAs considered here there are known weaknesses (for example in the treatment of long-range exchange) that may obscure the effect of the current dependence.For the case of NMR shielding constants a more detailed analysis of the significance of current effects and how these interplay with errors in a range of density functionals is presented in Ref. 46.We will now examine how these functionals perform when the magnetic field becomes much higher and has a stronger effect on the electronic structure.

B. The strong field regime: paramagnetic bonding
One approach to explore whether or not the inclusion of current effects via the modified kinetic energy density of Eq. ( 13) is physically reasonable is to increase the strength of the magnetic field.Lange et.al. 44 have recently performed full configuration-interaction (FCI) calculations at high field that have uncovered a new mechanism for chemical bonding in the presence of a strong magnetic field.This new bonding has been termed perpendicular paramagnetic bonding and occurs at field strengths similar to those found on some white dwarf stars.Since this work Murdin et al. 47 have shown that phosphorus and selenium doped silicon semiconductors can produce a viable laboratory analogue of free hydrogen 47 and helium 48 in strong magnetic fields.The description of these types of systems via quantum chemistry will require less computationally demanding approaches -and CDFT is one strong candidate for the simulation of these systems.
To investigate the performance of cB98, cTPSS and cTPSS(h) in strong magnetic fields potential energy profiles were calculated for H 2 , He 2 , NeHe and Ne 2 .In particular, we consider the 3 Σ + u (1σ g 1σ * u ) state of H 2 and the lowest 1 Σ + g states of He 2 , HeNe and Ne 2 .Each of these states is repulsive or weakly dispersion bound in the absence of a magnetic field but become more strongly bound when a field is applied.We note that only the 3 Σ + u (1σ g 1σ * u ) of H 2 is an overall ground state in the presence of the field.These states were compared against results from accurate CCSD(T) potential energy curves calculated using a recent non-perturbative implementation by Stopkowicz et al 29 in the London program.For comparison we have also generated similar profiles with standard LDA and GGA density functionals as well as with HF theory.The calculated potential energy curves for H 2 , He 2 , NeHe and Ne 2 are shown in Figures 3, 4, 5 and 6, respectively.Equilibrium bond lengths and dissociation energies were determined numerically and are presented in Table III.
The 3 Σ + u (1σ g 1σ * u ) state of the H 2 molecule in a perpendicular field was examined at the FCI level by Lange et al. 44 .The potential energy curves for this state are shown in Figure 3. HF strongly underbinds this state in comparison with the FCI data.In contrast LDA strongly overbinds, a tendency which is largely corrected by the PBE functional and further improved by the cTPSS and cTPSS(h) models.These trends are reflected in the equilibrium bond lengths and dissociation energies in Table III.
Although not highly accurate the cTPSS and cTPSS(h) models give a reasonable qualitative description of the potential energy curve.The empirically parameterized KT3 functional is interesting because it gives simultaneously a reasonable estimate of both the equilibrium bond length and dissociation energy.However, at intermediate separation an unphysical barrier is observed.For B98 an even more pronounced barrier is present and the potential energy curve is generally even less accurate than Hartree-Fock theory.This may suggested that heavily parameterized functional forms, determined to perform well at zero field, may not be the best candidates for use in strong-field CDFT studies.Examining the potential energy curves for He 2 in Figure 4 we see that HF tends to under-bind with a bond length of 3.30 a 0 compared with the CCSD(T) value of 2.98 a 0 .Similarly the HF dissociation energy of 0.218 mE h is much smaller than the corresponding CCSD(T) value of 1.259 mE h .For this small system we were able to compare the CCSD(T) results with FCI values (calculated in the slightly smaller aug-cc-pVTZ basis), as expected the agreement is excellent -the corresponding potential energy curves are essentially indistinguishable on the scale of Figure 4.For LDA we see a strong tendency to overbind giving much too short R e values and much too large estimates of D e .The GGA functionals PBE and KT3 show considerable improvement over LDA, however, they still strongly overbind.The improvement for the mGGA functionals is striking -in particular TPSS and TPSS(h) give a good qualitative description of the potential energy curve.The corresponding R e and D e values indicate that there still remains a tendency towards over binding but this is greatly reduced.
The cB98 functional tends to show more significant under binding.Here we note that the arguments used in the construction of cB98 and cTPSS are rather different.In particular, cB98 is an empirically parameterized functional (see the appendix), whereas cTPSS is constructed based on the satisfaction of known exact conditions.In this work we have used the parameters determined in Ref. 27 to define the cB98 form.These parameters were determined at zero field and from post-LDA calculationsas a result they may not be optimal for fully self-consistent calculations in the presence of a magnetic field.On the other hand the cTPSS functional is designed to satisfy selected constraints at zero field and it could be argued that in the presence of a field both the B98 and TPSS based functionals are open to further optimization, though this is beyond the scope of the present work.
The stability of the mGGA functionals is particularly evident in the the strong field regime when one compares the present results for He 2 with those for the VRG-based estimates in Figure 7 of Ref. 17.The VRG approaches led to very difficult SCF convergence and complex potential energy curves with a strong unphysical over binding.The mGGAs considered here are un-problematic in practical application and yield results surprisingly close to the CCSD(T) estimates.
Similar qualitative trends are observed for the NeHe and NeNe dimers in Figures 5 and 6.Again LDA and GGA functionals are not sufficiently accurate for practical use and the mGGA functionals provide a large improvement.The cTPSS and cTPSS(h) results remain impressivewith cTPSS(h) being consistently slightly more accurate than cTPSS.This trend is reflected in both the potential energy curves and the R e and D e values in Table III.

V. INTERPRETATION OF PARAMAGNETIC BONDING IN THE KS-CDFT FRAMEWORK
The mGGA CDFT functionals offer a computationally cheap correlated method for the examination of the exotic bonding mechanisms observed in a strong magnetic field.In many areas of chemistry the nature of bonding, chemical reactions, spectra, and properties of molecular species are interpreted qualitatively in terms of orbital g state of Ne2 in a magnetic field of 1 a.u.perpendicular to the bonding axis for a variety of methods with the aug-cc-pCVTZ basis set.
interactions.We now consider the extent to which information from KS-CDFT calculations can aid in simple interpretation of the perpendicular paramagnetic bonding interactions.
We begin by considering a molecular orbital analysis of the perpendicular paramagnetic bonding.KS-CDFT calculations provide a simple set of canonical molecular orbitals, which can be used to construct the electronic density via Here we note that the occupied KS orbitals ϕ i (r) can be complex in the presence of a magnetic field.In the second equality the density is expressed in terms of the one-particle density matrix D γζ and the basis functions ω γ (r).We will commence by considering how the molecular orbital energies associated with H 2 and the rare gas dimers change upon application of a magnetic field as the perpendicular paramagnetic bonding in Section V A evolves.Since the orbitals themselves can be complex in the presence of a field we then proceed in Section V B to analyze the bonding in terms of the changes in electronic density of Eq. ( 16) as a function of field, which is naturally a real observable quantity.
A. KS-CDFT molecular orbital analysis KS molecular orbitals have been widely used as an interpretive aid in chemical applications throughout the literature.The KS orbitals are defined to minimize the non-interacting kinetic energy and yield the physical electronic density via Eq.( 16).They also have appealing properties; for example the highest occupied MO energy is minus the first ionization potential (IP) 49,50 and the remaining orbital energies can be interpreted as Koopman's type approximations to higher IPs 51 .The extent to which these properties hold for general practical approximations has been a subject of debate in the literature 52,53 , as has the interpretation of KS virtual orbitals 54 owing to the role of the integer discontinuity 49 , which is missing from common approximations.However, from a practical standpoint it is widely accepted that interpretations based on occupied KS orbitals (to which we limit the following discussion) are theoretically justified and their utility has been borne out in many practical applications.
We now consider how the KS orbital energies change upon application of a perpendicular magnetic field of 1 a.u.Given that the cTPSS based models seem to be the most reliable of those studied in the present work we consider how the orbital energies from this functional change in Figures 7 and 8 for H 2 and He 2 , respectively.For the H 2 molecule in the 3 Σ + u (1σ g 1σ * u ) state we consider the energy of the occupied σ g and σ * u orbitals change upon bonding.In particular, we plot orbital energies in the absence of a field and in a perpendicular field of 1 a.u.relative to those of the atomic orbitals in the same field.We have plotted the orbital energies at an internuclear separation R e = 2.564 a.u.consistent with the cTPSS equilibrium bond length in the presence of a field.
We see that in the absence of a field the singly occupied σ g orbital is stabilized by 1.33 mE h , whilst the singly occupied σ * u is destabilized by 1.30 mE h relative to the 1s hydrogen orbitals.This is consistent with a net bond order of zero and a repulsive profile for the corresponding potential energy curve.In the perpendicular field of 1 a.u.we see that, relative to the hydrogen 1s orbital in the same field, the σ g orbital is stabilized by 63.21 mE h , whilst the σ * u orbital is destabilized by 34.97 mE h .This greater stabilization of the σ g orbital leads to a net bonding interaction, consistent with the analysis of Ref. 44.This illustrates that the KS orbitals can be a useful tool in rationalizing this exotic bonding phenomenon.
A similar analysis can be carried out for the lowest 1 Σ + g state of He 2 and is presented in Figure 8.In this plot we have separated the spin down and spin up orbitals and defined their energies relative to atomic orbitals of the same spin in the same field.Defined in this way the offset between the orbitals of different spin due to the Zeeman interaction is removed from the plot.Again the energies correspond to the cTPSS He 2 equilibrium internuclear separation R e = 2.864 a.u. in the presence of a perpendicular field.In the absence of a field we see that again the relative stabilization of σ g and destabilization of σ * u are approximately compensatory, whilst in the presence of a field the σ g orbital is more stabilized than the σ * u orbital is destabilized.It is also clear that the extra stabilization of the σ g orbital of ∼ 8 mE h is considerably less than the ∼ 28 mE h observed for H 2 .This is consistent with the strength of binding exhibited for these species in Figures 3 and 4 as well as in Table III.
Similar orbital energy diagrams may be constructed for the HeNe and NeNe systems, however, they become significantly more complex due to large differences between the orbital energies in HeNe and the splitting of the p-orbitals in NeNe.We therefore consider an alternative visualization of the bonding effects in these systems based on the charge and (physical) current density differences.

B. Electron density analysis
For each of the species He 2 , HeNe and Ne 2 , we have performed calculations of the electronic density and current density in the presence of varying perpendicular magnetic fields (B ⊥ ) using the cTPSS functional at the correspond- 8. Relative orbital energies for the occupied σg and σ * u CDFT molecular orbitals for He2 relative to the atomic 1s orbitals, with (blue) and without (red) a field of 1.0 a.u.perpendicular to the interatomic axis.
ing equilibrium geometries.Here we consider the density change ∆ρ B ⊥ (r) for each system relative to the isolated atoms in the same field, This difference density allows for a visualization of the nature of the paramagnetic bonding in these systems.In a similar manner one can consider the (gauge invariant) physical current density difference ∆j In Figure 9 we present plots of the density differences in Eqs. ( 17) and (18) for each of the species with B ⊥ = 1.0 a.u.The shading of the contours represents the buildup (red) or depletion of the charge density (blue), relative to two non-interacting atoms in the same field.In all three cases there is a clear build up of density between the atoms consistent with bonding.The charge density difference is elongated along the field, above and below the plane of the plots.The streamlines show the vector field associated with the current density difference.Paratropic circulations are clearly visible over the centre of the bonds where charge density accumulates, and diatropic circulations are visible in regions where the charge density is depleted.We have confirmed that for higher fields the alignment between para-/ dia-tropic circula-tions and charge accumulation / depletion becomes more pronounced.
This picture of the perpendicular paramagnetic bonding suggests that as the charge density is elongated along the field the constituent atoms may approach one another more closely.As they do so they experience a greater nuclear-nuclear repulsion, which may be screened by a rearrangement of the charge density towards the bond centre.This rearrangement is accompanied by consistent para-and dia-tropic current circulations.The cTPSS functional used in this work provides a simple, computationally cheap, route to perform analysis of the bonding encountered in the strong field regime.

VI. CONCLUSION
In this work we have implemented a previously detailed 23,24,34 modification to the kinetic energy density term in mGGA functionals to perform non-perturbative cDFT calculations in a stable, gauge invariant manner using the London program 18 .The modified mGGA functionals cB98, cTPSS and cTPSS(h) show a level of accuracy in predicting weak field magnetic properties that is competitive with existing GGA functionals without any additional fitting.The functionals cTPSS and cTPSS(h) show excellent prediction of perpendicular paramagnetic bonding behaviour, suggesting that the modification of the kinetic energy density is a viable route for the incorporation of current effects in standard mGGA densityfunctionals.In contrast to vorticity dependent forms previously studied 16,17 the functionals exhibited excellent numerical stability in the finite field setting without the need for delicate numerical regularizations.
Whilst the mGGA results show considerable promise in the high field regime, their performance in the weak field regime is perhaps disappointing -leading to only modest improvements of conventional GGA forms.To some extent this may be due to the fact that mGGAs contain many of the shortcomings associated with GGA forms.For example, the functionals still have potentials with incorrect asymptotic behaviour -particularly for the exchange contribution.Since the current effects at weak field are not dominant but rather add small corrections to the predominantly Coulombic exchange and correlation interactions then it may be necessary to further improve the underlying functional forms before the true impact of the current terms can be assessed in this regime.
We expect that this approach to include current dependence should play a central role in the future development of new CDFT functionals and many avenues are open for development.Obvious possibilities include the generalization of range-separated mGGA functionals 55,56 to obtain a better balance of errors between exchange, correlation and current contributions and the re-parameterization of functionals either empirically or via the consideration of alternative exact conditions in their construction.In the latter category we note that the presence of a magnetic This factor varies between 0 and 1 and vanishes in oneorbital regions, ensuring the functional is self-correlation free.For convenience in deriving fitted forms for the inhomogeneity correction factors g, Becke proposed the following transformation to a finite interval, where γ is a parameter to be determined and separate transformations are carried out for the exchange, oppositespin correlation and like-spin correlation respectively using the appropriate definitions of q as in Eqs (A3)-(A4).
Based on calculations of atomic exchange-correlation energies Becke proposed the values γ x,σ = 0.11, γ c,αβ = 0.14 and γ c,σσ = 0.16.The inhomogeneity corrections g were then determined by fitting a power series expansion of the form with m = 2 chosen to prevent unphysical over-fitting.This fitting was carried out using the G2 thermochemical dataset and basis set free, post-LSDA, calculations.The optimal parameters can be found in Ref. 27.
Here we use this parameterization directly but note that in future these parameters could be re-optimized based on self-consistent data.In addition an amount of Hartree-Fock exchange is included with weight c x = 0.1985.The resulting B98 functional therefore possesses a high degree of non-locality and may be classified as a hybrid mGGA functional with dependence not only on τ but also on the laplacian of the density ∇ 2 ρ.
The original definition of B98 utilized the zero-field exchange hole curvature of Eq. (10) in its definition.However, it was noted 27 that this form can be readily extended to include current effects via Eq.(11) and it is this avenue that we explore in the present work.Unless otherwise stated we employ this modified form throughout and denote it as cB98.

cTPSS
One of the most widely meta-GGA functionals is due to Tao, Perdew, Staroverov and Scuseria (TPSS) 26 .This functional is designed to satisfy exact constraints without empirical parameters and as such is an interesting candidate to study in the context of generalization to finite magnetic field strengths where much less is known about the performance of approximate functionals.In this functional the ratio Using the replacement in Eq. ( 13) leads to modifications in the exchange contribution via z in Eq. (A9) and in the correlation energy as shown in Eq. (A12).This modified form is noted cTPSS and is consistent with that used in the response implementation of Ref. 13.

2 FIG. 1 .
FIG. 1. Box-whisker plot of errors in (C)DFT molecular magnetizabilities relative to the extrapolated CCSD(T) values of Ref. 45.All calculations use the aug-cc-pCVTZ basis set.See the text for further details.

FIG. 2 .
FIG. 2. Box-whisker plot of errors in (C)DFT NMR shielding constants relative to the extrapolated CCSD(T) values of Ref. 43.All calculations use the aug-cc-pCVTZ basis set.See the text for further details.

FIG. 3 .
FIG.3.Potential energy curve for the 3 Σ + u (1σg1σ * u ) state of the H2 molecule in a magnetic field of 1 a.u.perpendicular to the bonding axis for a variety of methods with the aug-cc-pCVTZ basis set.

FIG. 4 .
FIG.4.Potential energy curve for the lowest 1 Σ + g state of He2 in a magnetic field of 1 a.u.perpendicular to the bonding axis for a variety of methods with the aug-cc-pCVTZ basis set.

FIG. 5 .FIG. 6 .
FIG.5.Potential energy curve for the lowest 1 Σ + g state of the NeHe dimer in a magnetic field of 1 a.u.perpendicular to the bonding axis for a variety of methods with the aug-cc-pCVTZ basis set.

FIG. 7 .
FIG.7.Relative orbital energies for the occupied σg and σ * u CDFT molecular orbitals for H2 relative to the atomic 1s orbitals, with (blue) and without (red) a field of 1.0 a.u.perpendicular to the interatomic axis.

TABLE I .
The test set of molecules for which accurate benchmark CCSD(T) data from Ref. 43 was available., He 2 , HeNe and Ne 2 .These non-perturbative calculations are assessed against CCSD(T) results computed using the implementation of Ref. 29 in the London program.

TABLE III .
Dissociation energies and equilibrium bond lengths for H2 and rare gas dimers in a 1 a.u.magnetic field perpendicular to the internuclear axis.The He 2 FCI calculations use the aug-cc-pVTZ basis set. a