Implementing and testing the AM05 spin density functional

(1)

Implementing and testing the AM05 spin

density functional

Ann E. Mattsson and Rickard Armiento

Linköping University Post Print

N.B.: When citing this work, cite the original article.

Original Publication:

Ann E. Mattsson and Rickard Armiento, Implementing and testing the AM05 spin density

functional, 2009, Physical Review B. Condensed Matter and Materials Physics, (79), 15,

155101.

http://dx.doi.org/10.1103/PhysRevB.79.155101

Copyright: American Physical Society

http://www.aps.org/

Postprint available at: Linköping University Electronic Press

(2)

Implementing and testing the AM05 spin density functional

Ann E. Mattsson1,

_*

_{and Rickard Armiento}2

1_{Multiscale Dynamic Materials Modeling MS 1322, Sandia National Laboratories, Albuquerque, New Mexico 87185-1322, USA} 2_{Physics Institute, University of Bayreuth, D-95440 Bayreuth, Germany}

共Received 4 December 2008; revised manuscript received 18 February 2009; published 2 April 2009兲 We show that the spin density generalization of the AM05 density functional 关R. Armiento and A. E. Mattsson, Phys. Rev. B 72, 085108共2005兲兴 predicts the correct ground spin state for iron, a system known to be heavily dependent on proper spin treatment. Using the fundamental assumptions in the subsystem functional scheme, we resolve an ambiguity in how to treat the separate spin densities in AM05 but also show that the other less preferred treatments give no significantly different numerical outcome of the iron body-centered-cubic and face-centered-body-centered-cubic test cases. Details and formulas are given to aid in the implementation of functionals in general, and the spin-resolved AM05 exchange-correlation potentials in particular, into different types of computer codes.

DOI:10.1103/PhysRevB.79.155101 PACS number共s兲: 71.15.Mb, 31.15.ej, 71.45.Gm, 75.50.Bb

I. INTRODUCTION

Due to its ability to treat a wide range of systems fairly accurately at relatively low computational cost, density-functional theory 共DFT兲共Refs. 1 and 2兲 and its spin

formulation3,4_{has become the foundation of most large scale}

quantum-mechanical simulations. The key to accuracy in DFT based simulations is the approximation made for the exchange-correlation density functional.

A number of recent studies have confirmed the excellent performance of the AM05 functional5 _{for nonmagnetic}

solids.6 _{Furthermore, Ropo et al.}7 _{recently investigated a}

limited set of properties also of magnetic solids using a spin-resolved version of AM05 共the xscst version discussed in Sec.II兲. Recent work by Haas et al.8_{confirms the findings in}

Refs. 6and7. The present paper presents the necessary de-tails in the derivation of spin-resolved AM05 and discusses its implementation into various types of computer codes. In addition, we test the spin-resolved functional for proper pre-diction of body-centered-cubic共bcc兲 and face-centered-cubic 共fcc兲 phases of iron. This is a classical test case, for which the local spin density approximation共LSDA兲 fails to give the correct ground spin state while, e.g., the Perdew, Burke, and Ernzerhof共PBE兲 functional9_{reproduces it correctly.}

In Sec.IIwe give an overview of the AM05 density func-tional from the subsystem funcfunc-tional viewpoint and we de-rive several forms of AM05 that naturally extends its use to systems where a spin-resolved electron density is needed for proper treatment. The best form to use is identified based on fundamental subsystem functional principles. SectionIII dis-cusses the general implementation of spin-resolved function-als in different types of codes currently in use, and we give specific formulas for the spin-resolved AM05. Section IV

gives the results of benchmark tests on iron with the positive result that spin-resolved AM05 correctly identifies the ground spin state. We summarize and discuss our findings in Sec. V.

II. AM05 AND SPIN

The AM05 functional5 _consists _of _a _subsystem

functional10,11_{for interior regions, a subsystem functional for}

edge regions, and an interpolation index that in each point in a system determines the ratio of edge and interior functionals to use. A subsystem functional is designed to properly treat a particular type of physical situation. AM05 thus combines a treatment based on the uniform electron gas, describing a situation where electrons are free to move in all directions, appropriate for interior regions, with a treatment of a very different edgelike situation, where electrons are confined by a surface, restricting their motion in one direction. AM05 is using only the density, n共r兲, and the magnitude of its gradi-ent, 兩ⵜn共r兲兩, to determine the contribution to the total exchange-correlation energy in each point, r, of the system, just as common generalized gradient approximation 共GGA兲 functionals do.

As the interior region functional, the local-density ap-proximation 共LDA兲共Ref. 2兲 with the Perdew-Wang 共PW兲

共Ref.12兲 parametrization of the Ceperly-Alder data13_for

cor-relation is used,⑀cLDA共n兲=⑀cPW共n兲:

⑀xcinterior_{关n兴 =}_⑀xLDA_{共n兲 +}_⑀cLDA_共n兲. _共1兲

As edge region functional, an exchange derived from the Airy Gas system5,14关the local Airy approximation 共LAA兲兴 is used together with a scaled PW LDA correlation:

⑀xcedge_{关n兴 =}_⑀xLAA_{共n,s兲 +}_␥⑀cLDA_共n兲 _␥_{= 0.8098,} _共2兲

where

s = s共n,兩ⵜn兩兲 = 兩ⵜn兩

2kF共n兲n

, kF共n兲 = 共3␲2n兲1/3 共3兲

is a dimensionless scaled gradient also used by many other GGA-type functionals, and we have suppressed the r depen-dence of the density; n = n共r兲.

The interpolation index,

X共s兲 = 1

1 +␣s2, ␣= 2.804, 共4兲

is only dependent on s, and determines the ration of edge vs interior functional to use in each point: X共0兲=1, X共s兲→0 when s→⬁, and

(3)

⑀xcAM05_{共n,s兲 =}_⑀xcinterior_{共n兲X共s兲 +}_⑀xcedge_{共n,s兲关1 − X共s兲兴. 共5兲}

Surface physics is included in AM05 by the use of the Airy gas system, and by the constants␣and␥determined from a fit to the jellium surface15 _RPA+ _共Ref. ₁₆_兲

exchange-correlation energies that are close to the recently published, to date most accurate, inhomogeneous Singwi-Tosi-Land-Sjölander energies.17

The nonspin-polarized total AM05 exchange-correlation energy, ExcAM05, is

Exc

AM05

=

冕

n共r兲⑀xcAM05„n共r兲,s共r兲…dr, 共6兲 =

冕

fxcAM05„n共r兲,s共r兲…dr. 共7兲

We note that the functional, or the exchange-correlation en-ergy per particle,⑀xc, in each point is multiplied by the den-sity in this point to obtain the exchange-correlation energy

