PHYSICAL REVIEW A 90, 022504 (2014) Comparative studies of density-functional approximations for light atoms in strong magnetic fields Wuming Zhu Department of Physics, Hangzhou Normal University, 16 Xuelin Street, Hangzhou, Zhejiang 310036, China Liang Zhang School of Information Science and Engineering, Hangzhou Normal University, 16 Xuelin Street, Hangzhou, Zhejiang 310036, China S. B. Trickey Quantum Theory Project, Department of Physics and Department of Chemistry, P.O. Box 118435, University of Florida, Gainesville, Florida 32611-8435, USA (Received 30 May 2014; published 11 August 2014) For a wide range of magnetic fields, 0 B 2000 a.u., we present a systematic comparative study of the performance of different types of density-functional approximations in light atoms (2 Z 6). Local, generalized-gradient approximation (GGA; semilocal), and meta-GGA ground-state exchange-correlation (xc) functionals are compared on an equal footing with exact-exchange, Hartree-Fock (HF), and current-densityfunctional-theory (CDFT) approximations. Comparison also is made with published quantum Monte Carlo data. Though all approximations give qualitatively reasonable results, the exchange energies from local and GGA functionals are too negative for large B. Results from the Perdew-Burke-Ernzerhof ground-state GGA and Tao-Perdew-Staroverov-Scuseria (TPSS) ground-state meta-GGA functionals are very close. Because of confinement, self-interaction error in such functionals is more severe at large B than at B = 0, hence selfinteraction correction is crucial. Exact exchange combined with the TPSS correlation functional results in a self-interaction-free (xc) functional, from which we obtain atomic energies of comparable accuracy to those from correlated wave-function methods. Specifically for the B and C atoms, we provide beyond-HF energies in a wide range of B fields. Fully self-consistent CDFT calculations were done with the Vignale-Rasolt-Geldart (VRG) functional in conjunction with the PW92 xc functional. Current effects turn out to be small, and the vorticity variable in the VRG functional diverges in some low-density regions. This part of the study suggests that nonlocal, self-interaction-free functionals may be better than local approximations as a starting point for CDFT functional construction and that some basic variable other than the vorticity could be helpful in making CDFT calculations practical. DOI: 10.1103/PhysRevA.90.022504 PACS number(s): 31.15.E− I. INTRODUCTION Density-functional theory (DFT) [1] has been enormously successful in electronic structure calculations of atoms, molecules, surfaces, and solids. The success was fostered by the availability of relatively accurate yet simple approximations to the exchange-correlation (xc) energy functional (and the corresponding potential in the Kohn-Sham equation). Examples are the local density and local spin-density approximations (LDA, LSDA) and generalized-gradient approximations (GGAs). For atoms and molecules, the literature of numerical studies comparing the accuracy of such approximations to that of wave-function or many-body-perturbation methods (socalled ab initio methods) is too large to cite with any fairness. For atoms subjected to intense magnetic fields (B > 1 hartree a.u. = 2.3505 × 105 T), however, comparatively little is known about the behavior and relative quality of various DFT approximations [2–5]. What is known for physical atoms is limited (as far as we can ascertain) to the LDA. Considerably more investigation has been done on quantum dots, but, as systems confined to two dimensions, they are not of direct relevance. Some time ago, two of us showed that various approximate functionals in DFT and current DFT (CDFT) [6] generate substantial discrepancies compared to exact results for the Hooke’s atom in an external field [7]. Since the density functional is, in principle, universal, those discrepancies are not 1050-2947/2014/90(2)/022504(14) encouraging, hence a critical issue is the extent to which such discrepancies occur in real three-dimensional (3D) atoms. That issue was sharpened by the subsequent discovery of peculiar v-representability properties for the two-electron system [8,9]. The present paper addresses the issue directly by systematic calculation of the energetics of light atoms (2 Z 6, Z = atomic number) over a very large range of B (0 B 2000 a.u.) with extremely rich, aspherical basis sets and several xc functionals, including one combination which apparently has gone unexplored heretofore. We take advantage of the availability of several Hartree-Fock (HF) [10–18] and correlated wave-function studies of at least a few light atoms in an external B field [19–27] and of quantum Monte Carlo (QMC) studies [28,29] for comparisons. In particular, comparison of our HF calculations with prior ones provides calibration of the accuracy of our techniques, while comparison with the QMC results gives a measure of absolute accuracy and trends in accuracy for the various DFT approximations. Though the manifestly gauge-invariant current CDFT was established more than two decades ago [6], CDFT calculations on 3D systems remain rare, with the majority being for vanishing or low-B-field cases [30–33]. There has been recent work on molecular magnetic properties in CDFT up to about B = 1 a.u. by Tellgren et al. [34]. In several ways their work is complementary to ours, thus some of their findings and technical details are discussed below. 022504-1 ©2014 American Physical Society WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) CDFT parametrizes the many-electron ground state via the electron number density n(r) and paramagnetic current density jp (r). The current-independent limit is ordinary DFT with appropriate explicit dependence on B, which is the second form of DFT studied here. Because it is a simple extension of existing zero-field DFT approximations to B = 0, that second form enables better comprehension of the challenges in CDFT calculations. As suggested above, some of the approximate xc functionals are well tested in the B = 0 case [35–37], but their performance for B = 0 is not. We did not include clusters and solids for several reasons: High-level correlated results for such systems are rare, there is the technical burden of developing accurate, fast multicenter integrals for the special basis functions suitable for high B fields, and the study of atoms avoids issues of gauge invariance and periodicity breaking. The essentials of the DFT and CDFT formulations are in the next section. Basis sets and numerical methods are addressed in Sec. III (and the Appendix), followed by results in Sec. IV, and a concluding summary in Sec. V. The xc contributions are + 1 1 . 2 i=j |ri − rj | (1) Here ri , mi , and ms,i are the space coordinate, magnetic quantum number, and spin z component for the ith electron, respectively. (As indicated already, Hartree atomic units = melectron = qelectron = 1 are used throughout unless stated explicitly; one Hartree a.u. of B field corresponds to 2.3505 × 105 T.) Without explicit exception, unpaired electrons always are taken as spin down in what follows, i.e., ms,i = − 12 for unpaired electrons. In CDFT, the variational minimization for the manyelectron ground state is mapped to a Kohn-Sham (KS) system which generates n(r) and jp (r). The system xc energy thus depends on both in general. The CDFT KS equation is B2 2 1 Z B (x + y 2 ) + (mi + 2ms,i ) + vH (r) − ∇2 − + 2 r 8 2 1 + vxc (r) + Axc (r) · ∇ φi (r) = i φi (r), (2) i with vH (r) = n(r ) dr . |r − r | δExc [n(r),jp (r)] δn(r) (4) Axc (r) = δExc [n(r),jp (r)] , δjp (r) (5) and while the densities are n(r) = |φi (r)|2 (6) i and jp (r) = 1 ∗ [φi (r)∇φi (r) − φi (r)∇φi∗ (r)]. 2i i The system total energy is Etot = Ts + J + Exc [n,jp ] + dr n(r) − (7) Z r B B2 2 2 (x + y ) + (mi + 2ms,i ), + 8 2 i II. METHODOLOGY Upon imposition of a uniform, static external magnetic field along the z direction, B = B zˆ , a central field atom at the origin becomes cylindrically symmetric. The system Hamiltonian commutes with rotations about the field direction, so the magnetic quantum number m is still good. Three other conserved quantities are the total spin S2 , its z component Sz , and spatial parity in z, hence a state can be labeled by those quantum numbers. In the Coulomb gauge, A(r) = 12 B × r, with A the vector potential, the system Hamiltonian then is 1 B Z B2 2 2 2 xi + yi + (mi + 2ms,i ) H= − ∇i − + 2 ri 8 2 i vxc (r) = (8) where Ts denotes the noninteracting (KS) kinetic energy and J is the Hartree (classical electron-electron repulsion) energy. In the spirit of linear response, it often is assumed that the Axc (r) contribution to Etot is small compared to the ordinary DFT Exc . The zeroth-order approximation to Exc [n,jp ] thus takes the same form as the xc functional in ordinary DFT, Exc [n], and the KS equation reduces to the familiar form except that explicit external field contributions appear. We refer to this approximation as naive B-DFT (nB-DFT). Observe that in the nB-DFT approximation, only the explicit dependence of Exc [n,jp ] upon jp (r) is neglected, that is, Exc [n,jp ] ≈ Exc [n,0]. Other external field dependencies, including the interaction between jp and the external B field, and the implicit dependence on jp (r), are unaltered. The diversity of approximate xc functionals for the B = 0 case is well documented. Following the Perdew-Schmidt Jacob’s ladder of functional complexity [38], the first two rungs, LDA and GGA, are local and semilocal functionals. The third, meta-GGA (mGGA), has an explicit orbital dependence from the KS kinetic energy density. We choose the PW92 [39], Perdew-Burke-Ernzerhof (PBE) [40], and Tao-Perdew-Staroverov-Scuseria (TPSS) [41] functionals as representatives of the respective rungs. As is well known, both local and semilocal functionals suffer from improper electron self-interaction. Among the correlation functionals for those three, only the TPSS correlation functional is self-interaction free [41,42]. Therefore, the combination of exact-exchange (EXX) and the TPSS correlation functional, Exc = ExXX + EcTPSS , (9) is a self-interaction-free xc functional. Exact exchange in DFT [43] is defined in terms of occupied KS orbitals {φi (r)} in the Fock integral, ∗ φi (r)φj∗ (r )φj (r)φi (r ) 1 XX dr dr . Ex [{φi (r)}] = − 2 i,j |r − r | (3) 022504-2 (10) COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . . PHYSICAL REVIEW A 90, 022504 (2014) Since ExXX is fully nonlocal, the combination in Eq. (9) corresponds to the fourth rung of the Jacob’s ladder or hyper-GGA (HGGA). Note that EcTPSS is not fully compatible with ExXX , since the TPSS correlation hole is rather local and misses the long-range part. Because the spatial extension of atomic systems usually is limited such that long-range hole components are negligible or nearly so [44], this combination should not be a problem for the present study. That appraisal is confirmed numerically (see below). However, one cannot assume that Eq. (9) also would be generally applicable to extended systems, such as molecules or solids. For a correlation functional fully compatible with exact exchange, we would need to move up to the fifth rung of the Jacob’s ladder. An example is the generalized random phase approximation (RPA) [45]. However, such RPA calculations would involve significant additional development, a task we leave to a future study. The implicit dependence of Exc on jp in the nB-DFT approximation does not matter for local and semilocal functionals, but does for orbital-dependent functionals. Equation (10) is a specific case. Since KS orbitals depend on jp , and ExXX in turn depends on those KS orbitals, there is an indirect dependence of ExXX on jp . Thus, the current contribution to the exchange energy is included at least partially in the nB-DFT approximation with exact exchange [see Eq. (9)]. An additional complication is that to stay strictly within the KS scheme therefore would require solution of the corresponding optimized effective potential (OEP) equations [46]. The OEP method and associated simplifying approximations have been generalized to the current-spin-DFT (CSDFT) case [47]. However, the numerical burden of the OEP equations is substantial, with the result that direct variation of Eq. (9) with respect to the orbitals frequently is used instead. Called “generalized KS” in the quantum chemistry literature, this approach gives one-electron equations of the same structure as the HF equations but with a local correlation potential added. Numerical results show that atomic energies obtained from OEP and generalized KS calculations are quite close for the X-only functional [46,48]. Thus we expect the energies also to be close for OEP and generalized KS calculations when Eq. (9) is used, because its correlation part is a mGGA. For computational simplicity, we therefore used the generalized KS procedure for orbital-dependent functionals. In contrast to the nB-DFT situation, there is little choice in xc functionals for CDFT. To our knowledge, the only one published for 3D systems that is in a form implementable for self-consistent solution of the KS equations is the local approximation due to Vignale, Rasolt, and Geldart (VRG) [6]. From the perturbative energy of a homogeneous electron gas (HEG) in uniform B and gauge-invariance arguments, Vignale and Rasolt gave the local approximation for the current contribution to the xc energy [6], VRG Exc [n(r),ν(r)] = g(r) |ν(r)|2 dr, (11) where jp (r) , n(r) χ (n(r)) kF − 1 . g(r) = g(n(r)) = 24π 2 χ0 (n(r)) ν(r) ≡ ∇ × (12) (13) The variable ν is the vorticity mentioned already, while χ and χ0 are the orbital magnetic susceptibilities for the interacting and noninteracting HEG, respectively. kF = (3π 2 n)1/3 is the usual local Fermi momentum. From computed susceptibility data over the range 1 rs 10 tabulated in Ref. [49], Lee, Colwell, and Handy (LCH) generated an analytical fit for the ratio [30], χ , (14) s≡ χ0 sLCH = (1.0 + 0.028rs )e−0.042rs (15) 3 1/3 with rs = ( 4πn ) the usual Wigner-Seitz radius. The resulting approximation is gLCH = kF (sLCH − 1) . 24π 2 (16) With various regularizations (cutoffs) which we discuss shortly, the VRG functional has been used to calculate magnetizabilities [30,34], nuclear shielding constants [31], and frequency-dependent polarizabilities [32] for small molecules, and ionization energies for atoms [33]. Because all those studies are for B → 0, they are inconclusive regarding the general, B 0 utility of the VRG functional. Except for our Hooke’s atom results, Ref. [7], little is known about the field-dependent behavior of finite systems for very large B that emerges from the VRG functional. An important technical issue is evident in Eq. (15), namely, the sensitivity of the predicted physics to the details of the low-density extrapolation. Orestes, Marcasso, and Capelle (OMC) proposed two other fits, both polynomial [33]. A fit by Tao and Perdew [50] actually gives the correct rs → 0 limit and physically plausible low-density behavior. Subsequently a refinement was given by Tao and Vignale [51]. As discussed in detail in Ref. [34], those two seem to differ little in practice. That reference also shows how the LCH, Tao-Perdew, and Tao-Vignale parametrizations tend to zero as rs → ∞, whereas the OMC and VRG forms do not. They also studied the consequences of the peculiar limiting behavior of the rather different Higuchi and Higuchi [52] functional. The relevant point here is that even at the comparatively low fields considered in Ref. [34], there is considerable variation in the numerical stability for that collection of low-density regularizations but not large shifts in the computed quantities. Thus we resorted to work we had done earlier on a different regularization [7]. From a study of many numerical experiments, we arrived at gcutoff = kF (c1 + c2 rs )e−αcutoff rs 24π 2 (17) for densities below ncutoff . The values of c1 and c2 follow from requiring a smooth connection between gcutoff and g(n) at the specified transition density ncutoff . In this work, we use ncutoff = 0.001 bohr−3 and αcutoff = 2.0 bohrs−1 , unless otherwise explicitly specified. The numerical value of ncutoff corresponds to rs = 6.2 bohrs, so that the regularization afforded by gcutoff could be tested and calibrated against the available data on 6.2 rs 10 bohrs. 022504-3 WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) where III. BASIS SETS AND NUMERICAL METHODS In molecular calculations, the KS equations usually are solved by expansion in a Gaussian-type orbital (GTO) basis. Less commonly, Slater-type orbitals are used, but they are not well suited for systems in a strong B field. Specifically Jones, Ortiz, and Ceperley used GTOs to calculate small atoms and molecules in strong fields, and found that basis elements with very large angular momentum (lmax up to 35) are required for reasonable numerical convergence [15,53]. Tellgren et al. [34] used rich, conventionally isotropic GTO basis sets and commented that their calculations thus were limited to maximum B ≈ 1 a.u. For vanishing B field, Porezag and Pederson introduced optimized GTO basis sets for DFT calculations [54]. The difficulty for our purposes is that their technique is based on total energy optimization and therefore requires specification of the theoretical framework (HF or DFT). Our intent is to compare the consequences of different approximations and methods with the same basis sets for all calculations in order to avoid artifacts. Here we choose anisotropic-GTO (AGTO) basis sets such as were introduced by Aldrich and Greene [55], and exploited extensively by Schmelcher and Cederbaum [56]. In cylindrical coordinates (ρ,z,φ), the basis functions are χj (ρ,z,φ) = Nj ρ nρj znzj e−αj ρ 2 −βj z2 imj φ e , where nρj = |mj | + 2kj , kj = 0,1, . . . , with mj = . . . , −2, −1,0,1,2, . . . , and nzj = πzj + 2lj ,lj = 0,1, . . . , with πzj = 0,1. AGTO basis sets have been used in several studies of atoms in strong B fields [20,21,24,25,27] because of their flexibility in describing elongation of the electron orbitals and densities along the B-field direction. With no unique way to determine the basis exponents, αj and βj , a physically sound prescription must be developed. Kappes and Schmelcher gave an optimization algorithm [57] which, however, must be executed for each combination of atomic configuration and field strength, a tedious, time-consuming task. Furthermore, a good initial guess is required [21]. An alternative is the use of nearly optimized basis sets. Kravchenko and Liberman (KL) [58] investigated one-electron systems, the hydrogen atom and the hydrogen molecular ion, and showed that systematically constructed AGTO basis sets could provide accuracy of 10−6 hartree or better. As detailed in the Appendix, we constructed KL-like highly optimized basis sets without needing full nonlinear optimization. Each is comprised of a sequence of subsets related by αj, = βj + (1 + μ )(B), μ2,3 = ±0.2, and, in practice, γ = B/Z 2 . The base sequence is μ=1 = 0. Others are μ = ±0.2,±0.4 for the second, third, fourth, and fifth sequences, respectively. For = 2,3, there are half as many exponents (with doubled spacing) as in the base sequence, while for = 4,5 there are one-fourth as many (with quadrupled spacing) as in the base sequence. Similar to the KL basis sets, we used even-tempered Gaussian (ETG) sequences for the longitudinal exponents βj . Again, see the Appendix. The kinetic energy, nuclear-attraction, overlap, and onecenter Coulomb repulsion matrix elements with respect to the AGTO basis all can be expressed in closed forms, though the Coulomb repulsion expression is lengthy [21]. The vorticity, ν(r), was calculated analytically. For CDFT, both xc potentials, vxc (r) and Axc (r), were evaluated numerically on a twodimensional mesh, with up to 500 points along each direction. The convergence of total energies with respect to basis set size and with respect to the number of mesh points was checked carefully. The self-consistent field (SCF) convergence criterion was 10−8 in total energy. In general the mesh convergence was no worse than a few units in the last digit of the tabulated data. j = 1,2,3, . . . , (18) μ1 = 0.0, b(γ ) = −0.16[tan−1 (γ )]2 + 0.77 tan−1 (γ ) + 0.74, (21) μ4,5 = ±0.4, (19) with αj,1 ≈ βj ∀j and (B) a parameter optimized for each B. The result of extensive numerical exploration (see Appendix) is B 4 βj −1/2 4 βj −2 (B) = + 1+ , 4 1+ 20 b(γ ) B b(γ ) B (20) IV. RESULTS AND DISCUSSION A summary of our results at four B magnitudes (42.54 a.u. = 1.0 × 107 T, 212.72 a.u. = 5.0 × 107 T, 425.43 a.u. = 1.0 × 108 T, 2127.2 a.u. = 5.0 × 108 T) with comparison to a broad selection of near-exact calculations, is given in Tables I and II. For illustration, a selection of detailed results is in Tables III and IV for the He and Li atoms, again with comparison to published data. Those four tables are the main focus of the discussion in this section. More extensive tabulations are in the Supplemental Material [59], which includes more electronic states for the He, Li, Be, B, and C atoms and the positive ions, Li+ , Be+ , and B+ , in field strengths 0 B 2000 a.u. In all cases, the electronic states are labeled according to their corresponding zero-field Hartree-Fock electronic configurations. A. Unrestricted Hartree-Fock calculation with AGTO basis We first did unrestricted Hartree-Fock (UHF) calculations mainly as a check on the basis set construction and for verification of numerical technique. The UHF energies are in excellent agreement with published values. For the He atom, Zhao and Stancil expanded the 1s 2 electron orbitals in a B-spline basis and used quadruple-precision calculations to obtain very accurate HF energies for 0 B 100 a.u. [17]. Our results differ from theirs by no more than a couple of μ-hartree. For other electronic states, we compare with the data by Jones, Ortiz, and Ceperley [20], who also used AGTO basis sets. Frequently our energies are slightly lower, perhaps because of better basis optimization. Comparison data for atoms and ions with Z 3 are from the studies by Ivanov and Schmelcher [16], who used two-dimensional (2D) mesh methods. Differences usually are less than 0.1 m-hartree, and the overall agreement is quite satisfactory. Exceptions are the 1s 2 2s 2 state of the Be atom 022504-4 COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . . PHYSICAL REVIEW A 90, 022504 (2014) TABLE I. Atomic energies for He, Li, Be, B, and C atoms in B fields by different methods (energies in hartree, B in tesla). Present work is “HF, pres.” plus the four xc functionals PW92 (LDA), PBE (GGA), TPSS (MGGA), and exact exchange plus TPSS correlation (HGGA). Comparison data are for various QMC methods denoted as in Refs. [28] Tables I–IV for “RPDQMC” through “MCPH3” and [29] for FPDQMC2. Z B (T) HF, pres. LDA GGA MGGA HGGA RPDQMC FPDQMC VQMC HFFEM 2DHF MCPH3 FPDQMC2 2 1.0 × 107 −9.6969 −9.6304 −9.8332 −9.7988 −9.7216 −9.735 −9.735 −9.617 −9.393 −9.697 −9.603 5.0 × 107 −16.928 −17.151 −17.618 −17.535 −16.957 −17.00 −17.01 −16.88 −16.72 −16.93 −16.78 1.0 × 108 −21.314 −21.873 −22.544 −22.425 −21.344 −21.41 −21.41 −21.28 −21.15 −21.31 −21.19 5.0 × 108 −35.349 −37.790 −39.345 −39.074 −35.381 −35.51 −35.54 −35.35 −35.24a −35.35 −35.18 −9.714 −16.968 −21.374 −35.489 3 1.0 × 107 5.0 × 107 1.0 × 108 5. × 108 −19.61 −35.01 −44.61 −76.36 −19.891 −35.420 −45.107 −77.054 4 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 −33.013 −32.916 −33.239 −33.126 −33.076 −33.18 −33.15 −32.50 −31.14 −33.01 −32.70 −59.395 −59.860 −60.571 −60.354 −59.471 −59.61 −59.57 −59.02 −58.14 −59.40 −58.80 −76.184 −77.301 −78.313 −78.036 −76.264 −76.47 −76.44 −75.89 −75.11 −76.18 −75.56 −132.68 −137.56 −139.86 −139.37 −132.76 −133.40 −133.25 −132.55 −131.97 −132.68 −131.78 −33.059 −59.504 −76.344 −133.067 5 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 −48.961 −48.853 −49.240 −49.116 −49.044 −49.17 −49.17 −47.96 −45.53 −48.96 −48.77 −88.602 −89.189 −90.015 −89.765 −88.702 −88.97 −88.86 −87.94 −86.32 −88.60 −87.83 −114.26 −115.65 −116.82 −116.50 −114.37 −114.73 −114.62 −113.74 −112.34 −114.26 −113.37 −202.10 −208.14 −210.78 −210.25 −202.22 −203.04 −203.00 −201.79 −200.83 −202.10 −201.24 −49.021 −88.744 −114.467 −202.605 6 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 −67.580 −67.467 −67.921 −67.796 −67.684 −67.88 −67.95 −65.85 −62.00 −67.58 −68.13 −122.61 −123.32 −124.26 −124.00 −122.73 −123.07 −122.93 −121.49 −118.92 −122.61 −121.57 −158.75 −160.42 −161.74 −161.40 −158.88 −159.42 −159.16 −157.80 −155.67 −158.75 −157.58 −284.26 −291.45 −294.41 −293.86 −284.41 −285.39 −285.32 −283.96 −282.20 −284.26 −282.78 −67.655 −122.783 −159.010 −284.874 a −19.860 −35.345 −44.997 −76.781 −19.777 −35.691 −45.841 −80.471 −20.040 −36.283 −46.688 −82.417 −19.957 −36.122 −46.475 −82.010 −19.904 −35.397 −45.052 −76.840 −19.92 −35.50 −45.20 −77.28 −19.93 −35.47 −45.16 −77.21 −19.58 −35.19 −44.83 −76.73 −18.99 −34.76 −44.50 −76.44 −19.86 −35.35 −45.00 −76.78 An obvious sign error in the original is corrected here. and B+ ion, and the 1s 2 2s 2 2p−1 state of the B atom, especially in very large B fields [59]. Large discrepancies (from 0.5 mhartree at B = 1 a.u. to 1.9 hartree at B = 2000 a.u. for the Be atom, for example) seem plausibly to be attributable to different allowed spatial symmetries for the electron wave function. They adopted asymmetrical wave functions for 2s 2 electrons with respect to the z = 0 plane [16], whereas we required definite z parity, πz = 0,1. Excluding those three instances, we can compare our HF energies with theirs evenhandedly, because they did not find symmetry breaking in any other cases. Usually our HF energies are slightly higher than theirs, with differences no more than a few units in the last quoted digit. Presumably those differences arise from numerical mesh error in their calculations and basis set truncation error in ours. On variational grounds, the basis set truncation error always is positive, so our data are an upper bound for the HF energies. In contrast, there can be small magnitude but negative numerical errors from Ivanov and Schmelcher’s 2D HF mesh method. Compare, for example, their H atom results with those from Ref. [60]. Hence, one may speculate that the true HF energies lie between our data and theirs. The fact that our basis-based calculation is comparable to the best numerical results in accuracy suggests strongly that our construction should result in sufficiently large and flexible basis sets for describing atoms in a wide range of B fields for a range of theoretical methods. B. Results of nB-DFT approximations The PW92, PBE, and TPSS functionals (LDA, GGA, and MGGA, respectively) were implemented and tested in the context of the nB-DFT approximation discussed above. Results again are listed in Tables I, III, and IV, and in the Supplemental Material [59]. Regarding context, there are fewer published DFT calculations for atoms in B than for HF. Comparison of nonzero B-field DFT results moreover is handicapped by the different magnetic field values for which various authors present their results, and by the different xc functionals implemented. The functional due to Jones [2], which is at the level of LDA, was used by Neuhauser, Koonin, and Langanke [11], and by Braun [4]. The simple Dirac X-only functional was used by Relovsky and Ruder [3], and by Braun [4]. Our results for the Dirac X functional agree well with Braun’s. Medin and Lai considered only the high-field case. Based on the adiabatic approximation [5], they chose the Danz-Glasser X functional and a strong-field-limit expression for the C functional [61]. Both are local density approximations which are somewhat crude in comparison to modern ground-state, B = 0 xc functionals. And, it is not clear that a well-behaved low-field limit of that xc functional exists or, if it does, that it recovers the ordinary LDA. As far as we can tell, there is no prior application of density-gradient-dependent functionals to atoms in a strong B field. In the broader context, since DFT calculations in principle include electron correlation, it also is appropriate to compare our nB-DFT results with those from correlated wave-function methods [19–27]. We address some of the more predictable outcomes first. Inspection of Tables I, III, and IV shows that total atomic energies from several different methods differ from a few hundredths of a hartree at low B up to a few hartrees at very large B field. To illustrate the point, Table II recasts the results in 022504-5 WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) TABLE II. Deviations from reference atomic energies (last column, FPDQMC2) for He, Li, Be, B, and C atoms in B fields by different methods (energies in hartree, B in tesla). Present work is “HF, pres.” plus the four xc functionals PW92 (LDA), PBE (GGA), TPSS (MGGA), and exact exchange plus TPSS correlation (HGGA). Comparison data are for various QMC methods denoted as in Refs. [28] Tables I–IV for “RPDQMC” through “VQMC” and [29] for FPDQMC2. Z B (T) HF, pres. LDA GGA MGGA HGGA RPDQMC FPDQMC VQMC 2 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 0.017 0.040 0.060 0.140 0.084 −0.183 −0.499 −2.301 −0.119 −0.650 −1.170 −3.856 −0.085 −0.567 −1.051 −3.585 −0.008 0.011 0.030 0.108 −0.021 −0.032 −0.036 −0.021 −0.021 −0.042 −0.036 −0.051 0.097 0.088 0.094 0.139 −9.714 −16.968 −21.374 −35.489 3 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 0.031 0.075 0.110 0.273 0.114 −0.271 −0.734 −3.417 −0.149 −0.863 −1.581 −5.363 −0.066 −0.702 −1.368 −4.956 −0.013 0.023 0.055 0.214 −0.029 −0.080 −0.093 −0.226 −0.039 −0.050 −0.053 −0.156 0.311 0.230 0.277 0.324 −19.891 −35.420 −45.107 −77.054 4 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 0.046 0.109 0.160 0.387 0.143 −0.356 −0.957 −4.493 −0.180 −1.067 −1.969 −6.793 −0.067 −0.850 −1.692 −6.303 −0.017 0.033 0.080 0.307 −0.121 −0.106 −0.126 −0.333 −0.091 −0.066 −0.096 −0.183 0.559 0.484 0.454 0.517 −33.059 −59.504 −76.344 −133.067 5 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 0.060 0.142 0.207 0.505 0.168 −0.445 −1.183 −5.535 −0.219 −1.271 −2.353 −8.175 −0.095 −1.021 −2.033 −7.645 −0.023 0.042 0.097 0.385 −0.149 −0.226 −0.263 −0.435 −0.149 −0.116 −0.153 −0.395 1.061 0.804 0.727 0.815 −49.021 −88.744 −114.467 −202.605 6 1.0 × 107 5.0 × 107 1.0 × 108 5.0 × 108 0.075 0.173 0.260 0.614 0.188 −0.537 −1.410 −6.576 −0.266 −1.477 −2.730 −9.536 −0.141 −1.217 −2.390 −8.986 −0.029 0.053 0.130 0.464 −0.225 −0.287 −0.410 −0.516 −0.295 −0.147 −0.150 −0.446 1.805 1.293 1.210 0.914 −67.655 −122.783 −159.010 −284.874 Table I in the form of deviations of calculated values for various approximations with respect to the FPDQMC2 values from Ref. [29]. (There are not enough tabulated near-exact results for a similar recasting of Tables III and IV, but study shows similar results when comparison is possible.) Even within DFT there is similar variation for different xc functionals. All give similar, qualitatively correct results. The changes in energies with increasing B field also are shown in Figs. 1 and 2 for the 1s2p−1 state of He and the 1s2p−1 3d−2 state of Li, respectively. Clearly, the TPSS (meta-GGA) and PBE (GGA) results are very close; the two curves are essentially indistinguishable. Moving from LDA to GGA significantly improves the accuracy of atomic energies at low B field, but shifting from GGA to MGGA (hence including the kinetic energy density as an additional ingredient) does not offer a comparable gain. Also see tables in the Supplemental Material [59]. Published results for correlated wave-function methods (WFMs) are only for 0 B 100 a.u. for He and Li, and for the Be atom, 0 B 10 a.u. From Figs. 1 and 2, it is hard to tell which methods of those we are comparing agree better with those data points (indicated by “+”). All curves essentially are superimposed upon the quoted data. It is easy to see why different methods do not give drastically different outcomes. Equations (1) and (2) show that the B field enters the system Hamiltonians in the same way for all the methods. The field only modifies the single-electron operators, with no effect on the two-electron operator. Effects of a B field on the electron-electron interaction can occur only indirectly through changing the many-electron wave function or, in DFT, the KS orbitals. Since all methods handle the oneelectron part exactly and it contributes a large portion of the total energy, results from different methods are not expected FPDQMC2 to differ too much as long as the two-electron energy is treated reasonably. Figure 3 illustrates these points by giving separated one-electron and two-electron energies in the He atom triplet state, 1s2p−1 . The curves for one-electron energies, labeled by E1e , are nearly indistinguishable by method. Differences arise mainly from the two-electron energies, E2e , precisely because of differing evaluations of the Hartree-plus-xc energy, EHxc . Detailed examination of the contributions to EHxc leads, however, to what seem to be previously undocumented results. Of course, both the Hartree term and the X term are much larger (in magnitude) than the C energy. In HF calculations, electron self-interaction energy ESI is canceled exactly by the selfexchange part of Ex and Ec is totally neglected. In DFT calculations, with either local or semilocal approximate xc functionals, the cancellation is incomplete, causing the so-called self-interaction error. (We discuss only integer self-interaction error, which is unexplored in the large B-field context, and leave noninteger self-interaction error (SIE) [62] to separate study.) Nevertheless, some gradient-dependent functionals do behave remarkably well in some circumstances. An example is the He atom low-field ground state 1s 2 . Define the selfinteraction energy in terms of the KS orbital self-repulsions as 1 2 i N ESI := |φi (r)φi (r )|2 dr dr . |r − r | (22) Exact self-interaction cancellation for a two-electron singlet is Ex [n]/J [n] = − 12 (or ESI /Ex = −1). The PW92 LDA X functional achieves 87% cancellation for this He state at B = 0, compared to more than 99% for B < 0.4 a.u. for both the PBE and TPSS X functionals. This approximate cancellation breaks down completely when B > 10 a.u. (see Fig. 4). Both 022504-6 COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . . PHYSICAL REVIEW A 90, 022504 (2014) TABLE III. Atomic energies of the helium atom in B fields by different methods. (Energies in hartree, B in a.u. The xc functionals are PW92 (LDA), PBE (GGA), and TPSS (mGGA), respectively. HGGA comprises exact exchange with the TPSS correlation functional.) State B (a.u.) HF, present 0 0.08 0.1 0.5 0.8 1 2 5 8 10 20 50 80 100 800 1000 −2.861679 −2.860417 −2.859709 −2.814450 −2.746839 −2.688884 −2.289144 −0.532442 1.591275 3.110634 11.319611 38.143906 66.092087 85.004179 770.54396 968.44540 0 0.08 0.1 0.5 0.8 1 2 5 8 10 20 50 80 100 800 1000 −2.131347 −2.236463 −2.259234 −2.615549 −2.830207 −2.959686 −3.502049 −4.617248 −5.400409 −5.829510 −7.427702 −10.264491 −12.101321 −13.076652 −26.126547 −28.032093 0 0.08 0.1 0.5 0.8 1 5 8 10 50 80 100 800 1000 −2.055211 −2.166315 −2.187305 −2.500874 −2.687529 −2.800387 −4.276634 −4.987052 −5.378085 −9.455332 −11.154700 −12.058706 −24.229300 −26.015249 HF, lit.a LDA GGA MGGA HGGA −2.83445 −2.83308 −2.83232 −2.78378 −2.71240 −2.65177 −2.23932 −0.45536 1.68598 3.21416 11.44959 38.29203 66.23035 85.13000 770.0003 967.7239 −2.89294 −2.89160 −2.89085 −2.84347 −2.77363 −2.71423 −2.30882 −0.54618 1.57640 3.09325 11.28092 38.01980 65.87939 84.73367 768.7686 966.3349 −2.90966 −2.90837 −2.90764 −2.86151 −2.79320 −2.73491 −2.33511 −0.58565 1.52752 3.03926 11.20842 37.91648 65.75751 84.60239 768.5391 966.0951 −2.90483 −2.90358 −2.90289 −2.85816 −2.79116 −2.73364 −2.33581 −0.58245 1.53956 3.05817 11.26517 38.08790 66.03580 84.94788 770.4910 968.3932 −2.08231 −2.18772 −2.21064 −2.56811 −2.78140 −2.90948 −3.44367 −4.54168 −5.31616 −5.74199 −7.33773 −10.20901 −12.09497 −13.10498 −27.19100 −29.33698 −2.13734 −2.24267 −2.26552 −2.62156 −2.83493 −2.96349 −3.50238 −4.61775 −5.40773 −5.84283 −7.47698 −10.42916 −12.37538 −13.41981 −28.12546 −30.38698 −2.13855 −2.24357 −2.26631 −2.62180 −2.83575 −2.96467 −3.50386 −4.61440 −5.39922 −5.83143 −7.45563 −10.39149 −12.32654 −13.36467 −27.95970 −30.20117 −2.13508 −2.24045 −2.26330 −2.62091 −2.83658 −2.96675 −3.51208 −4.63229 −5.41799 −5.84823 −7.44960 −10.28980 −12.12807 −13.10402 −26.15762 −28.06336 −2.00274 −2.11480 −2.13603 −2.45423 −2.64415 −2.75888 −4.24997 −4.96603 −5.36068 −9.51569 −11.27323 −12.21479 −25.35462 −27.35677 −2.05764 −2.17116 −2.19274 −2.51315 −2.70305 −2.81783 −4.32273 −5.05093 −5.45305 −9.70756 −11.51580 −12.48636 −26.15001 −28.25025 −2.06192 −2.17478 −2.19618 −2.51430 −2.70317 −2.81753 −4.31868 −5.04363 −5.44360 −9.66842 −11.46230 −12.42486 −25.95352 −28.02904 −2.05915 −2.17102 −2.19218 −2.50754 −2.69480 −2.80795 −4.28822 −5.00052 −5.39247 −9.47577 −11.17646 −12.08100 −24.25357 −26.03943 WFMb 1s 2 1s2p−1 1s3d−2 a b −2.861679996 −2.860417861 −2.859709376 −2.814450946 −2.746839677 −2.688884848 −2.289144423 −0.532445132 1.591274097 3.110633781 11.319608967 38.143903320 66.092085756 85.004177725 −2.2353 −2.8301 −5.4000 −12.1011 −26.1264 −2.1659 −2.6871 −4.9866 −11.1540 −24.2283 −2.903351 −2.902083 −2.855859 −2.787556 −2.729508 −2.329780 −0.574877 3.064582 11.267051 38.076320 84.918313 −2.133149 −2.238504 −2.620021 −2.835619 −2.965504 −3.508911 −4.625491 −5.839475 −7.440556 −10.28410 −13.10478 −2.055635 −2.166519 −2.502362 −2.689916 −2.803296 −4.284050 −5.387931 −9.476057 −12.088566 For the 1s 2 state, data are from Ref. [17]; for other states, data are from Ref. [20]. Data are from Ref. [21]. the local and semilocal exchange functionals overestimate |Ex | when B is very large, thereby overcompensating for the self-interaction error. (Even for B = 0, spurious compensation behavior from simple xc functionals has been known for a long time. LDA xc gets the H atom total energy roughly right by having a completely spurious nonzero Ec [63].) Figure 5 displays the ratio |Ex /ESI | for the fully spinpolarized state of the He atom, 1s2p−1 . Since the antisymmetric wave function makes the two electrons well separated 022504-7 WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) TABLE IV. Atomic energies of the lithium atom in B fields by different methods. (Energies in hartree, B in a.u. The chosen xc functionals for LDA, GGA, and MGGA calculations are PW92, PBE, and TPSS functionals, respectively. For HGGA, we use exact exchange and TPSS correlation functional.) HF, lit.a LDA GGA MGGA HGGA −7.43275 −7.46856 −7.47740 −7.40878 −7.19621 −6.08810 −5.90114 −3.35784 3.49120 71.80766 939.55235 −7.43275 −7.46857 −7.47741 −7.40879 −7.19621 −6.08811 −5.90113 −3.35777 3.49120 71.807 939.54 −7.34328 −7.37939 −7.39479 −7.33924 −7.14208 −6.04813 −5.86217 −3.32762 3.50491 71.67426 938.1123 −7.46217 −7.49852 −7.51699 −7.46832 −7.28261 −6.21946 −6.03738 −3.54382 3.21595 71.03954 936.0146 −7.48906 −7.52548 −7.54402 −7.49671 −7.31483 −6.26245 −6.08170 −3.60153 3.13744 70.88976 935.8673 −7.48223 −7.51840 −7.52833 −7.45885 −7.24582 −6.13916 −5.95240 −3.41095 3.43600 71.74951 939.4964 −7.477766 −7.517154 −7.528055 −7.458550 −7.244919 −6.136918 −5.949297 −3.406556 −7.4763360 −7.4777957 −7.5122102 −7.5137817 −7.5216127 −7.5235946 −7.4529046 −7.2397460 −7.36507 −7.44174 −7.58787 −7.66652 −7.66245 −6.94229 −6.79515 −4.61775 1.70566 68.17349 930.84309 −7.36509 −7.44176 −7.58790 −7.66653 −7.66246 −6.94230 −6.79517 −4.61777 1.70565 68.1735 930.84308 −7.27909 −7.35647 −7.50718 −7.58789 −7.58026 −6.83784 −6.68788 −4.48455 1.86892 68.30834 929.6385 −7.39752 −7.47501 −7.62468 −7.70476 −7.69981 −6.97522 −6.82790 −4.65417 1.64374 67.80872 927.9366 −7.42270 −7.49961 −7.64670 −7.72652 −7.72379 −7.00375 −6.85670 −4.68388 1.61495 67.78076 927.9620 −7.41506 −7.49225 −7.64052 −7.72182 −7.72297 −7.01341 −6.86726 −4.69801 1.61654 68.06848 930.7310 −7.407126 −7.484773 −7.634547 −7.716679 −7.715709 −7.002346 −6.855410 −4.684076 −7.4088037 −7.4097907 −7.4858382 −7.4869343 −7.6341245 −7.6362483 −7.7151944 −7.7137805 State B (a.u.) HF, present WFM 1b WFM 2c WFM 3d 1s 2 2s 0 0.1 0.5 1 2 5 5.4 10 20 100 1000 1s 2 2p−1 0 0.1 0.5 1 2 5 5.4 10 20 100 1000 1s2p−1 3d−2 0 0.1 1 2 5 10 100 1000 −5.9448544 −5.8555577 −3.4020661 3.4446412 71.7573135 −6.8503299 −6.8361629 -4.6811073 1.6398139 −5.08377 −5.08379 −5.00614 −5.09830 −5.09921 −5.09177 −5.142319 −5.32138 −5.32140 −5.24405 −5.33574 −5.33547 −5.33166 −5.341030 −6.57079 −6.57081 −6.48534 −6.58053 −6.57729 −6.58846 −6.582361 −7.52002 −7.52003 −7.42928 −7.52682 −7.52159 −7.54072 −7.530125 −9.57693 −9.57694 −9.47135 −9.58275 −9.56775 −9.60391 −9.591769 −11.93900 −11.93902 −11.82087 −11.95879 −11.92698 −11.97189 −11.957294 −27.01926 −27.0192 −27.07687 −27.47894 −27.36026 −27.06791 −60.05888 −60.0589 −62.02394 −63.34180 −63.0414 −60.1164 a Data are from Ref. [16]. Data are from Ref. [24]. c Data are from Ref. [25]. d Data are from Ref. [23]. b in space, the true two-electron exchange contributes only a small fraction of Ex . About 90% or more of Ex is simply to counter ESI . Here again, both the PBE and TPSS exchange functionals do reasonably well when B is small, but fail badly for larger B. HF exchange exhibits nonmonotonic behavior in Ex /ESI , with near saturation at B = 10 a.u. All three approximate xc functionals checked here lack this characteristic. Instead, they grossly overshoot the X energy in magnitude at large B. Hence, ESI cannot be expected to be canceled numerically by a local or semilocal X functional when B is significant (even if that cancellation were to be from some mixture of X and C as in the B = 0 case). Because of the field-induced transverse compression of the density, the self-interaction problem appears in particularly severe form in nB-DFT calculations. Since ExXX , Eq. (10) and EcTPSS both are free from electron self-interaction by design, the HGGA combination of the two, Eq. (9), also is self-interaction free. Limitations of this combination were addressed in Sec. II, but they do not pose a problem to our present study on magnetized atoms and ions. HGGA results are presented in Tables I, III, and IV, and in the Supplemental Material [59]. For Z 4, our HGGA results agree better with correlated wave-function results [19–27] than those from other DFT approximations we investigated. For Z 5, we did not find correlated wave-function results with which to compare. Table I, however, shows that our HGGA results are the only ones from DFT which come quite close to the best available QMC calculations for the atomic energies of the B and C atoms in large B fields. The comparative computational cost advantage of HGGA is obvious. Detailed insight into the excellent behavior of this HGGA functional comes from comparison with the QMC study in Ref. [28]. Those authors observed that density-functional approximations cannot guarantee in all cases to “produce an 022504-8 COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . . PHYSICAL REVIEW A 90, 022504 (2014) FIG. 1. (Color online) Atomic total energies for the He atom fully spin-polarized state 1s2p−1 in magnetic fields from Hartree-Fock and DFT calculations with different functionals. FIG. 3. (Color online) One-electron and two-electron energies for the He atom fully spin-polarized state 1s2p−1 in magnetic fields from Hartree-Fock and different DFT exchange-correlation functionals. upper bound on the ground-state energy in magnetic fields . . . .” Such failures are examples of the well-known issue of possible loss of N representability in an approximate xc functional. However, comparison of the HGGA results with the released phase diffusion QMC (RPDQMC) results of Ref. [28] show that Etot,HGGA (B) never is below Etot,RPDQMC (B) (see Table I). The salient point is that while the key features of HGGA do not guarantee N representability, they certainly make it plausible. By construction Ex in HGGA is variationally exact. While EcTPSS is not guaranteed to be variationally exact, it is self-interaction-free, hence does not rely upon spurious overly strong (too negative) C to address SIE. Furthermore, Ec is a small fraction of ExXX , so one would not expect EcTPSS to drive the total energy below the exact energy. The numerical evidence is wholly consistent with this, so operationally there is evidence that HGGA is N representable. All of the foregoing DFT calculations rest on the nB-DFT approximation, with current effects in Exc neglected. The excellent agreement of our HGGA results with those from the most sophisticated correlated WFMs and QMC calculations reveals two messages. One is that the current dependence of Exc apparently can be neglected safely over a wide range of experimentally accessible B fields, at least for atomic systems. The nB-DFT approximation is a useful approach in which the major B-field effects are included. Another message is that the development of ground-state DFT functionals with B = 0, especially the nonlocal HGGA functionals which treat the self-interaction explicitly has in fact made accurate B > 0 DFT calculations available, at least for light atoms. This suggests that nonlocal DFT functionals may be a better starting point for devising a full-fledged CDFT functional which would also encompass the current effect in Exc . FIG. 2. (Color online) Atomic total energies for the Li atom fully spin-polarized state 1s2p−1 3d−2 in magnetic fields from HartreeFock and DFT calculations with different functionals. FIG. 4. (Color online) The ratio of exchange energy Ex over electron self-interaction energy ESI for the He atom singlet state 1s 2 in magnetic fields from Hartree-Fock and different DFT exchange functionals. The ratio is exactly −1 for Hartree-Fock calculation. 022504-9 WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) TABLE V. CDFT corrections to the LDA PW2 results within VRG approximation (parameters ncutoff = 0.001a0−3 , αcutoff = 2a0−1 VRG are used for the cutoff function, Exc in hartree). Atom FIG. 5. (Color online) The ratio of exchange energy Ex over electron self-interaction energy ESI for the He atom fully spinpolarized state 1s2p−1 in magnetic fields from Hartree-Fock and different DFT exchange functionals. C. Results for the CDFT VRG approximation For CDFT, there does not appear to have been any fully self-consistent calculations on atoms in large B fields based on the VRG functional [Eq. (11)]. There is related work on the perturbative implementation of the VRG functional by Orestes, Marcasso, and Capelle [33] for atomic ionization energies in vanishing B field. Gross and collaborators compared spin-DFT and CSDFT in several different systems, including light atoms (for the spurious energy splittings between current-carrying and zero-current ground states) [47], magnetic and nonmagnetic solids for the spin-orbit splittings, and orbital magnetic moments [64], both at B = 0 and the X-only functional ExXX implemented via the OEP procedure. They found little differences in the outcome from including jp dependence or not. For a two-dimensional quantum dot (2D-QD) in B, they added self-interaction-corrected LSDA correlation energies post hoc to the bare ExXX results to compare with QMC energies. Again they found little effect from the current density [65]. Similar results with the VRG functional recently were presented by Tellgren et al. [34] for molecular magnetizabilities, hypermagnetizabilities, and NMR shielding constants for B up to 1 a.u. An important formal point is that the VRG functional VRG was constructed as an additive term Exc to provide the current contribution to Exc at the LDA level of refinement. That is, VRG is supposed to be used in conjunction with an LDA functional. The VRG functional also is negative definite by design. Since LDA energies can be either above or below accurate results in both low- and high-field regimes (see tables), the VRG correction cannot always be a proper correction to bring LDA values closer to exact total energies. Table V lists representative results for the PW92 + VRG approximation for the fully spin-polarized states of the He and Li atoms. An estimate of the effects of jp contributions is to evaluate the VRG functional using the LDA Kohn-Sham orbitals. Results from that post-SCF procedure are listed in VRG . They are very close to SCF Table V as non-SCF Exc State He 1s2p−1 Li 1s2p−1 3d−2 B (a.u.) VRG Non-SCF Exc VRG SCF Exc 0 0.5 1 5 10 100 −0.0022 −0.0045 −0.0077 −0.036 −0.074 −0.81 −0.0021 −0.0047 −0.0081 0 2 5 10 −0.0070 −0.027 −0.065 −0.129 −0.0071 −0.029 VRG Exc whenever self-consistency can be reached. Evidently full self-consistency in CDFT calculations has little effect on the current correction term and, hence, the total energy. Since SCF convergence problems occur when B is sufficiently large because of the divergent vorticity ν, we have used non-SCF VRG Exc values in those cases for which self-consistency could not be reached. Similar observations were reported by Helbig et al. in 2D-QD, except that they used LSDA for the correlation energy [65]. This behavior is consistent with the variational nature of DFT. A small change in the density near the self-consistent density only results in a second-order energy change. Presumably the change in density from addition of current effects is small. Table V shows that VRG corrections depend strongly upon particular atomic configurations, but the overall values are small in magnitude compared to the variation among HF, LDA, and gradient-dependent DFT results when the B field is not too large (below roughly 1 a.u.). Within each configuration, the correction increases with increasing B field. For B < 5 a.u., the VRG correction is rather small compared with the discrepancy between the LDA and WFM results and does not bring them significantly closer. At B = 10 a.u., the VRG functional fortuitously gives a good correction to the Li atom 1s2p−1 3d−2 state, but at B = 100 a.u., LDA itself gives an accurate energy for the He atom 1s2p−1 state. At that point, the VRG correction can only degrade the results. Overall, there is no obviously systematic behavior. As a by-product of numerical issues, we were forced to study the vorticity ν, which turns out to be a rather difficult variable to handle. For example, in Fig. 6 we display various quantities in CDFT along z and radially along ρ for the He atom 1s2p−1 state in B = 1 a.u. Because jp (0,z) = 0, it is not displayed. However, ν is not zero along z. Instead, it diverges at the two poles of the atom and so does the VRG xc energy density gLCH (n) |ν|2 . Notice that n(r) decays rapidly along z. At z = 3 bohrs, it is smaller than 10−4 bohr−3 , but the VRG xc energy density gLCH (n) |ν|2 becomes increasingly large. It seems rather peculiar that the largest corrections from the VRG xc functional is at those places where both the electron density and the current density are nearly zero. Further analysis shows that this strange behavior is rooted in the choice of ν as a basic variable in the VRG functional. 022504-10 COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . . PHYSICAL REVIEW A 90, 022504 (2014) we view the CDFT corrections for total energies as subject to caution; they are semiquantitative at best. V. DISCUSSION AND CONCLUSIONS FIG. 6. (Color online) Various quantities (electron density n, paramagnetic current density jp , vorticity ν, and the magnitude of the current correction to the exchange-correlation energy density, |gν 2 |, in the VRG functional) for the He atom 1s2p−1 state in B = 1 a.u. All quantities are evaluated from the LDA KS orbitals and plotted along the z and ρ axes (cylindrical coordinates). One may argue that at the places of such a low density the prefactor g(n) is not well defined in practice because gLCH (n) (or other fitted functions) only used data points with rs 10 bohrs [49]. This is indeed a problem, but unfortunately we do not know the form of g(n) in the low-density limit. The picture of HEG at a very low electron density in a strong B field may be inappropriate anyway. To address this divergence problem numerically, we introduced a rapidly decaying cutoff function gcutoff . Details were given in the methodology section. However, there is no obvious, physically grounded set of criteria for choosing the cutoff parameters ncutoff and αcutoff . We tested different parameter sets by using the He atom 1s2p−1 state in B = 1 a.u. Different cutoff parameters generate quite different results see Table VI). This outcome is distinct from the findings in Ref. [34] presumably for two reasons. Our study focuses on total energies, theirs on response (e.g., magnetizability). Second, we have gone to much higher fields. The test results shown in Table VI are at B = 1 a.u., the highest field used in Ref. [34] but at the low end of the range we have examined. Comparison with the He 1s2p−1 LDA energy at B = 1 a.u. (−2.909 48 hartree; recall Table III) shows the undesirable effect upon total energies of the dependence of the VRG correction upon regularization. Because of that dependence, TABLE VI. Effect of cutoff parameters on CDFT VRG corrections for the He atom 1s2p−1 state in magnetic field B = 1 a.u. (energies in hartree). ncutoff (units of a0−3 ) 0.005 0.001 0.001 0.0001 0.00001 αcutoff (units of a0−1 ) VRG Non-SCF Exc 2.0 2.0 1.0 2.0 2.0 −0.004 −0.008 −0.010 −0.025 −0.064 It is perhaps unsurprising that the VRG functional fails when applied to atomic systems in large B. The functional was developed from the study of the moderately dense to dense HEG in low B, for which Landau orbitals were used as approximations. This physical situation is quite different from a finite system such as an atom. First, n(r) and jp (r) vary considerably within an atom, and the low-density regions (rs > 10 a0 ) are non-negligible. Secondly, there is not a direct relationship between jp (r) and the external B field as there is for the HEG. The question whether the electron gas remains homogeneous after imposing a substantial B field is even unclear. If the field were to induce some form of crystallization, the basic picture on which the VRG functional is based is lost. As Pittalis et al. pointed out, the B field can induce derivative discontinuities in the xc energy density. Such discontinuities make practical calculations extremely difficult [47]. Our analysis and numerical studies also suggest the picture of Landau orbitals used for the HEG may not be applicable at all for the atomiclike systems. Unlike the LDA, also based on the HEG, the VRG functional may be too simple to encompass the essential physics of the atomic systems. Alternatives are admittedly difficult to imagine. A more fundamental question is whether ν is a suitable basic variable in gauge-invariant CDFT as Vignale and Rasolt required [6]. While appealing on purely theoretical grounds, our numerical results on atomic systems in B suggest ν is, at best, an awkward choice. It is ν as the basic variable in the VRG functional which drives unphysical results and the need for cutoffs. While ν appears to be a misbehaved quantity for atomic systems, jp is still a well-behaved variable, but it suffers from gauge dependence. Recently, Tellgren et al. reexamined the possibility of using (n,jp ) as basic variables, and found that the gauge problem can be circumvented by using modified conjugated potentials [66]. Perhaps one should take a step back from ν to jp as a basic variable to conceive viable CDFT functionals. We note, however, that this is not the route suggested in Ref. [67]. Pan and Sahni also formulated CDFT by using the electron density n and physical current density j [68], providing another interesting idea, though with some restrictions [69], the consequences of which remain to be explored (see, for example, Ref. [70]). There is a larger issue. To a considerable extent, the present study calls into question the utility of CDFT (in any form) itself. Exact exchange combined with the TPSS correlation functional yields a self-interaction-free HGGA functional, Eq. (9). Energies from this functional used as an nB-DFT approximation agree quite well with results from correlated WFMs and QMC. That agreement indicates that the current correction to Exc from CDFT is rather small over a wide range of B. As a result, the nonlocal HGGA functional may be both a practical alternative to CDFT (at least over the range of fields which we have considered) as well as a better starting point than LDA to construct full-fledged CDFT functionals. One final remark is pertinent to the results presented here. Neither relativistic effects nor finite nuclear mass effects 022504-11 WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) have been considered. Those effects can be important for matter in superstrong fields (B > 104 a.u.), in which regime the adiabatic approximation will be applicable. But for the field strengths covered in this paper, both effects should be negligible. ACKNOWLEDGMENTS Helpful discussions with Gao Xianlong (Zhejiang Normal University), Wenjian Liu (Peking University), X.-Y. Pan (Ningbo University), and Jian Wang (Huzhou University) are acknowledged with thanks. This work was funded by Zhejiang Provincial Natural Science Foundation of China under Grant No. LY13A050002 (W.Z. and L.Z.), in part by the National Natural Science Foundation of China Grant No. 11274085 (W.Z.), and in part by U.S. Department of Energy Grant No. DE-SC-0002139 (S.B.T.). In early portions of this work both W.Z. and S.B.T. were supported by U.S. National Science Foundation Grant No. DMR-0218957. APPENDIX The KL basis sets start with a base sequence of orbital exponents defined as those nearest to spherical, i.e., closest to αj = βj ∀j . The derived sequences compress the transverse (radial) functions (the αj ) relative to the corresponding longitudinal (axial) ones (βj ). Thus, the KL basis sets have several (one to five) exponent sequences of the form μ2,3 = ±0.2, μ4,5 = ±0.4, (A1) with αj,1 ≈ βj ∀j and KL a parameter to be optimized at one fixed B. Note that in each sequence the differences between the transverse (αj ) and longitudinal exponents (βj ) are the same for all members of that sequence; neither μ nor KL is indexed by j . To keep the basis size within reason, the second and third sequences ( = 2,3) are half the length (with doubled spacing) of the first sequence, while the fourth and fifth sequences ( = 4,5) are one-quarter as long (with quadrupled spacing) as the base sequence. We reworked the KL scheme as follows. First was some numerical experimentation on the H atom. At B = 0 a.u., an even-tempered Gaussian sequence of length Nb = 16 and rule of formation [71], αj = βj ≡ pq j , j = 1,2, . . . ,Nb , ln p = a ln(q − 1) + a , ln(ln q) = b ln Nb + b , a = 0.3243, b = −0.4250, := βj + (1 + μ1 )(B), (A2) a = −3.6920, b = 0.9280, gives Etot = −0.499 999 92 hartree. However, at B = 10 a.u., it gives a 24% error compared to the known value Etot (B = 10) = −1.747 797 163 714 hartree [60]. The task therefore is to discover as near-optimum asymmetry between the transverse and longitudinal manifolds as feasible. Within only the basic KL sequence and the Nb = 16 spherical set, we searched in the parameter space αj to minimize the total energy of the hydrogen atom at B = 10 a.u. (A3) where b(γ ) = −0.16[tan−1 (γ )]2 + 0.77 tan−1 (γ ) + 0.74. (A4) is the reduced field strength for an In general γ = effective nuclear charge Zeff . For the innermost electrons, Zeff is close to the bare nuclear charge. For other electrons, it is close to the nuclear charge minus the number of inner shell electrons (screening effects). In practice, however, we found that nominal Zeff values (bare nuclear charge), such as just given, are good enough. The expressions for αj , Eq. (A3), and b(γ ), Eq. (A4), came from fits to the optimized energies of the H isoelectronic series (H, He+ , Li++ , Be+++ , C5+ , and O7+ in reduced fields γ = 0.1, 0.5, 1, 2, 5, 10, 20, 50, 100, 200, 500, 800, 1000, 2000, and 4000. Those optimizations also started from the Nb = 16 basis set. Equation (A3) with μ1 = 0 defines the first sequence. Other sequences are defined analogously with the KL basis set, μ = ±0.2, ± 0.4 for the second, third, fourth, and fifth sequences (with increased spacing as mentioned in the main text), respectively. Similar to the KL basis sets, we used ETG sequences for the longitudinal exponents βj ; recall Eq. (A2). Since the magnetic field does not change the confinement along the z direction, we expect it to have little effect upon ETG exponents. In cases for which the electron density has diffuse, nonzero orbital angular momentum contributions, extrapolation of the tempering to include a small number of negative j values is used. Testing confirmed that the resulting basis sets work as well for multielectron atoms as for one-electron systems. Our criterion for basis set error is that it be less than 1 μHartree in the whole B range investigated (0 B 2000 a.u.) for the hydrogen atom, with respect to more accurate algebraic results [60]. We find that the HF energies for multielectron atoms from our basis sets are nearly indistinguishable from full numerical two-dimensional mesh results [16]. 2 B/Zeff αj, = βj + (1 + μ )KL , μ1 = 0.0, and found that the optimized difference αj − βj is not constant at fixed B as the KL sequences suggest. Particularly for the smaller αj values, the fractional difference (αj − βj )/αj is quite large, as much as ≈ 0.95. The behavior is quite understandable. A small-exponent basis function samples a large volume far from the nucleus, where the B field overpowers the nuclear attraction such that the distortion from the field-free spherical symmetry will be relatively larger than for the region sampled by larger-exponent basis functions. In the limit βj → 0, which can be thought of either as the large B limit or zero nuclear charge, the electron wave function is a Landau orbital, with an exponential parameter αj ≡ aB = B/4 (aB is the scale parameter in the Landau orbital). The opposite limit, βj → ∞, corresponds to B = 0 a.u., for which spherical symmetry is restored, αj = βj . Thus, a convenient choice for orbital exponent asphericity is a scaling of aB . By consideration of such large and small B-field limits and use of a nonlinear fit to calculated one-electron system results, we reached a prescription for nearly optimized exponents, namely, B 4 βj −1/2 4 βj −2 αj = βj + + 1+ 4 1+ 20 b(γ ) B b(γ ) B 022504-12 COMPARATIVE STUDIES OF DENSITY-FUNCTIONAL . . . PHYSICAL REVIEW A 90, 022504 (2014) [1] P. Hohenberg and W. Kohn, Phys. Rev. 136, B864 (1964); W. Kohn and L. J. Sham, ibid. 140, A1133 (1965). [2] P. B. Jones, Phys. Rev. Lett. 55, 1338 (1985). [3] B. M. Relovsky and H. Ruder, Phys. Rev. A 53, 4068 (1996). [4] M. Braun, Phys. Rev. A 65, 033415 (2002). [5] Z. Medin and D. Lai, Phys. Rev. A 74, 062507 (2006). [6] G. Vignale and M. Rasolt, Phys. Rev. Lett. 59, 2360 (1987); ,Phys. Rev. B 37, 10685 (1988); G. Vignale, M. Rasolt, and D. J. W. Geldart, Adv. Quantum Chem. 21, 235 (1990). [7] Wuming Zhu and S. B. Trickey, J. Chem. Phys. 125, 094317 (2006). [8] M. Taut, P. Machon, and H. Eschrig, Phys. Rev. A 80, 022517 (2009). [9] E. H. Lieb and R. Schrader, Phys. Rev. A 88, 032516 (2013). [10] P. Pr¨oschel, W. R¨osner, G. Wunner, H. Ruder, and H. Herold, J. Phys. B 15, 1959 (1982). [11] D. Neuhauser, S. E. Koonin, and K. Langanke, Phys. Rev. A 36, 4163 (1987). [12] D. Lai, E. E. Salpeter, and S. L. Shapiro, Phys. Rev. A 45, 4832 (1992). [13] H. Ruder, G. Wunner, H. Herold, and F. Geyer, Atoms in Strong Magnetic Fields (Springer-Verlag, Berlin, 1994). [14] J. S. Heyl and L. Hernquist, Phys. Rev. A 58, 3567 (1998); A. Thirumalai and J. S. Heyl, ibid. 79, 012514 (2009). [15] M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys. Rev. A 54, 219 (1996). [16] M. V. Ivanov and P. Schmelcher, Phys. Rev. A 57, 3793 (1998); ,60, 3558 (1999); ,61, 022505 (2000); ,Eur. Phys. J. D 14, 279 (2001); ,J. Phys. B 34, 2031 (2001). [17] L. B. Zhao and P. C. Stancil, Phys. Rev. A 77, 035401 (2008). [18] D. Engel and G. Wunner, Phys. Rev. A 78, 032515 (2008). [19] A. Scrinzi, Phys. Rev. A 58, 3879 (1998). [20] M. D. Jones, G. Ortiz, and D. M. Ceperley, Phys. Rev. A 59, 2875 (1999). [21] W. Becken, P. Schmelcher, and F. K. Diakonos, J. Phys. B 32, 1557 (1999); W. Becken and P. Schmelcher, ibid. 33, 545 (2000); ,Phys. Rev. A 63, 053412 (2001). [22] X. Wang, J. Zhao, and H. Qiao, Phys. Rev. A 80, 053425 (2009). [23] X. Guan and B. Li, Phys. Rev. A 63, 043413 (2001). [24] O.-A. Al-Hujaj and P. Schmelcher, Phys. Rev. A 70, 033411 (2004). [25] X. Wang and H. Qiao, Phys. Rev. A 75, 033421 (2007). [26] X. Guan, B. Li, and K. T. Taylor, J. Phys. B 36, 2465 (2003). [27] O.-A. Al-Hujaj and P. Schmelcher, Phys. Rev. A 70, 023411 (2004). [28] S. B¨ucheler, D. Engel, J. Main, and G. Wunner, Phys. Rev. A 76, 032501 (2007). [29] C. Schimeczek, S. Boblest, D. Meyer, and G. Wunner, Phys. Rev. A 88, 012509 (2013). [30] S. M. Colwell and N. C. Handy, Chem. Phys. Lett. 217, 271 (1994); A. M. Lee, S. M. Colwell, and N. C. Handy, ibid. 229, 225 (1994). [31] A. M. Lee, N. C. Handy, and S. M. Colwell, J. Chem. Phys. 103, 10095 (1995). [32] S. M. Colwell, N. C. Handy, and A. M. Lee, Phys. Rev. A 53, 1316 (1996); A. G. Ioannou, S. M. Colwell, and R. D. Amos, Chem. Phys. Lett. 278, 278 (1997). [33] E. Orestes, T. Marcasso, and K. Capelle, Phys. Rev. A 68, 022105 (2003). [34] E. I. Tellgren, A. M. Teale, J. W. Furness, K. K. Lange, U. Ekstr¨om, and T. Helgaker, J. Chem. Phys. 140, 034101 (2014). [35] B. G. Johnson, P. M. W. Gill, and J. A. Pople, J. Chem. Phys. 98, 5612 (1993). [36] I. H. Lee and R. M. Martin, Phys. Rev. B 56, 7197 (1997). [37] S. Kurth, J. P. Perdew, and P. Blaha, Int. J. Quantum Chem. 75, 889 (1999). [38] J. P. Perdew and K. Schmidt, Density Functional Theory and Its Application to Materials, edited by V. Van Doren, C. Van Alsenoy, and P. Geerlings, AIP. Conf. Proc. No. 577 (AIP, Melville, NY, 2001), p. 1. [39] J. P. Perdew and Y. Wang, Phys. Rev. B 45, 13244 (1992). [40] J. P. Perdew, K. Burke, and M. Ernzerhof, Phys. Rev. Lett. 77, 3865 (1996); ,78, 1396 (1997). [41] J. Tao, J. P. Perdew, V. N. Staroverov, and G. E. Scuseria, Phys. Rev. Lett. 91, 146401 (2003). [42] O. A. Vydrov, G. E. Scuseria, J. P. Perdew, A. Ruzsinszky, and G. I. Csonka, J. Chem. Phys. 124, 094108 (2006). [43] A. G¨orling, Phys. Rev. Lett. 83, 5459 (1999). [44] S. K¨ummel and L. Kronik, Rev. Mod. Phys. 80, 3 (2008). [45] L. A. Constantin, J. M. Pitarke, J. F. Dobson, A. Garcia-Lekue, and J. P. Perdew, Phys. Rev. Lett. 100, 036401 (2008). [46] J. D. Talman and W. F. Shadwick, Phys. Rev. A 14, 36 (1976). [47] S. Pittalis, S. Kurth, N. Helbig, and E. K. U. Gross, Phys. Rev. A 74, 062511 (2006). [48] J. B. Krieger, Y. Li, and G. J. Iafrate, Phys. Rev. A 45, 101 (1992). [49] G. Vignale, M. Rasolt, and D. J. W. Geldart, Phys. Rev. B 37, 2502 (1988). [50] J. Tao and J. P. Perdew, Phys. Rev. Lett. 95, 196403 (2005). [51] J. Tao and G. Vignale, Phys. Rev. B 74, 193108 (2006). [52] K. Higuchi and M. Higuchi, Phys. Rev. B 74, 195122 (2006). [53] M. D. Jones, G. Ortiz, and D. M. Ceperley, Int. J. Quantum Chem. 64, 523 (1997). [54] D. Porezag and M. R. Pederson, Phys. Rev. A 60, 2840 (1999). [55] C. Aldrich and R. L. Greene, Phys. Status Solidi B 93, 343 (1979). [56] P. Schmelcher and L. S. Cederbaum, Phys. Rev. A 37, 672 (1988). [57] U. Kappes and P. Schmelcher, J. Chem. Phys. 100, 2878 (1994). [58] Yu. P. Kravchenko and M. A. Liberman, Int. J. Quantum Chem. 64, 513 (1997). [59] See Supplemental Material at http://link.aps.org/supplemental/ 10.1103/PhysRevA.90.022504 for extensive tabulations of HF and DFT calculated data for the He, Li, Be, B, and C atoms and the positive ions, Li+ , Be+ , and B+ . [60] Y. P. Kravchenko, M. A. Liberman, and B. Johansson, Phys. Rev. A 54, 287 (1996). [61] R. W. Danz and M. L. Glasser, Phys. Rev. B 4, 94 (1971). [62] Weitao Yang, P. Mori-Sanchez, and A. J. Cohen, J. Chem. Phys. 139, 104114 (2013). [63] S. B. Trickey, Phys. Rev. Lett. 56, 881 (1986). [64] S. Sharma, S. Pittalis, S. Kurth, S. Shallcross, J. K. Dewhurst, and E. K. U. Gross, Phys. Rev. B 76, 100401(R) (2007). [65] N. Helbig, S. Kurth, S. Pittalis, E. R¨as¨anen, and E. K. U. Gross, Phys. Rev. B 77, 245106 (2008). [66] E. I. Tellgren, S. Kvaal, E. Sagvolden, U. Ekstr¨om, A. M. Teale, and T. Helgaker, Phys. Rev. A 86, 062506 (2012). 022504-13 WUMING ZHU, LIANG ZHANG, AND S. B. TRICKEY PHYSICAL REVIEW A 90, 022504 (2014) [67] K. Higuchi, M. Miyasita, and M. Higuchi, Int. J. Quantum Chem. 110, 2286 (2010). [68] X.-Y. Pan and V. Sahni, Int. J. Quantum Chem. 110, 2833 (2010); ,J. Phys. Chem. Sol. 73, 630 (2012); ,Phys. Rev. A 86, 042502 (2012); V. Sahni and X.-Y. Pan, ibid. 85, 052502 (2012). [69] X.-Y. Pan and V. Sahni, Int. J. Quantum Chem. 114, 233 (2014). [70] A. Laestadius and M. Benedicks, Int. J. Quantum Chem. 114, 782 (2014). [71] M. W. Schmidt and K. Ruedenberg, J. Chem. Phys. 71, 3951 (1979). 022504-14
© Copyright 2024 ExpyDoc