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

###### Abstract

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. (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 Vignale and Rasolt (1987, 1988); Vignale, Rasolt, and Geldart (1988), where it was recognised that the exchange–correlation functionals must depend not only on the electronic density, , but also the paramagnetic current density 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 Pan and Sahni (2010, 2012, 2013) suggested that the physical current density , rather than , 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 Tellgren et al. (2012); Vignale, Ullrich, and Capelle (2013). Furthermore, the work of Tellgren et al. Tellgren et al. (2012) showed how CDFT may be brought into Lieb’s convex-conjugate formulation of DFT Lieb (1983), 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 by Lieb and Schrader Lieb and Schrader (2013) and Tellgren et al. Tellgren, Kvaal, and Helgaker (2014) have also addressed the issue of -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. Lee, Handy, and Colwell (1995) and for excitation energies at the meta generalized gradient approximation (mGGA) level by Bates and Furche Bates and Furche (2012). 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. Pittalis et al. (2006, 2007) in the context of the optimized effective potential method, Zhu and Trickey Zhu, Zhang, and Trickey (2014) for atomic systems and our own implementation Tellgren et al. (2014) for general atomic and molecular species.

A number of challenges arise when implementing quantum chemical methods for molecules in magnetic fields, the London program Tellgren, Soncini, and Helgaker (2008) has been specifically designed to address these and is utilized throughout the present work. In particular, London atomic orbitals London (1937); Ditchfield (1974, 1976) 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 Maximoff and Scuseria (2004); Tao (2005); Bates and Furche (2012), 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 Dobson (1993). Becke has already suggested the use of this approach to produce a current dependent generalisation of the Becke–Roussel Becke (1996) functional. Recently, Bates and Furche Bates and Furche (2012) have also explored a similar generalization of the Tao-Perdew-Staroverov-Scuseria (TPSS) functional Tao et al. (2003) to calculate excitation energies via response theory.

We will consider current dependent extensions of the B98 Becke (1998), TPSS Tao et al. (2003) and TPSS(h) Staroverov et al. (2003a) 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 a.u. In the present work we make use of the recent implementation by Stopkowicz et al. Stopkowicz et al. (2015) 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.1 we explore the performance of mGGA functionals for calculating molecular properties in the weak field regime accessible via linear response theory; in Section IV.2 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, , which may be represented in terms of the vector potential

(1) |

where is an arbitrary gauge origin. The London program makes use of London atomic orbitals London (1937); Ditchfield (1974, 1976) 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 is a standard Gaussian-type orbital centred at position . 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 non-interacting system of electrons with the same charge and current densities as the physical interacting system: . The KS equations can be written as

(3) |

where the KS potentials are defined as

(4) |

and

(5) |

The first three terms in the scalar potential 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 . In addition to the vector potential due to the external field as defined in Eq. (1) the KS vector potential contains an exchange–correlation contribution defined as

(6) |

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 Vignale and Rasolt (1988). 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

(7) |

with

(8) |

as proposed by Vignale, Rasolt and Geldart (VRG) Vignale, Rasolt, and Geldart (1988). The original VRG form for the exchange–correlation energy was parameterized using Monte Carlo simulations of the high density limitVignale and Rasolt (1987). A number of re-parameterisations for this form have been suggested based on accurate calculations in the high-density regime Lee, Handy, and Colwell (1995); Orestes, Marcasso, and Capelle (2003); Tao and Perdew (2005); Tao and Vignale (2006). Higuchi and Higuchi Higuchi et al. (2007) (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 Zhu, Zhang, and Trickey (2014); Tellgren et al. (2014). 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

(9) |

in their construction, where 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 , arises in the work of DobsonDobson (1992, 1993) 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 Becke (1983) and is commonly considered Becke and Roussel (1989); Becke (1996, 1998) up to the quadratic term

(10) |

This expansion can be generalised to non-zero field and the curvature term becomes Dobson (1992, 1993); Becke (1996)

(11) |

where is the paramagnetic current density

(12) |

Comparing Eqs. (10) and (11) it is possible to identify a correction to the conventional that is gauge invariant and may be utilized in mGGA functionals,

(13) |

The use of Eq. (13) has been put forward many times in the literature. Becke suggested its use in the Becke–Roussel model Becke (1996) 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 Becke (1998). It has also been suggested for use to generalize the TPSS functional Tao et al. (2003) by Tao Tao (2005). Recently Bates and Furche Bates and Furche (2012) 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 Dunning (1989); Kendall et al. (1992). All DFT calculations have been performed using the London Tellgren, Soncini, and Helgaker (2008) program. This code utilizes the XCFun library Ekström et al. (2010) 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. Staroverov et al., 2003b. 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) Perdew, Burke, and Ernzerhof (1996), and Keal-Tozer-3 (KT3) Keal and Tozer (2004) 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 1. Errors for these quantities are presented relative to the benchmark data of Ref. Teale et al., 2013. Results are also compared with those including corrections from the Tao-Perdew parameterization Tao and Perdew (2005) of the Vignale-Rasolt-Geldart functional Vignale and Rasolt (1988), taken from Ref. Tellgren et al., 2014. In the strong field regime we consider the prediction of perpendicular paramagnetic bonding Lange et al. (2012) in a field strength of 1 a.u. perpendicular to the internuclear axes of H, He, HeNe and Ne. These non-perturbative calculations are assessed against CCSD(T) results computed using the implementation of Ref. Stopkowicz et al., 2015 in the London program.

HF | CO | N | HO | HCN | HOF | LiH |

NH | HCO | CH | CH | AlF | CHF | CH |

FCCH | FCN | HS | HCP | HFCO | HCO | LiF |

NO | OCS | HCO | PN | SO |

## Iv Results

### iv.1 The weak field regime: magnetic properties

We commence by considering the molecular magnetizabilities and NMR shielding constants of the 26 small molecules in Table 1. The magnetizability tensor elements, , are defined as

(14) |

where and label cartesian components of the tensor and magnetic field. The NMR shielding tensor for a given nucleus is defined by

(15) |

where is the nuclear magnetic moment of nucleus . 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. Tellgren et al., 2014, 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.2, 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. Teale et al., 2013. 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. Teale et al., 2013; Lutnaes et al., 2009.

Magnetizability | NMR Shielding | |||||
---|---|---|---|---|---|---|

(10) | (ppm) | |||||

ME | MAE | SD | ME | MAE | SD | |

HF | -3.06 | 6.35 | 7.29 | -15.18 | 21.40 | 40.97 |

LDA | 5.01 | 9.18 | 10.86 | -24.81 | 24.85 | 30.00 |

PBE | 6.75 | 8.81 | 9.03 | -19.66 | 19.78 | 21.46 |

PBE+VRG(TP) | 7.85 | 9.61 | 9.44 | -20.33 | 20.46 | 22.43 |

KT3 | 8.18 | 8.95 | 7.83 | -6.53 | 8.94 | 13.13 |

KT3+VRG(TP) | 9.18 | 9.83 | 8.19 | -7.45 | 9.18 | 13.37 |

B97-2 | 5.46 | 5.84 | 5.96 | -16.34 | 16.48 | 20.60 |

cB98 | 0.52 | 4.84 | 6.58 | -12.44 | 12.66 | 17.79 |

cTPSS | 7.13 | 7.51 | 6.76 | -14.14 | 14.35 | 15.61 |

cTPSS(h) | 6.41 | 6.51 | 6.00 | -14.33 | 14.52 | 16.68 |

The errors in the calculated magnetizabilities and NMR shielding constants are summarised in Table 2 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 2 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. Lutnaes et al., 2009 for these type of functionals. The B97-2 functional gives slightly reduced errors and this is consistent with previous conclusions Lutnaes et al. (2009) 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 Tellgren et al. (2014) 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 Tao and Perdew (2005), 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. Tellgren et al., 2014; Zhu, Zhang, and Trickey, 2014.

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 Bates and Furche (2012). 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. Reiman et al., 2015. We will now examine how these functionals perform when the magnetic field becomes much higher and has a stronger effect on the electronic structure.

### iv.2 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.Lange et al. (2012) 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. Murdin et al. (2013) have shown that phosphorus and selenium doped silicon semiconductors can produce a viable laboratory analogue of free hydrogen Murdin et al. (2013) and helium Litvinenko et al. (2014) 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, He, NeHe and Ne. In particular, we consider the state of H and the lowest states of He, HeNe and Ne. 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 of H 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 alStopkowicz et al. (2015) 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, He, NeHe and Ne are shown in Figures 3, 4, 5 and 6, respectively. Equilibrium bond lengths and dissociation energies were determined numerically and are presented in Table 3.

/ | / mE | |||||||
---|---|---|---|---|---|---|---|---|

H | He | NeHe | Ne | H | He | NeHe | Ne | |

HF | 2.708 | 3.296 | 3.543 | 3.773 | 2.340 | 0.218 | 0.482 | 0.978 |

LDA | 2.374 | 2.550 | 2.846 | 3.062 | 14.81 | 8.231 | 11.817 | 19.730 |

PBE | 2.514 | 2.810 | 3.080 | 3.294 | 5.985 | 2.463 | 3.938 | 7.128 |

KT3 | 2.511 | 2.852 | 3.121 | 3.342 | 4.637 | 2.661 | 3.937 | 6.527 |

cB98 | 2.640 | 3.203 | 3.472 | 3.676 | 1.328 | 1.011 | 1.514 | 2.050 |

cTPSS | 2.564 | 2.864 | 3.124 | 3.346 | 5.255 | 1.307 | 2.344 | 4.577 |

cTPSSh | 2.558 | 2.879 | 3.154 | 3.379 | 5.263 | 1.245 | 2.130 | 3.990 |

CCSD(T) | 2.578 | 2.977 | 3.248 | 3.487 | 4.554 | 1.259 | 2.217 | 4.016 |

FCI | 2.578 | 2.975 | - | - | 4.554 | 1.271 | - | - |

The He FCI calculations use the aug-cc-pVTZ basis set.

The state of the H molecule in a perpendicular field was examined at the FCI level by Lange et al. Lange et al. (2012). 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 3. 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 in Figure 4 we see that HF tends to under-bind with a bond length of compared with the CCSD(T) value of . Similarly the HF dissociation energy of mE is much smaller than the corresponding CCSD(T) value of mE. 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 values and much too large estimates of . 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 and 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. Becke, 1998 to define the cB98 form. These parameters were determined at zero field and from post-LDA calculations – as 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 with those for the VRG-based estimates in Figure 7 of Ref. Tellgren et al., 2014. 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 impressive – with cTPSS(h) being consistently slightly more accurate than cTPSS. This trend is reflected in both the potential energy curves and the and values in Table 3.

## 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 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

(16) |

Here we note that the occupied KS orbitals 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 and the basis functions . We will commence by considering how the molecular orbital energies associated with H and the rare gas dimers change upon application of a magnetic field as the perpendicular paramagnetic bonding in Section V.1 evolves. Since the orbitals themselves can be complex in the presence of a field we then proceed in Section V.2 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.

### v.1 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) Perdew et al. (1982); Almbladh and Von Barth (1985) and the remaining orbital energies can be interpreted as Koopman’s type approximations to higher IPs Bartlett et al. (2005). The extent to which these properties hold for general practical approximations has been a subject of debate in the literature Baerends and Gritsenko (1997); Stowasser and Hoffmann (1999), as has the interpretation of KS virtual orbitals van Meer, Gritsenko, and Baerends (2014) owing to the role of the integer discontinuity Perdew et al. (1982), 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 and He, respectively. For the H molecule in the state we consider the energy of the occupied and 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 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 orbital is stabilized by mE, whilst the singly occupied is destabilized by mE 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 orbital is stabilized by mE, whilst the orbital is destabilized by mE. This greater stabilization of the orbital leads to a net bonding interaction, consistent with the analysis of Ref. Lange et al., 2012. 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 state of He 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 equilibrium internuclear separation a.u. in the presence of a perpendicular field. In the absence of a field we see that again the relative stabilization of and destabilization of are approximately compensatory, whilst in the presence of a field the orbital is more stabilized than the orbital is destabilized. It is also clear that the extra stabilization of the orbital of mE is considerably less than the mE observed for H. This is consistent with the strength of binding exhibited for these species in Figures 3 and 4 as well as in Table 3.

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.

### v.2 Electron density analysis

For each of the species He, HeNe and Ne, we have performed calculations of the electronic density and current density in the presence of varying perpendicular magnetic fields () using the cTPSS functional at the corresponding equilibrium geometries. Here we consider the density change for each system relative to the isolated atoms in the same field,

(17) |

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

(18) |

In Figure 9 we present plots of the density differences in Eqs. (17) and (18) for each of the species with 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 circulations 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 detailedDobson (1992, 1993); Tao (2005) 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 Tellgren, Soncini, and Helgaker (2008). 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 density-functionals. In contrast to vorticity dependent forms previously studied Tellgren et al. (2014); Zhu, Zhang, and Trickey (2014) 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 Goll et al. (2009); Peverati and Truhlar (2011) 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 field causes compression of the electronic density in two dimensions perpendicular to the field and elongation along it. As a result non-uniform coordinate scaling relations may be one example of conditions that may provide a powerful tool in further development. The London program Tellgren, Soncini, and Helgaker (2008) provides a powerful platform for this work since CDFT approaches can be calibrated against accurate ab initio data for a range of field strengths, where experimental data may be scarce. In the future we hope that these relatively inexpensive CDFT approaches may then be applied to the study of larger systems such as those in Refs. Murdin et al., 2013; Litvinenko et al., 2014.

## Appendix A meta-GGAs in CDFT

We present the key working equations defining the B98 and TPSS functionals, indicating how the the modified of Eq. (13) enters each functional.

### a.1 cB98

The B98 functional Becke (1998) has a general construction that is similar to the popular B97 functional form Becke (1997), however, instead of using reduced spin-density gradients it makes use of an inhomogeneity parameter

(19) |

where the exchange hole curvature for the uniform electron gas takes the simple form

(20) |

and is defined in Eq. (10). This inhomogeneity factor, , controls the enhancement or attenuation of the exchange and correlation energy over the uniform gas values.

The exchange component of the functional takes the form

(21) |

and the opposite- and same-spin components of the correlation energy take the forms

(22) |

(23) |

The functions , and are dimensionless inhomogeneity correction factors depending on and . In addition the same-spin correlation energy contains a self-correlation correction (SCC) factor

(24) |

This factor varies between and and vanishes in one-orbital regions, ensuring the functional is self-correlation free. For convenience in deriving fitted forms for the inhomogeneity correction factors , Becke proposed the following transformation to a finite interval,

(25) |

where is a parameter to be determined and separate transformations are carried out for the exchange, opposite-spin correlation and like-spin correlation respectively using the appropriate definitions of as in Eqs (21)–(22). Based on calculations of atomic exchange–correlation energies Becke proposed the values , and . The inhomogeneity corrections were then determined by fitting a power series expansion of the form

(26) |

with 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 . 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 .

The original definition of B98 utilized the zero-field exchange hole curvature of Eq. (10) in its definition. However, it was noted Becke (1998) 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.

### a.2 cTPSS

One of the most widely used meta-GGA functionals is due to Tao, Perdew, Staroverov and Scuseria (TPSS)Tao et al. (2003). 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

(27) |

plays a key role as a dimensionless inhomogeneity parameter, along with

(28) |

Note that throughout this work we use the definition of in Eq. (9), which does not include the factor of 1/2 commonly employed. The exchange functional then takes the form

(29) |

where the precise details of the form chosen for the enhancement factor can be found in Ref. 26. The correlation energy takes the form

(30) |

where hartree.

###### Acknowledgements.

A. M. T. is grateful for support from the Royal Society University Research Fellowship scheme. We are grateful for access to the University of Nottingham High Performance Computing Facility. This work was supported by the Norwegian Research Council through the CoE Centre for Theoretical and Computational Chemistry (CTCC) Grant No. 179568/V30 and the Grant No. 171185/V30 and through the European Research Council under the European Union Seventh Framework Program through the Advanced Grant ABACUS, ERC Grant Agreement No. 267683.## References

- Vignale and Rasolt (1987) G. Vignale and M. Rasolt, Phys. Rev. Lett. 59, 2360 (1987).
- Vignale and Rasolt (1988) G. Vignale and M. Rasolt, Phys. Rev. B 37, 10685 (1988).
- Vignale, Rasolt, and Geldart (1988) G. Vignale, M. Rasolt, and D. J. W. Geldart, Phys. Rev. B 37, 1988 (1988).
- Pan and Sahni (2010) X. Pan and V. Sahni, Int. J. Quantum Chem. 110, 2833 (2010).
- Pan and Sahni (2012) X.-Y. Pan and V. Sahni, J. Phys. Chem. Solids 73, 630 (2012).
- Pan and Sahni (2013) X.-Y. Pan and V. Sahni, Int. J. Quantum Chem. 113, 1424 (2013).
- Tellgren et al. (2012) E. I. Tellgren, S. Kvaal, E. Sagvolden, U. Ekström, A. M. Teale, and T. Helgaker, Phys. Rev. A 86, 062506 (2012).
- Vignale, Ullrich, and Capelle (2013) G. Vignale, C. a. Ullrich, and K. Capelle, Int. J. Quantum Chem. 113, 1422 (2013).
- Lieb (1983) E. H. Lieb, Int. J. Quantum Chem. 24, 243 (1983).
- Lieb and Schrader (2013) E. H. Lieb and R. Schrader, Phys. Rev. A 88, 032516 (2013).
- Tellgren, Kvaal, and Helgaker (2014) E. I. Tellgren, S. Kvaal, and T. Helgaker, Phys. Rev. A 89, 012515 (2014).
- Lee, Handy, and Colwell (1995) A. M. Lee, N. C. Handy, and S. M. Colwell, J. Chem. Phys. 103, 10095 (1995).
- Bates and Furche (2012) J. E. Bates and F. Furche, J. Chem. Phys. 137, 164105 (2012).
- Pittalis et al. (2006) S. Pittalis, S. Kurth, N. Helbig, and E. K. U. Gross, Phys. Rev. A 74, 062511 (2006).
- Pittalis et al. (2007) S. Pittalis, S. Kurth, S. Sharma, and E. K. U. Gross, J. Chem. Phys. 127, 124103 (2007).
- Zhu, Zhang, and Trickey (2014) W. Zhu, L. Zhang, and S. B. Trickey, Phys. Rev. A 90, 022504 (2014).
- Tellgren et al. (2014) E. I. Tellgren, A. M. Teale, J. W. Furness, K. K. Lange, U. Ekström, and T. Helgaker, J. Chem. Phys. 140, 034101 (2014).
- Tellgren, Soncini, and Helgaker (2008) E. I. Tellgren, A. Soncini, and T. Helgaker, J. Chem. Phys. 129, 154114 (2008).
- London (1937) F. London, J. Phys. Rad. 8, 397 (1937).
- Ditchfield (1974) R. Ditchfield, Mol. Phys. 27, 789 (1974).
- Ditchfield (1976) R. Ditchfield, J. Chem. Phys. 65, 3123 (1976).
- Maximoff and Scuseria (2004) S. N. Maximoff and G. E. Scuseria, Chem. Phys. Lett. 390, 408 (2004).
- Tao (2005) J. Tao, Phys. Rev. B 71, 205107 (2005).
- Dobson (1993) J. F. Dobson, J. Chem. Phys. 98, 8870 (1993).
- Becke (1996) A. D. Becke, Can. J. Chem. 74, 995 (1996).
- Tao et al. (2003) J. Tao, J. Perdew, V. Staroverov, and G. Scuseria, Phys. Rev. Lett. 91, 146401 (2003).
- Becke (1998) A. D. Becke, J. Chem. Phys. 109, 2092 (1998).
- Staroverov et al. (2003a) V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. Perdew, J. Chem. Phys. 119, 12129 (2003a).
- Stopkowicz et al. (2015) S. Stopkowicz, J. Gauss, K. K. Lange, E. I. Tellgren, and T. Helgaker, submitted (2015).
- Orestes, Marcasso, and Capelle (2003) E. Orestes, T. Marcasso, and K. Capelle, Phys. Rev. A 68, 022105 (2003).
- Tao and Perdew (2005) J. Tao and J. Perdew, Phys. Rev. Lett. 95, 196403 (2005).
- Tao and Vignale (2006) J. Tao and G. Vignale, Phys. Rev. B 74, 193108 (2006).
- Higuchi et al. (2007) K. Higuchi, M. Miyasita, M. Kodera, and M. Higuchi, J. Magn. Magn. Mater. 310, 1065 (2007).
- Dobson (1992) J. F. Dobson, J. Phys. Condens. Matter 7877 (1992).
- Becke (1983) A. D. Becke, Int. J. Quantum Chem. 23, 1915 (1983).
- Becke and Roussel (1989) A. D. Becke and M. R. Roussel, Phys. Rev. A 39, 3761 (1989).
- Dunning (1989) T. H. Dunning, J. Chem. Phys. 90, 1007 (1989).
- Kendall et al. (1992) R. Kendall, R. Kendall, T. Jr., Dunning, T. Jr., Dunning, R. Harrison, and R. Harrison, J. Chem. Phys. 96, 6796 (1992).
- Ekström et al. (2010) U. Ekström, L. Visscher, R. Bast, and A. J. Thorvaldsen, J. Chem. Theory Comput. 6, 1971 (2010).
- Staroverov et al. (2003b) V. N. Staroverov, G. E. Scuseria, J. Tao, and J. P. Perdew, J. Chem. Phys. 119, 12129 (2003b).
- Perdew, Burke, and Ernzerhof (1996) J. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996).
- Keal and Tozer (2004) T. W. Keal and D. J. Tozer, J. Chem. Phys. 121, 5654 (2004).
- Teale et al. (2013) A. M. Teale, O. B. Lutnæ s, T. Helgaker, D. J. Tozer, and J. Gauss, J. Chem. Phys. 138, 024111 (2013).
- Lange et al. (2012) K. K. Lange, E. I. Tellgren, M. R. Hoffmann, and T. Helgaker, Science 337, 327 (2012).
- Lutnaes et al. (2009) O. B. Lutnaes, A. M. Teale, T. Helgaker, D. J. Tozer, K. Ruud, and J. Gauss, J. Chem. Phys. 131, 144104 (2009).
- Reiman et al. (2015) S. Reiman, U. Ekström, S. Stopkowicz, A. M. Teale, A. Borgoo, and T. Helgaker, Phys. Chem. Chem. Phys. (submitted) (2015).
- Murdin et al. (2013) B. N. Murdin, J. Li, M. L. Y. Pang, E. T. Bowyer, K. L. Litvinenko, S. K. Clowes, H. Engelkamp, C. R. Pidgeon, I. Galbraith, N. V. Abrosimov, H. Riemann, S. G. Pavlov, H.-W. Hübers, and P. G. Murdin, Nat. Commun. 4, 1469 (2013).
- Litvinenko et al. (2014) K. L. Litvinenko, M. Pang, J. Li, E. Bowyer, H. Engelkamp, V. B. Shuman, L. M. Portsel, a. N. Lodygin, Y. a. Astrov, S. G. Pavlov, H.-W. Hübers, C. R. Pidgeon, and B. N. Murdin, Phys. Rev. B 90, 1 (2014).
- Perdew et al. (1982) J. P. Perdew, R. G. Parr, M. Levy, and J. L. Balduz, Phys. Rev. Lett. 49, 1691 (1982).
- Almbladh and Von Barth (1985) C. O. Almbladh and U. Von Barth, Phys. Rev. B 31, 3231 (1985).
- Bartlett et al. (2005) R. J. Bartlett, I. Grabowski, S. Hirata, and S. Ivanov, J. Chem. Phys. 122, 034104 (2005).
- Baerends and Gritsenko (1997) E. J. Baerends and O. V. Gritsenko, J. Phys. Chem. A 101, 5383 (1997).
- Stowasser and Hoffmann (1999) R. Stowasser and R. Hoffmann, J. Am. Chem. Soc. 121, 3414 (1999).
- van Meer, Gritsenko, and Baerends (2014) R. van Meer, O. V. Gritsenko, and E. J. Baerends, J. Chem. Theory Comput. 10, 4432 (2014).
- Goll et al. (2009) E. Goll, M. Ernst, F. Moegle-Hofacker, and H. Stoll, J. Chem. Phys. 130, 234112 (2009).
- Peverati and Truhlar (2011) R. Peverati and D. G. Truhlar, J. Phys. Chem. Lett. 2, 2810 (2011).
- Becke (1997) A. D. Becke, J. Chem. Phys. 107, 8554 (1997).