density, fxc.

The exchange-correlation energy densities obtained from the separate subsystem functionals in Eqs. 共1兲 and 共2兲 can

unambiguously be made spin-resolved by using the spin scaling for the total exchange energy:

Ex关n↑,n↓兴 =

1

2共Ex关2n↑兴 + Ex关2n↓兴兲, 共8兲 and replacing the nonspin PW LDA correlation density with its spin-resolved form, fcLDA共n↑, n↓兲=共n↑+ n↓兲⑀cLDA共n↑, n↓兲.

However, a slight ambiguity is introduced by the interpola-tion index, a unique feature of AM05, which necessitate a more elaborate exploration than is usual for GGAs in general.

If one lets the mixing between the spin-resolved versions of the interior and edge subsystem functionals still be deter-mined by the total density and its scaled gradient via s = s共n,兩ⵜn兩兲 and X共s兲, as in Eq. 共5兲, one obtains a version we

will refer to as xtctt. This and other upcoming acronyms are constructed to indicate whether the t_{គotal or separate sគpin} densities determine the interpolation index for ex_{គchange and} c_{គorrelation while the last letter indicates if the translation} from correlation energy per particle to correlation energy density is via total or spin densities. The xtctt thus denotes that the t_{គotal density determines the interpolation for both} ex_{គchange and cគorrelation, and the correlation energy per} par-ticle is multiplied by the total density to obtain the correla-tion energy density:

fxc xtctt_共n ↑,n↓,s↑,s↓,s兲 =

再

1 2关fx LDA_共2n ↑兲 + fxLDA共2n↓兲兴 + fLDAc 共n↑,n↓兲

冎

X共s兲 +

再

1 2关fx LAA_共2n ↑,s↑兲 + fxLAA共2n↓,s↓兲兴 +␥fcLDA共n↑,n↓兲

冎

⫻关1 − X共s兲兴, 共9兲

where f_xLDA共n兲=n⑀xLDA共n兲, f_xLAA共n,s兲=n⑀xLAA共n,s兲, and, using Eq. 共3兲, we have defined

s_␯= s共2n_␯,兩ⵜ共2n_␯兲兩兲 ␯=↑,↓. 共10兲 However, this version leads to a total AM05 exchange energy that does not obey the exact exchange spin-scaling relation, Eq. 共8兲.

A version of spin-polarized AM05 that obeys the ex-change spin-scaling relation can be created by instead using separate indices for the spin-up 关X共s_↑兲兴 and spin-down 关X共s↓兲兴 exchange energy densities. This is the xsctt version

共exគchange interpolation is determined from the separate sគpin densities兲: fxc xsctt_共n ↑,n↓,s↑,s↓,s兲 = 1 2关fx LDA_共2n ↑兲X共s↑兲 + fx LDA_共2n ↓兲X共s↓兲兴 + fc LDA_共n ↑,n↓兲X共s兲 +1 2兵fx LAA_共2n ↑,s↑兲关1 − X共s↑兲兴 + fx LAA_共2n ↓,s↓兲关1 − X共s↓兲兴其 +␥fc LDA_共n ↑,n↓兲关1 − X共s兲兴. 共11兲

The xsctt version of AM05 fulfills the spin-scaling relation, Eq.共8兲, but violates a fundamental principle of the subsystem

functional scheme, the compatibility between exchange and correlation. Compatibility in this context requires that in each point in a system the exchange energy density needs to be combined with a correlation energy density based on the same model system as exchange. Since the amount of interior vs edge exchange and correlation in the xsctt version is determined with different indices, the total contribution in a point is not derived from the same model system. To restore the compatibility we need to use the same indices we use for exchange also for correlation. This version is xscst:

关1 − X共s↓兲兴. 共12兲

(4)

Since the two separate interpolation indices used in xscst are constructed from the respective spin-up and spin-down densities, the extent to which the physics in a specific point is interpreted as edgelike or interiorlike is in the general case different for the two spin channels, i.e., X共s_↑兲⫽X共s_↓兲. A physical motivation for such separate environments over the joint treatment in xtctt is seen in, e.g., a strongly magnetic situation where one should expect the surface that confines spin-up and spin-down electrons to be placed differently.

The separate indices introduce a correlation energy per

particle that is not the same for the spin-up and spin-down

electrons, reflecting the possibility that up and spin-down electrons in a point actually can experience different physical environments: the spin-up electrons could be at a surface while the spin-down electrons could be freely mov-ing as in an interior region. This is different from the case of most other functionals where the correlation energy per par-ticle is usually determined by the total density and the rela-tive spin polarization in a point, resulting in spin-up and spin-down particles contributing the same correlation energy per particle to the correlation energy density and the total correlation energy. The compatibility foundation of the sub-system functional scheme thus makes it natural to instead translate the correlation energy per particle into a correlation energy density by multiplying the separate up and spin-down parts of the correlation energy per particle with the spin-up/spin-down densities, respectively. We call this ver-sion for xscss: fxc xscss_共n ↑,n↓,s↑,s↓兲 =关n_↑⑀xLDA共2n_↑兲 + n_↑⑀cLDA共n_↑,n_↓兲兴X共s_↑兲 +关n_↑⑀xLAA共2n_↑,s_↑兲 + n_↑␥⑀cLDA共n_↑,n_↓兲兴关1 − X共s_↑兲兴 +关n_↓⑀xLDA共2n_↓兲 + n_↓⑀cLDA共n_↑,n_↓兲兴X共s_↓兲 +关n_↓⑀xLAA共2n_↓,s_↓兲 + n_↓␥⑀cLDA共n_↑,n_↓兲兴关1 − X共s_↓兲兴. 共13兲 By considering a fully spin-polarized system, that is, a system where one of the spin densities is zero everywhere, we can further clarify this argument. In such a system the total correlation energy should be determined only by the nonzero spin density profile and the nonzero spin density based interpolation. However, since the correlation energy per particle is nonvanishing also for the nonexisting elec-trons, in xscst 关Eq. 共12兲兴, the multiplication with the total

density will give a nondesired contribution to the total cor-relation energy while in xscss 关Eq. 共13兲兴 the multiplication

with the zero density will result in no contribution. The xscss version is the final version of the AM05 spin functional and the arguments based on fundamental subsystem functional principles given here for this choice over other possible choices resolve the apparent ambiguity of how to handle the interpolation index in a spin-resolved version of AM05.

As seen in Eq.共13兲, the spin AM05 functional consists of

a subsystem functional for interior regions and a subsystem functional for edge regions. The up density and spin-down density form separate density profiles, and separate indices are used for the separate density profiles to determine

the ratio of interior vs edge functional to use.

Using that ⑀xLAA共n,s兲=⑀xLDA共n兲FxLAA共s兲共see Ref. 5兲, the

AM05 exchange-correlation energy density in Eq. 共13兲 can

be written on the form

fxcAM05共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩兲 = n_↑⑀xLDA共2n_↑兲Hx共s_↑兲 + n_↓⑀xLDA共2n_↓兲Hx共s_↓兲 + n↑⑀cLDA共n_↑,n_↓兲Hc共s↑兲 + n↓⑀cLDA共n↑,n↓兲Hc共s↓兲, 共14兲 where Hx共s兲 = X共s兲 + 关1 − X共s兲兴Fx LAA_共s兲, _共15兲 Hc共s兲 = X共s兲 + 关1 − X共s兲兴␥. 共16兲

III. IMPLEMENTING AM05

It is straightforward to implement the AM05 exchange-correlation energy directly from Ref.5and Eq.共14兲, the only

minor obstacle being the need for a subroutine for calculat-ing the Lambert W function18_{used in F}

x

LAA_{共s兲. At the AM05}

website19 _{a subroutine for calculating the AM05}

exchange-correlation energy is provided共including subroutines for the required Lambert W function and the LSDA correlation in the parametrization of Perdew and Wang12,13_{兲. The 35 line}

subroutine for the Lambert W function is also provided as supplemental material to this paper.20

In this section and Appendix B, however, we will explic-itly give all formulas needed for implementing AM05, and for understanding the subroutines we provide at the AM05 website.19 _{In particular, we will give all necessary formulas}

for the exchange-correlation potentials needed for self-consistent DFT implementations.

The most elaborate part of an AM05 implementation is the LAA refinement factor, F_xLAA共s兲. From its definition in Eq. 8 in Ref.5 we obtain

3 4

冊

1/3 Z共s兲, 共20兲 Z共s兲 = s

冋

W„␹共s兲… ␹共s兲

册

2/3 = constant⫻ 关W„␹共s兲…兴2/3, 共21兲

(5)

␹共s兲 = s3/2

2

冑

6, 共22兲

where c = 0.7168, and Z共s兲 is a normalized function 共it is easy to verify that Z共s兲/s→1 when s→0兲 defined by Eq. 共21兲

which contains the Lambert W function and that it is conve-nient to use in implementations of AM05.

However, in self-consistent DFT calculations, not only is the exchange-correlation energy density needed but also the exchange-correlation potential for the spin-up electrons, the functional derivative of the exchange-correlation energy with respect to spin-up density, and its spin-down counter part:

Vxc,␯= ␦Exc

␦n_␯ ␯=↑,↓. 共23兲

There are several different schemes for calculating the exchange-correlation potentials within a code, and in this section we will give the formulas needed for the White and

Bird21_{scheme, and the fully assembled exchange-correlation}

potentials共that we call the traditional scheme兲. A short over-view of the different schemes for calculating the exchange-correlation potentials, including information on how to trans-fer the White and Bird expressions here given into the input needed for the Pople, Gill, and Johnson22 scheme, and de-tailed derivations, applicable also for other GGA-type func-tionals, are given in Appendix B.

Generally for a GGA-type functional the White and Bird scheme needs the five quantities ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲

⳵n↑ , ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵n↓ , ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn↑兩 , ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn↓兩 , and ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn兩 . Since AM05

does not use 兩ⵜn兩, we immediately see that

⳵fxc

AM05_共n

↑,n↑,兩ⵜn↑兩,兩ⵜn↓兩兲

⳵兩ⵜn兩 = 0. 共24兲

The remaining four quantities are

⳵fxc AM05_共n ↑,n↑,兩ⵜn↑兩,兩ⵜn↓兩兲 ⳵n_␯ a =vxLDA共2n␯_a兲Hx共s_␯_a兲 − 4 3⑀x LDA_共2n ␯a兲s␯a ⳵Hx共s␯_a兲 ⳵s_␯ a +v_c,_␯ a LDA_共n ↑,n↓兲 1 共n_↑+ n_↓兲关n↑Hc共s↑兲 + n↓Hc共s↓兲兴 +⑀cLDA共n_↑,n_↓兲 n␯b 共n_↑+ n_↓兲关Hc共s␯a兲 − Hc共s␯b兲兴 − 4 3⑀c LDA_共n ↑,n↓兲s␯a ⳵Hc共s␯_a兲 ⳵s_␯ a , 共25兲 and ⳵fxcAM05共n↑,n↑,兩ⵜn↑兩,兩ⵜn↓兩兲 ⳵兩ⵜn␯a兩 =⑀x LDA_共2n ␯a兲 2kF,␯_a ⳵Hx共s␯_a兲 ⳵s_␯ a +⑀c LDA_共n ↑,n↓兲 2kF,␯_a ⳵Hc共s␯_a兲 ⳵s_␯ a , 共26兲

where␯a=↑ or ↓, and␯b denotes the opposite spin to␯a. We have defined kF,␯= kF共2n␯兲, and vx

LDA_共2n

␯兲 is the LDA exchange

potential, andv_c,LDA_␯ 共n_↑, n_↓兲 is the LDA correlation potential for the␯spin. The explicit formulas for the derivatives of Hx共s兲 and Hc共s兲 are given in Appendix B.

In the traditional scheme, in addition to s_␯,␯=↑ ,↓, given in Eq. 共10兲, four other dimensionless density derivatives are used:

t_␯= ⵜ 2_共2n ␯兲 共2kF,␯兲2共2n␯兲 and u_␯=ⵜ共2n␯兲 · ⵜ兩ⵜ共2n␯兲兩 共2kF,␯兲3共2n␯兲2 . 共27兲

The full explicit exchange-correlation potentials are

V_xc,_␯ a AM05 =vx LDA_共2n ␯a兲

冋

Hx共s␯a兲 − s␯a ⳵Hx共s␯_a兲 ⳵s_␯ a

册

+⑀xLDA共2n_␯ a兲

再

冉

4 3s␯a 2 − t_␯ a

冊

冊

⳵ ⳵s_␯ a

冋

1 s_␯ a ⳵Hc共s_␯_a兲 ⳵s_␯ a

册

冎

+关⑀cLDA共n_↑,n_↓兲 − v_c,_␯ a LDA_共n ↑,n↓兲兴 n_␯ b n_↑+ n_↓

冋

Hc共s␯a兲 − Hc共s␯b兲 − s␯a ⳵Hc共s␯_a兲 ⳵s_␯ a

册

+关⑀cLDA共n_↑,n_↓兲 − v_c,_␯ b LDA_共n ↑,n↓兲兴ⵜn↑ ·ⵜn_↓ 兩ⵜn_␯_a兩2 n_␯ a n_↑+ n_↓s␯a ⳵Hc共s␯_a兲 ⳵s_␯ a . 共28兲

(6)

Again,␯a=↑ or ↓, and␯bdenotes the opposite spin to␯a. A detailed derivation of Eq.共28兲, and formulas for the

deriva-tives of Hx共s兲 and Hc共s兲, are given in Appendix B. Note that

even though the AM05 functional itself is not dependent on the total density gradient, this gradient is indeed present in the full correlation potential via the quantity ⵜn_↑·ⵜn_↓ =共兩ⵜn兩2₋_兩ⵜn

↑兩2−兩ⵜn↓兩2兲/2. This have no bearing on the

implementation of AM05 since the magnitude of the gradient of the total density is used in ordinary PBE共Ref.9兲

correla-tion while the magnitude of the gradients of the separate spin-up and spin-down densities are used in PBE exchange. As is the case with AM05 potentials obtained by the White and Bird, and the Pople, Gill, and Johnson schemes, AM05 potentials via the traditional scheme can thus be readily implemented in every code already containing PBE. Note that if exchange and correlation are treated separately in a code, it could be more convenient to include the last part of second line in Eq.共28兲 into the exchange part than the

cor-relation part since ordinary GGA exchange关see last part of first line in Eq.共28兲兴 already uses the needed density

deriva-tive quantities t_␯ and u_␯.

IV. AM05 RESULTS FOR Fe

We have investigated the spin version of AM05 by per-forming nonspinresolved and spin-resolved calculations for

the bcc and fcc phases of iron, usingVASP.23_{As is described}

in Refs. 6 and 24, in VASP 5 existing projector augmented wave共PAW兲 potentials can be used together with functionals they are not created with. Note that this is not generally true for other implementations of the PAW potentials and other types of pseudopotentials but must be tested thoroughly from case to case. Further details of our calculations are given in Appendix A.

Figures 1 and 2 show that minimal differences are ob-tained in AM05 results using LDA or PBE PAW potentials. It is clearly seen that AM05 obtains the correct ground state, the ferromagnetic bcc phase. Interestingly also the differ-ences between the four different spin versions are very small, as seen in Figs.3and4. In particular, there are no differences in lattice constants, bulk moduli, and magnetizations, in the bcc ground state, see TableI. The main difference is instead the energy difference between states with widely different magnetizations, such as nonmagnetic fcc and magnetic bcc, a property that is hard to derive from experiments. This will be investigated and discussed in future publications.

10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom bcc FM fcc NM fcc FM bcc NM 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom bcc FM fcc NM fcc FM bcc NM Iron: AM05xscss

on PBEfull and LDA dashed PAWpotentials

FIG. 1. 共Color online兲 Energy vs volume for non-magnetic 共NM兲 and ferro-magnetic 共FM兲 phases of bcc and fcc iron, calcu-lated with AM05.

9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom _{fcc FM} _{bcc FM} 9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom _{fcc FM} _{bcc FM} Iron: AM05xscss

on PBEfull and LDA dashed PAW potentials

FIG. 2. 共Color online兲 Magnetization vs volume for the FM phases of bcc and fcc iron, calculated with AM05.

10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom bcc FM fcc NM fcc FM bcc NM 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom bcc FM fcc NM fcc FM bcc NM 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom bcc FM fcc NM fcc FM bcc NM 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom bcc FM fcc NM fcc FM bcc NM

Iron: Comparison between AM05xscss and xsctt, xscst, xtctt, from bottom to top on PBE PAW potentials

FIG. 3.共Color online兲 Energy vs volume for NM and FM phases of bcc and fcc iron, calculated with the AM05 functional and three intermediate spin versions. All spin versions give the same results when the magnetization is zero, when they all collapse to the unique nonspin version of AM05.

9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom fcc FM bcc FM 9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom fcc FM bcc FM 9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom fcc FM bcc FM 9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom fcc FM bcc FM Iron: Comparison between AM05xscss and xsctt, xscst, xtctt,

from left to right on PBE PAW potentials

FIG. 4. 共Color online兲 Magnetization vs volume for the FM phases of bcc and fcc iron, calculated with the AM05 functional and three intermediate spin versions.

(7)

Finally, in Figs. 5 and 6 we compare the AM05 results with those obtained with LDA and the PBE 共Ref. 9兲

func-tionals. AM05 gives the same correct ferromagnetic bcc ground state as PBE while LDA gives the wrong spin state. Thus, the inclusion of surface physics into AM05 corrects one of the largest deficiencies of LDA. However, the PBE lattice constant, bulk modulus, and magnetization are closer to the experimental values than the corresponding AM05 re-sults, see TableI. This observation is further discussed in the next section.

V. CONCLUSIONS

We have discussed several possible versions of the spin density generalized AM05 functional. One of these versions,

xscss, is found to extend the underlying subsystem functional

framework of AM05 in a natural way, and is therefore pre-ferred. As a benchmark calculation we have chosen the ground state of Fe. The different versions of spin AM05 per-form very similarly and all give the correct spin ground state. Inclusion of surface physics into AM05 thus corrects a major failure of LDA, which gives the wrong spin state for this system. However, the lattice constants are not in as good agreement with experiment as PBE. Since nonspin AM05 generally is more accurate than PBE for nonmagnetic solid-state systems,6,8 _{there are two possible reasons for this. It}

might be that Fe just is a problematic case for AM05. But another possibility is that the excellent performance of AM05 for nonspin-polarized systems does not directly gen-eralize to a similar performance for spin-resolved cases when using spin-AM05. Either due to that the underlying spin for-mulation of the more approximate correlation part of the Airy Gas based edge functional in Eq.共2兲 is not adequate. Or

that magnetic materials are of a different class than nonmag-netic materials 共containing more localized electrons兲, which the interior共uniform electron gas based兲 and edge 共Airy gas surface system兲 subsystem functionals included in AM05 cannot adequately describe. The lattice constant and bulk moduli results of Ropo et al.7 _{indicate that the deficiency is}

not only for Fe but more general for 3d metals. If so, the minimal differences in results between the different versions of AM05 discussed in Sec.II, shown in Figs. 3 and4, and TableI, indicate that a change in the correlation of the edge

functional in Eq. 共2兲 would not result in a large change in

lattice constant, and thus not cure this deficiency. However, other properties than lattice constants, bulk moduli, and mag-netizations should also be considered, in particular properties that more closely probe differences in energy between non-spin and non-spin phases of solids. The most striking difference between the PBE and AM05 results shown in Fig. 5 is the slope of the spin-resolved fcc curve, and properties probing this should be investigated. Future applications of spin-AM05 will show which of the scenarios discussed above is the correct one. If AM05 turns out to generally be less accu-rate for spin-polarized systems than for nonspin-polarized systems, the next functional constructed according to the subsystem functional scheme should address this deficiency. TABLE I. Lattice constants, bulk moduli, and magnetizations,

for the ground-state bcc phase, obtained with the spin version of AM05, and the three intermediate versions. Experimental values are taken from Ref.25.

Version Lattice constant 共Å兲 Bulk modulus 共GPa兲 Magnetization 共␮B兲 xtctt 2.783 222 2.127 xscst 2.784 220 2.140 xsctt 2.785 220 2.143 AM05共xscss兲 2.786 218 2.148 Experiment 2.86 168 2.22 PBE 2.83 185 2.20 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom Iron: PBE bcc FM fcc NM fcc FM bcc NM 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom Iron: AM05xscss on PBE PAW potentials

bcc FM fcc NM fcc FM bcc NM 10 11 12 13 Volume Å3_atom 0.1 0.1 0.2 0.3 0.4 0.5 0.6 Energy eVatom Iron: LDA bcc FM fcc NM fcc FM bcc NM

FIG. 5.共Color online兲 Energy vs volume for NM and FM phases of bcc and fcc iron, calculated using LDA, AM05, and PBE.

(8)

ACKNOWLEDGMENTS

We thank T. R. Mattsson for fruitful discussions. R.A. acknowledges support from the Alexander von Humboldt Foundation and the German-Israeli Foundation. Sandia is a multiprogram laboratory operated by Sandia Corporation, a Lockheed Martin Company, for the United States Depart-ment of Energy’s National Nuclear Security Administration under contract DE-AC04-94AL85000.

APPENDIX A: DETAILS OF THE CALCULATIONS (CF. REF.26)

We implemented the spin version of AM05 and the three intermediate versions into version 5.1.40 ofVASP, which

al-ready had nonspin AM05 implemented. We used cubic cells for both the bcc共two atoms兲 and fcc 共four atoms兲 phases, and 16⫻16⫻16 k points in the Monkhorst-Pack scheme.27_The

PAW potentials are standard, the PBE one is labeled “PAW-_PBE Fe 06Sep2000,” and the LDA one “PAW Fe 03Mar1998.” A cutoff energy of 350 eV was used in all calculations. For the spin-resolved calculations an initial magnetic moment of 4␮Bon each site was used. The lattice constants and bulk moduli presented in TableI, and the cor-responding PBE results are calculated by a fit of seven energy/volume pairs, spaced evenly⫾10% around the equi-librium volume, to the Murnaghan28 _{equation of state. The}

magnetization was subsequently calculated at the obtained minimum-energy point by a separate VASP 5calculation that also served as a control so that the energy obtained in this point indeed was lower than in the previously calculated points. Our PBE results compares well with the PW91 re-sults of Ref. 25, as do our PBE and xscst results for lattice constant and bulk moduli with those of Ropo et al.共Ref.7兲.

The calculations for some of the intermediate spin versions converged slowly in the range where the magnetization changes rapidly with lattice constant. The results in this re-gion are sensitive to initial conditions; this can be seen in Fig.4in that some of the curves are jagged. Even though the magnetization in this region might be less accurate than in other regions, the general trend is not affected.

APPENDIX B: IMPORTANT FORMULAS AND BACKGROUND INFORMATION FOR IMPLEMENTING

SPIN AM05 INTO DFT CODES

The exchange-correlation energy for a spin-polarized sys-tem can be written as

Exc=

冕

drfxc„n_↑共r兲,n_↓共r兲,ⵜn_↑共r兲,ⵜn_↓共r兲…, 共B1兲

where n_↑共r兲 and n_↓共r兲 are the densities of spin-up electrons and spin-down electrons at r, respectively, dr = drxdrydrz is

the volume element, and fxc= fx+ fc is the

exchange-correlation energy density. fxcis generally a functional of the

spin-up and spin-down densities, fxc= fxc关n↑共r兲,n↓共r兲兴. We

here restrict ourselves to the case where the functional de-pendency is only through the densities and gradients of the densities since this is the case for AM05 关see Eq. 共14兲兴 and

most other GGA-type functionals.

In this appendix we will explicitly describe how to obtain the potentials given in Eq. 共23兲, needed for self-consistent

calculations. In addition, in some codes, the stress tensor is calculated and a functional routine also needs to output some ingredients for this calculation. The stress tensor is

␴␣␤=␦␣␤Exc关n共r兲兴 +

兺

␰=↑,↓

冕

drVxc,␰共r兲

冏

⳵n_␰共␧ˆ,r兲 ⳵␧␣␤

冏

␧ˆ=0 −

兺

␰=↑,↓

冕

dr⳵fxc共n↑共r兲,n↓共r兲,ⵜn↑共r兲,ⵜn↓共r兲兲 ⳵

共

⳵n␰共r兲 ⳵r␤

兲

⳵n_␰共r兲 ⳵r_␣ . 共B2兲 The terms on the first line of this expression are handled as

9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom Iron: PBE fcc FM bcc FM 9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom Iron: AM05 xscss on PBE PAW potentials

fcc FM bcc FM 9 10 11 12 13 Volume Å3_atom 0.5 1.0 1.5 2.0 2.5 3.0 Magnetization ΜBatom Iron: LDA fcc FM bcc FM

FIG. 6. 共Color online兲 Magnetization vs volume for the FM phases of bcc and fcc iron, calculated using LDA, AM05, and PBE.

(9)

in any LDA calculation, only replacing Exc, Vxc,↑, and Vxc,↓

with the corresponding GGA quantities. They will not be further discussed here. The terms on the second line consti-tutes a GGA correction and we will explain how to obtain the needed functional related ingredients for this correction.

The exchange-correlation potential, Eq. 共23兲, for the

exchange-correlation energy in Eq.共B1兲 is Vxc,␯共r兲 =

⳵fxc„n↑共r兲,n↓共r兲,ⵜn↑共r兲,ⵜn↓共r兲… ⳵n_␯共r兲

−ⵜ ·⳵fxc„n↑共r兲,n↓共r兲,ⵜn↑共r兲,ⵜn↓共r兲…

⳵ⵜn␯共r兲 . 共B3兲

The first term in this potential is straightforward to obtain but the second term can be treated in several ways. In the tradi-tional scheme this term is expanded until derivatives of the densities and derivatives of fxcare separated. This means that

for a given density the full Vxc,␯can be calculated within the

functional routine. In contrast, the White and Bird, and Pople, Gill, and Johnson schemes only partially expand this term and Vxc,␯needs to be assembled in a routine outside of

the functional subroutine.

Let us focus on the second term in Eq.共B3兲 for the

mo-ment. It is a scalar product between two vector quantities,ⵜ and⳵fxc共n↑共r兲,n↓共r兲,ⵜn↑共r兲,ⵜn↓共r兲兲

⳵ⵜn␯共r兲 . In a plane-wave basis this second

treated in the same way in the White and Bird, and the Pople, Gill, and Johnson schemes, the exchange-correlation poten-tial 共or Fock matrices兲 needs to be assembled outside of a functional subroutine and that this subroutine, instead of the full V_xc,↑ and V_xc,↓, needs to output ⳵fxc共n↑,n_⳵n↓,ⵜn↑,ⵜn↓兲

↑ and

⳵fxc共n↑,n↓,ⵜn↑,ⵜn↓兲

⳵ⵜn↑ , and the corresponding spin-down quantities,

where we from now on will omit the spatial argument of the densities.

However, ⳵fxc共n↑,n↓,ⵜn↑,ⵜn↓兲

⳵ⵜn↑ is a vector and together with its

spin-down counterpart gives a total of six different scalar

quantities to calculate. Due to the fact that functionals of GGA type for symmetry reasons only depend on the gradi-ents of the densities through their absolute values, 兩ⵜn_↑兩, 兩ⵜn_↓兩, and 兩ⵜn兩, where n=n_↑+ n_↓, only three scalar quantities are actually needed. Note that since 兩ⵜn兩2₌_兩ⵜn

↑兩2+兩ⵜn↓兩2

+ 2ⵜn_↑·ⵜn_↓⫽兩ⵜn_↑兩2+兩ⵜn_↓兩2, the dependency on the two vector quantities ⵜn_↑ and ⵜn_↓, needs to be replaced by a dependency on three scalar quantities. It is customary to use 兩ⵜn_↑兩, 兩ⵜn_↓兩, and 兩ⵜn兩 in the traditional and the White and Bird schemes while兩ⵜn_↑兩2_,_兩ⵜn

↓兩2, andⵜn↑·ⵜn↓are used in

the Pople, Gill, and Johnson scheme. One can show that

⳵fxc共n↑,n↓,ⵜn↑,ⵜn↓兲 ⳵ⵜn_␯_a = 2⳵fxc共n↑,n↓,兩ⵜn↑兩 2_,_兩ⵜn ↓兩2,ⵜn↑·ⵜn↓兲 ⳵兩ⵜn␯a兩 2 ⵜ n␯a +⳵fxc共n↑,n↓,兩ⵜn↑兩 2_,兩ⵜn ↓兩2,ⵜn↑·ⵜn↓兲 ⳵ⵜn_↑·ⵜn_↓ ⵜ n␯b 共B6兲 =⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn␯a兩 ⵜn␯a 兩ⵜn␯a兩 +⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn兩 ⵜn 兩ⵜn兩, 共B7兲 where the first equality 关Eq. 共B6兲兴 gives the quantity that

is usually seen in the Pople, Gill, and Johnson scheme and the second equality 关Eq. 共B7兲兴 is seen in the White

and Bird scheme, and␯a=↑ or ↓, and␯bdenotes the opposite spin from ␯a. _兩ⵜnⵜn↑

↑兩,

ⵜn↓

兩ⵜn↓兩, and

ⵜn

兩ⵜn兩 are unit vectors that are

well defined even when 兩ⵜn_↑兩, 兩ⵜn_↓兩, or 兩ⵜn兩→0, respec-tively, and those factors can be handled elsewhere in the code, and ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲

⳵兩ⵜn↑兩 ,

⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲

⳵兩ⵜn↓兩 , and

⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲

⳵兩ⵜn兩 can be given out from the functional

sub-routine. However, there is no consensus in codes about this matter and care needs to be taken in order to make sure that the right quantities are output.

The AM05 subroutine provided at the AM05 website19

gives output for the White and Bird scheme:

⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵n↑ , ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵n↓ , ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn↑兩 , ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn↓兩 , and ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲

⳵兩ⵜn兩 . If instead the Pople, Gill, and Johnson

output is needed it can be obtained from these quantities by using that 2⳵fxc共n↑,n↓,兩ⵜn↑兩 2_,兩ⵜn ↓兩2,ⵜn↑·ⵜn↓兲 ⳵兩ⵜn_↑兩2 = 1 兩ⵜn↑兩 ⳵fxc共n_↑,n_↓,兩ⵜn_↑兩,兩ⵜn_↓兩,兩ⵜn兩兲 ⳵兩ⵜn↑兩 + 1 兩ⵜn兩 ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn兩 , 共B8兲

(10)

2⳵fxc共n↑,n↓,兩ⵜn↑兩 2_,兩ⵜn ↓兩2,ⵜn↑·ⵜn↓兲 ⳵兩ⵜn↓兩2 = 1 兩ⵜn↓兩 ⳵fxc共n_↑,n_↓,兩ⵜn_↑兩,兩ⵜn_↓兩,兩ⵜn兩兲 ⳵兩ⵜn↓兩 + 1 兩ⵜn兩 ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn兩 , 共B9兲 ⳵fxc共n↑,n↓,兩ⵜn↑兩2,兩ⵜn↓兩2,ⵜn↑·ⵜn↓兲 ⳵ⵜn↑·ⵜn↓ = 1 兩ⵜn兩 ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn兩 . 共B10兲

This same output can also be used in calculations of the GGA correction for the stress tensor in Eq.共B2兲. The terms

on the second line of Eq. 共B2兲 constitutes the GGA

correc-tion and it can be shown that

兺

␰=↑,↓ ⳵fxc共n_↑,n_↓,ⵜn_↑,ⵜn_↓兲 ⳵

共

⳵n␰ ⳵r␤

兲

⳵n_␰ ⳵r_␣ =

兺

␰=↑,↓2 ⳵fxc共n_↑,n_↓,兩ⵜn_↑兩2,兩ⵜn_↓兩2,ⵜn_↑·ⵜn_↓兲 ⳵兩ⵜn␰兩2 ⳵n_␰ ⳵r_␤ ⳵n_␰ ⳵r_␣ +⳵fxc共n↑,n↓,兩ⵜn↑兩 2_,兩ⵜn ↓兩2,ⵜn↑·ⵜn↓兲 ⳵ⵜn↑·ⵜn↓ ⫻

冉

⳵n↑ ⳵r_␤ ⳵n_↓ ⳵r_␣+ ⳵n_↓ ⳵r_␤ ⳵n_↑ ⳵r_␣

冊

共B11兲 =

兺

␩=↑,↓,0 ⳵fxc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 ⳵兩ⵜn_␩兩 1 兩ⵜn_␩兩 ⳵n_␩ ⳵r_␤ ⳵n_␩ ⳵r_␣, 共B12兲 where n0= n = n↑+ n↓, and we again have given the Pople,

Gill, and Johnson 关Eq. 共B11兲兴, and the White and Bird 关Eq.

共B12兲兴 relevant formulas. We recognize the same derivatives

of fxcas appear in Eqs.共B6兲 and 共B7兲.

In the traditional scheme the second term in Eq.共B3兲 is

further expanded using Eq.共B7兲:

. 共B13兲 It is straightforward to show that

ⵜ · ⵜn␩ 兩ⵜn_␩兩= ⵜ2_n ␩ 兩ⵜn_␩兩− ⵜn_␩·ⵜ兩ⵜn_␩兩 兩ⵜn_␩兩2 for ␩=↑,↓,0, 共B14兲 and the two last terms in Eq. 共B13兲 can also be further

ex-panded by using that ⵜg共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩,兩ⵜn兩兲 =ⵜn_↑⳵g ⳵n_↑+ⵜn↓ ⳵g ⳵n_↓+ⵜ兩ⵜn↑兩 ⳵g ⳵兩ⵜn_↑兩+ⵜ兩ⵜn↓兩 ⳵g ⳵兩ⵜn_↓兩 +ⵜ兩ⵜn兩 ⳵g ⳵兩ⵜn兩. 共B15兲

However, the final general formula would be very cumber-some and we thus now leave generality and instead make use of specific dependencies existing in GGA-type functionals, focusing on the AM05 functional.

The exchange terms in AM05 are only dependent on ei-ther spin-up or spin-down densities. This is also the case for GGA exchange in general. The dependency of the correlation part of the spin-dependent AM05 functional is fc

= fc共n_↑, n_↓,兩ⵜn_↑兩兲+ fc共n_↑, n_↓,兩ⵜn_↓兩兲 while, for example, PBE

ⵜn_␩·ⵜn_␰ 兩ⵜn_␩兩 ⳵2_f_共n ↑,n↓,兩ⵜn␩兩兲 ⳵n_␰⳵兩ⵜn_␩兩 +ⵜn_␩·ⵜ兩ⵜn_␩兩 ⳵ ⳵兩ⵜn_␩兩

冉

1 兩ⵜn_␩兩 ⳵f共n_↑,n_↓,兩ⵜn_␩兩兲 ⳵兩ⵜn_␩兩

冊

, 共B16兲

(11)

when␩=␯or 0, and 0 if␩is the opposite spin to␯. We have here suppressed the index denoting exchange and correla-tion. For the exchange, one of the terms in the sum over spin up and spin down will vanish since the exchange part of the functional only depends on one spin and thus the cross-spin derivative will be zero. The density gradient ⵜn_␩·ⵜn_␰ is a function only of 兩ⵜn_␩兩, ␩=↑, ↓,0, which can be made ex-plicit by using thatⵜn_↑·ⵜn_↓=共兩ⵜn兩2₋_兩ⵜn

↑兩2−兩ⵜn↓兩2兲/2.

We notice that derivatives of the functional and deriva-tives of the densities no longer are intermingled but sepa-rated and if the density derivativesⵜ2n_␩andⵜn_␩·ⵜ兩ⵜn_␩兩 are handled into a functional routine in addition to the already required␩_↑,␩_↓, and兩ⵜn_␩兩, the two exchange-correlation po-tentials can be assembled. Note that even if the quantity ⵜn␩·ⵜ兩ⵜn␩兩 looks complicated it is easily shown that

ⵜn_␩·ⵜ兩ⵜn_␩兩 = 1 兩ⵜn_␩兩

兺

i=1 3

兺

j=1 3 ⳵n_␩ ⳵ri ⳵n_␩ ⳵rj ⳵2_n ␩ ⳵ri⳵rj , 共B17兲 where r₁= rx, r2= ry, and r3= rz.

In the following more detailed manipulations, we will only treat AM05, leaving, for example, the derivation of the PBE correlation potentials as an exercise. For ease we divide up the AM05 exchange-correlation energy density in Eq. 共14兲 in separate exchange and correlation parts:

fx共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩兲 =1 2关fx LDA_共2n ↑兲Hx共s_↑兲 + fx LDA_共2n ↓兲Hx共s_↓兲兴, 共B18兲

where fxLDA共n兲=n⑀xLDA共n兲 is the LDA exchange energy

den-sity, and fc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩兲 = fc LDA_共n ↑,n↓兲

冉

n_↑ 共n↑+ n↓兲Hc共s↑兲 + n_↓ 共n↑+ n↓兲Hc共s↓兲

冊

, 共B19兲 where fcLDA共n↑, n↓兲=共n↑+ n↓兲⑀cLDA共n↑, n↓兲 is the LDA

correla-tion energy density. Using that ⳵s共n,ⵜn兲 ⳵n = − 4 3 s n, 共B20兲 ⳵s共n,ⵜn兲 ⳵兩ⵜn兩 = 1 2kFn , 共B21兲

we obtain from Eq. 共B18兲

⳵fx共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩兲 ⳵n_␯ =⳵fx LDA_共2n ␯兲 ⳵共2n␯兲 Hx共s␯兲 − fx LDA_共2n ␯兲4 3 s_␯ 2n_␯ ⳵Hx共s␯兲 ⳵s_␯ =vx LDA_共2n ␯兲Hx共s_␯兲 − 4 3⑀x LDA_共2n ␯兲s␯⳵ Hx共s␯兲 ⳵s_␯ , 共B22兲 and ⳵fx共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩兲 ⳵兩ⵜn␯兩 = fx LDA_共2n ␯兲 1 2kF,␯共2n␯兲 ⳵Hx共s␯兲 ⳵s_␯ =⑀x LDA_共2n ␯兲 2kF,␯ ⳵Hx共s␯兲 ⳵s_␯ , 共B23兲 where k_F,_␯= kF共2n␯兲, vx LDA_共2n

␯兲 is the LDA spin-␯ exchange

+ fc LDA_共n ↑,n↓兲 2n_␯ a 共n_↑+ n_↓兲 ⳵Hc共s␯_a兲 ⳵共2n␯a兲 =v_c,_␯ a LDA_共n ↑,n↓兲 1 共n_↑+ n_↓兲关n↑Hc共s↑兲 + n↓Hc共s↓兲兴 +⑀c LDA_共n ↑,n↓兲 n_␯ b 共n_↑+ n_↓兲关Hc共s␯a兲 − Hc共s␯b兲兴 −4 3⑀c LDA_共n ↑,n↓兲s␯a ⳵Hc共s␯_a兲 ⳵s_␯ a , 共B24兲

where␯bis the opposite spin to␯a, and

⳵fc共n↑,n↓,兩ⵜn↑兩,兩ⵜn↓兩兲 ⳵兩ⵜn_␯兩 = fcLDA共n↑,n↓兲 2n_␯ 共n_↑+ n_↓兲 ⳵Hc共s_␯兲 ⳵兩ⵜ共2n_␯兲兩= ⑀cLDA_共n ↑,n↓兲 2kF,␯ ⳵Hc共s_␯兲 ⳵s_␯ , 共B25兲

(12)

where v_c,LDA_␯ 共n_↑, n_↓兲 is the LDA spin-␯ correlation potential and⑀cLDA共n_↑, n_↓兲 is the LDA correlation energy per particle.

The quantities in Eqs. 共B22兲–共B25兲 gives the input

needed; ⳵fxc ⳵n␯a = ⳵fx ⳵n␯a + ⳵fc ⳵n␯a

, giving Eq. 共25兲, and ⳵fxc

⳵兩ⵜn␯a兩 = ⳵fx

⳵兩ⵜn␯a兩 + ⳵fc

⳵兩ⵜn␯a兩

, giving Eq.共26兲, in the White and Bird and the Pople,

Gill, and Johnson schemes for constructing the exchange-correlation spin-␯apotential outside of the functional subrou-tine 关see, however, the note below Eqs. 共B6兲 and 共B7兲兴.

For the traditional scheme, starting with exchange, we obtain from Eq.共B16兲, Eqs. 共B22兲 and 共B23兲, and the

defi-nitions in Eqs.共10兲 and 共27兲:

+ s ⳵H共s兲 ⳵s , 共B27兲 and, from Eqs.共B3兲, 共B22兲, and 共B26兲, we arrive at the final

form of the exchange part of the spin-␯exchange-correlation potential: Vx,␯=vxLDA共2n␯兲

冋

Hx共s_␯兲 − s_␯ ⳵Hx共s␯兲 ⳵s_␯

册

+⑀x LDA_共2n ␯兲 ⫻

再

冉

4 3s␯ 2 − t_␯

冊

1 s_␯ ⳵Hx共s␯兲 ⳵s_␯ +

冉

4 3s␯ 3 − u_␯

冊

⳵ ⳵s_␯

冋

1 s_␯ ⳵Hx共s␯兲 ⳵s_␯

册

冎

. 共B28兲

For correlation, by a similar procedure starting from Eqs. 共B24兲 and 共B25兲 and again using Eq. 共B16兲, we get

. 共B29兲

Using Eqs.共B3兲, 共B24兲, 共B27兲, and 共B9兲, we arrive at V_c,_␯ a=vc,␯a LDA_共n ↑,n↓兲

冋

Hc共s␯_a兲 − s␯_a ⳵Hc共s␯_a兲 ⳵s_␯ a

册

+⑀cLDA共n_↑,n_↓兲

再

冉

4 3s␯a 2 _{− t} ␯a

冊

1 s_␯ a ⳵Hc共s␯_a兲 ⳵s_␯ a +

冉

4 3s␯a 3 _{− u} ␯a

冊

⳵ ⳵s_␯ a

冋

1 s_␯ a ⳵Hc共s␯_a兲 ⳵s_␯ a

册

冎

+关⑀cLDA共n_↑,n_↓兲 − v_c,_␯ a LDA_共n ↑,n↓兲兴 n_␯ b n_↑+ n_↓

冋

Hc共s␯a兲 − Hc共s␯b兲 − s␯a ⳵Hc共s␯_a兲 ⳵s_␯ a

册

+关⑀cLDA共n_↑,n_↓兲 − v_c,_␯ b LDA_共n ↑,n↓兲兴ⵜn↑ ·ⵜn_↓ 兩ⵜn␯a兩 2 n_␯ a n_↑+ n_↓s␯a ⳵Hc共s␯_a兲 ⳵s_␯ a . 共B30兲

The total spin-␯aexchange-correlation potential, V_xc,_␯

a= Vx,␯a+ Vc,␯a, is given in Eq.共28兲.

For implementation into codes also the various s derivatives of Hx共s兲 and Hc共s兲 in Eqs. 共B22兲–共B25兲, 共B28兲, and 共B30兲 are

needed. We will here adhere closely to the notation used in the subroutines provided at the AM05 web page.19_{We have}

1 s ⳵Hx共s兲 ⳵s =关1 − X共s兲兴 1 s ⳵Fx LAA_共s兲 ⳵s +关1 − Fx LAA_共s兲兴1 s ⳵X共s兲 ⳵s , 共B31兲 1 s ⳵Hc共s兲 ⳵s =关1 −␥兴 1 s ⳵X共s兲 ⳵s , 共B32兲 and

(13)

=关1 −␥兴 1 s ⳵ ⳵s

冉

冊

, 共B34兲 where 1 s ⳵X共s兲 ⳵s = − 2␣关X共s兲兴 2_, _共B35兲 1 s ⳵ ⳵s

冉

冊

= 8␣ 2_{关X共s兲兴}3_, _共B36兲

and X共s兲 is given in Eq. 共4兲. The derivatives of FxLAA共s兲 are more elaborate. Using the definitions in Eqs. 共17兲 and 共18兲 we

obtain 1 s ⳵FxLAA共s兲 ⳵s = 1 s c2sD共s兲 − 共cs2_{+ 1}_兲dD共s兲 ds 关D共s兲兴2 = U共s兲 关D共s兲兴2= c 关D共s兲兴2

冋

2 −共1 + cs 2_兲dA共s兲 ds −共1 − cs 2_兲A共s兲 s

册

, 共B37兲

where A共s兲 is given in Eq. 共19兲 and dA共s兲/ds is given below, in Eq. 共B41兲. Lastly, using the definition of U共s兲 on the first line

of Eq.共B37兲, we find 1 s ⳵ ⳵s

Finally, the s derivatives of Z共s兲 can be derived from its definition in Eq. 共21兲 by using that

dW共x兲 dx = 1 x W共x兲 1 + W共x兲, 共B43兲 which together with the definition of ␹共s兲 in Eq. 共22兲 gives

dW„␹共s兲… ds = 3 2 1 s W„␹共s兲… 1 + W„␹共s兲…. 共B44兲 The final ingredients needed for the implementation of the AM05 potentials are thus

dZ共s兲 ds = Z共s兲 s 1 1 + W„␹共s兲…, 共B45兲 and, noting that the s times the second derivative of Z共s兲 at the first line in Eq.共B42兲 only will appear together with the