How close are the Slater and Becke-Roussel potentials in solids?

The Becke–Roussel (BR) potential [Phys. Rev. A 1989, 39, 3761] was proposed as an approximation to the Slater potential, which is the Coulomb potential generated by the exact exchange hole. In the present work, a detailed comparison between the Slater and BR potentials in solids is presented. It is shown that the two potentials usually lead to very similar results for the electronic structure; however, in a few cases, e.g., Si, Ge, or strongly correlated systems like NiO, the fundamental band gap or magnetic properties can differ markedly. Such differences should not be neglected when the computationally expensive Slater potential is replaced by the cheap semilocal BR potential in approximations to the exact-exchange Kohn–Sham potential, such as the one proposed by Becke and Johnson [J. Chem. Phys. 2006, 124, 221101].

Introduction

The Kohn–Sham (KS) version of density functional theory[1,2] represents a favorable compromise between accuracy and computational cost for the calculation of the structural and electronic properties of molecules and solids.[3−5] In particular, the use of a semilocal approximation for the exchange-correlation (xc) energy Exc = ∫εxc(r)d3r and the multiplicative KS xc potential vxc(r) = δExc/δρ(r) (rungs one to three of Jacob’s ladder[6]) allows one to treat systems containing thousands of atoms, which is out of reach for methods using nonlocal Hartree–Fock (HF) exchange and post-HF many-body methods (e.g., perturbation theories or random-phase approximation). The semilocal approximations lead to fast calculations since the evaluation of εxc(r) and vxc(r) at a point r requires the value of the electron density ρ(r), and eventually its derivatives and the kinetic-energy density t(r), only at that same point r. In the local density approximation[2,7,8] (LDA, first rung of Jacob’s ladder), εxc(r) depends only on ρ(r), whereas in the generalized gradient approximation[9,10] (GGA, second rung), εxc(r) depends additionally on ∇ρ(r). A further dependence on t(r) and/or ∇2ρ(r) leads to the meta-GGA functionals (see, e.g., ref (11)) of the third rung.

These semilocal methods lead to sufficiently accurate results in many circumstances, but this is certainly not systematically the case. For instance, the standard semilocal xc potentials (the focus of this work) are particularly inaccurate for systems with strongly correlated electrons such that the results are very often qualitatively wrong.[12] Improved results can be obtained by calculating the exchange part vx of the KS xc potential exactly (see, e.g., ref (13)), which can be done by solving the equations of the optimized effective potential method applied to the exact exchange energy (EXX-OEP).[14,15] However, this leads to calculations that are computationally expensive and prone to instabilities, in particular because of the use of unoccupied orbitals. Nevertheless, EXX-OEP can be very useful for the construction of more accurate semilocal approximations for the exchange potential. Accurate KS correlation potentials vc derived from many-body theories can also be calculated (see, e.g., refs (16) and (17)); however, such calculations are even more complex than for EXX.

Related to this, we will consider the exchange potential of Becke and Roussel[18] (BR) that was proposed as an approximation to the Slater potential.[19] The Slater potential, which is the Coulomb potential due to the HF exchange hole, is often used as the first term in (beyond semilocal) approximations to the EXX-OEP like the Krieger–Li–Iafrate approximation[20] (see ref (21) for a summary of such methods). The construction of the BR potential starts with the modeling of the exchange hole by using the electron density and its first two derivatives as well as the kinetic energy density.[18] This leads to a potential that belongs to the family of semilocal meta-GGA methods. Up to now, comparisons between the Slater and BR potentials have been done by Becke and co-workers[18,22] for spherical atoms (from He to Cd), Heßelmann and Manby[23] for molecules containing elements of the first and second periods, and Karolewski et al.[24] on the Be atom and C6H8 molecule. In these studies, it was shown that the BR potential is a fairly good approximation to the Slater potential. However, it is not clear if this conclusion would remain valid in the case of more complicated systems, e.g., molecules or solids with heavier elements or magnetic systems. A more thorough comparison between the Slater and BR potentials is also motivated by the extensive use of the BR potential (instead of the much more expensive Slater potential) as the first term in the exchange potential proposed by Becke and Johnson[22] (BJ), which was proposed as an approximation to the EXX-OEP in atoms. Furthermore, a semilocal and accurate approximation to the Slater potential could also lead to a cheap replacement to the nonlocal HF exchange energy for total energy calculations.[5,25−30] Our goal for the present work is to obtain more insight into the differences between the Slater and BR potentials in solids and their consequences on the calculated quantities like the KS fundamental band gap (defined as the conduction band minimum CBM minus the valence band maximum VBM) or magnetic moment.

The article is organized as follows. Section gives a brief presentation of the Slater and BR potentials. Then, the results are presented and discussed in Section , and a summary of the results is given in Section .

Methods and Computational Details

As a simplification of the HF method, Slater[19] (S) proposed replacing the HF exchange potential vx,HF, which is different for each orbital ψ, by a common orbital-averaged potential (|ψ|2/ρσ are the weights)where Nσ is the number of occupied orbitals for the spin σ. From the expression of the HF exchange energywhereis the HF exchange hole (the depletion of spin-σ density at r′ for a reference spin-σ electron at r due to the Pauli exclusion principle), we can see that the Slater potential can also be interpreted as the Coulomb potential generated by the HF exchange hole:The Slater potential has been used as the leading term in various approximations to the EXX-OEP (see refs (21 and 31−33) and references therein), whereas the remaining part vx,σresp = vx,σEXX-OEP – vx,σS is often called the response (resp) term. It is important to mention that the Slater potential does not reduce to the correct limit[2,7,8] for the homogeneous electron gas given by vx,σLDA = −(6ρσ/π)1/3, but to (3/2)vx,σLDA instead.[19] Note that α(3/2)vx,σLDA, where α is a parameter (see, e.g., refs (34) and (35)), is also known as the Slater (or Xα) potential.

It is clear from eqs and 2 that the calculation of the Slater potential is still rather expensive and leads to a computational cost that is basically the same as the HF exchange energy. Therefore, in ref (18), Becke and Roussel proposed a semilocal form for the potential generated by the exchange hole. The derivation begins with the modeling of the exchange hole with semilocal quantities, namely, ρ and its two derivatives and the positive-definite kinetic energy density tσ = (1/2)∑ ∇ψ*·∇ψ. Then, the potential generated by this approximate exchange hole is given byThe function xσ in eq can be calculated either by solving at each point of space the nonlinear equation[18]wherewithor by using the analytical interpolation formula for xσ proposed in ref (27). After xσ is calculated, bσ in eq is given byThe value of γ in eq should be 1 in principle; however, it was shown that the potential generated by the exchange hole in the homogeneous electron gas (i.e., (3/2)vx,σLDA) is recovered for γ = 0.8.[18] More recently, it has been reported that other values for γ may lead to better agreement with the Slater potential or the EXX-OEP (when combined with some response term vx,σresp).[23,36]

The BJ potential,[22] which readshas attracted considerable attention in recent years (see, e.g., refs (24, 31, and 36−44)). The second term in eq was proposed as a semilocal approximation to the response term vx,σresp and reduces to –(1/2)vx,σLDA at the limit of a constant electron density such that vx,σBJ reduces to vx,σLDA at this limit if γ = 0.8 in eq . In ref (22), it was shown that the BJ potential visually resembles very much the EXX-OEP in spherical atoms, independently of which potential (Slater or BR) was used for the first term. However, non-negligible differences between the EXX-OEP and BJ eigenvalue spectra have also been reported.[36,37,41,45] For large molecules and periodic solids, it is obviously much more advantageous to use BR in eq in order to make the BJ potential fully semilocal, as has been done, for instance, for the modified BJ potential that has been used for the calculation of band gaps in solids.[38,46−51] Other implementations of the BJ potential or one of its variants using the BR potential have been reported in refs (45 and 52−59), whereas the Slater potential was used for the works in refs (31 and 39−44). In the following, the acronyms BJS and BJBR will refer to eq with the Slater and BR potentials, respectively.

For the comparison of the Slater and BR potentials, we will consider the set of nonmetallic solids listed in Table along with their space group and geometrical parameters. The solids were chosen such that various types of bonding are represented: ionic (e.g., MgO), covalent (e.g., Si), and van der Waals (rare gases). MnO, FeO, CoO, and NiO are antiferromagnetic (the ferromagnetic planes are stacked along the [111] direction[60]), whereas all other solids are nonmagnetic. The differences between the Slater and BR potentials will be measured by comparing the results obtained for the electronic structure, total energy, electric-field gradient (EFG) in Se and Cu2O, and magnetic moment in MnO, FeO, CoO, and NiO.

Table 1

Space Group, Geometrical Parameters, and Core Electrons of the Solids Considered in This Worka

solid	space group	geometrical parameters	core states
Ne	Fm3m	a = 4.470	1s
Ar	Fm3m	a = 5.260	[Ne]
Kr	Fm3m	a = 5.598	[Ar]3d
Xe	Fm3m	a = 6.130	[Kr]
C	Fd3m	a = 3.568	1s
Si	Fd3m	a = 5.430	[Ne]
Ge	Fd3m	a = 5.652	[Ar]
Se	P3121	a = 4.366, c = 4.954,  xSe = 0.225	[Ar]
BN	F43m	a = 3.616	B: 1s; N: 1s
SiC	F43m	a = 4.358	Si: [Ne]; C: 1s
GaAs	F43m	a = 5.648	Ga: [Ar]; As: [Ar]
InP	F43m	a = 5.869	In: [Kr]; P: [Ne]
CdS	F43m	a = 5.818	Cd: [Ar]3d4s; S: [Ne]
LiH	Fm3m	a = 4.084
LiCl	Fm3m	a = 5.106	Cl: [Ne]
BeO	P63mc	a = 2.694, c = 4.384,  zO = 0.378	Be: 1s; O: 1s
MgO	Fm3m	a = 4.207	Mg: 1s; O: 1s
CsF	Fm3m	a = 6.030	Cs: [Kr]; F: 1s
BaO	Fm3m	a = 5.523	Ba: [Kr]4d; O: 1s
PbS	Fm3m	a = 5.936	Pb: [Xe]4f; S: [Ne]
ScN	Fm3m	a = 4.500	Sc: [Ne]; N: 1s
SrTiO3	Pm3m	a = 3.905	Sr: [Ar]3d; Ti: [Ne]; O: 1s
MnO	Fm3m, R3m	a = 4.445	Mn: [Ne]3s; O: 1s
FeO	Fm3m, R3m	a = 4.334	Fe: [Ne]3s; O: 1s
CoO	Fm3m, R3m	a = 4.254	Co: [Ne]3s; O: 1s
NiO	Fm3m, R3m	a = 4.171	Ni: [Ne]3s; O: 1s
ZnO	P63mc	a = 3.258, c = 5.220,  zO = 0.382	Zn: [Ar]; O: 1s
Cu2O	Pn3m	a = 4.267	Cu: [Ne]3s; O: 1s
CeO2	Fm3m	a = 5.411	Ce: [Kr]4d; O: 1s

The lattice parameters are in Å, and the internal parameters are in internal units. For MnO, FeO, CoO, and NiO, the antiferromagnetic order leads to a lowering of the symmetry (second indicated space group). The last column shows the electrons that were considered as core electrons, for which a fully relativistic treatment (i.e., spin-orbit coupling included) is used.

The Slater and BR potentials will be compared both with and without the second term of the BJ potential (see eq ). Furthermore, in order to give a better idea of how important the discrepancies between Slater and BR are, the results obtained with other exchange-only potentials will also be shown. These potentials are the LDA[2,7,8] and various GGAs, namely, the ones from Perdew et al.[10] (PBE), Engel and Vosko[61] (EV93), and Armiento and Kümmel[62] (AK13). Note that the derivative discontinuity Δxc,[3,63−65] which accounts for the difference between the KS band gap (CBM minus VBM) and the experimental band gap (ionization potential minus electron affinity), is zero for the LDA, PBE, and EV93 potentials, but it is nonzero for the others.[39,62] We also mention that in ref (58) the BJ and AK13 potentials were (rigidly) shifted such that they go to zero in the asymptotic region far from the nuclei. Such a shift has an effect on properties (e.g., ionization potential) that is calculated by using the absolute value of the eigenvalues and not just the differences between them. Such a shift would not change the results for the properties considered in our work.

The calculations were done with the all-electron code WIEN2k,[66] which is based on the full-potential linearized augmented plane-wave method,[67,68] to solve the KS equations. The implementation of the Slater potential has been done without any approximation, and the sums over the occupied orbitals in eq include both the band and core electrons. As was done previously for the implementation of the HF method (see ref (69)), the pseudocharge method,[70,71] combined with the technique from refs (72) and (73) to treat the Coulomb singularity, has been used to calculate the Slater potential. The computational parameters, like the size of the basis set, the number of k-points for the integrations of the Brillouin zone, or those specific to the calculation of the Slater potential, were chosen such that the results are well-converged. For instance, the band gaps should be converged within ∼0.02 eV. The core electrons (indicated in Table ) were treated fully relativistically (i.e., including spin-orbit coupling), whereas a scalar-relativistic treatment[74] was used for the valence electrons.

Results and Discussion

The calculations with the BR potential were done using γ = 0.8 in eq , which as mentioned above, is the value that leads to (3/2)vx,σLDA in the limit of a constant electron density (as with the Slater potential). The results obtained with other values for γ will be discussed later in the article. Before discussing the results, we mention that in solids with an open shell of strongly correlated electrons it is usually possible to stabilize more than one solution, which is the case for FeO and CoO among the solids in our test set. As has been already shown in ref (75) for FeO, EV93 leads to two solutions of nearly equal total energy. One solution has a nonzero fundamental band gap, whereas the other is metallic. For AK13, BJS, and BJBR(0.8) (as well as LDA+U), the solution with a gap is more stable (using the HF total energy for BJS and BJBR(0.8) orbitals). For these potentials, only the results for the insulating state will be considered for the discussion below. With the other potentials (LDA, PBE, Slater, and BR(0.8)), only the metallic state could be stabilized. In the case of CoO, the most stable state is insulating for all potentials except LDA and BR(0.8). For both FeO and CoO, we have observed that the occupation matrix of the 3d orbitals for the insulating state is rather similar among all potentials, but it is different from the occupation matrix of the metallic state.

The results that will be discussed in this work are shown in Tables –5 as well as graphically in Figures S1–S35 of the Supporting Information, except for those obtained with the Slater and BR potentials in order to use a reasonable scale for the results for the core states and total energy because the results obtained with these two potentials are very different from the others.

Table 2

Fundamental Band Gap (in eV) Calculated with Different Exchange-Only Potentialsa

solid	LDA	PBE	EV93	AK13	BR(0.8)	S	BJBR(0.8)	BJS	EXX-OEP
Ne	10.80	11.00	10.72	20.07	15.89	15.58	13.20	12.89	14.15b, 14.79c
Ar	7.82	8.38	8.95	15.16	9.99	9.95	9.21	9.23	9.61b, 9.65c
Kr	6.39	7.02	7.74	12.80	8.31	8.20	7.57	7.49	7.87b
Xe	5.45	6.07	6.82	10.68	7.03	6.88	6.40	6.34	6.69b
C	4.00	4.46	4.60	4.78	4.64	4.72	4.31	4.51	4.57d
Si	0.35	0.80	1.12	1.60	0.69	1.00	0.71	1.11	1.18d
Ge	0	0	0.47	0.71	0	0	0.25	0.67	0.89c
Se	0.97	1.12	1.35	1.56	0.99	1.08	1.17	1.26
BN	4.20	4.89	5.25	5.68	5.39	5.43	4.80	4.96	5.43d
SiC	1.19	1.69	1.89	2.18	2.14	2.23	1.67	1.87	2.29c
GaAs	0.30	0.47	1.05	1.43	0	0	0.76	1.25	1.72c
InP	0.41	0.65	1.24	1.76	0	0.14	0.82	1.37
CdS	0.74	1.15	1.86	2.86	0.93	1.39	1.36	1.75
LiH	2.52	3.52	4.26	6.16	2.58	2.58	3.50	3.51
LiCl	5.88	6.63	7.60	9.80	6.66	6.70	6.77	6.84
BeO	7.41	8.05	8.62	9.40	8.97	8.97	8.31	8.39
MgO	4.57	5.09	5.54	6.68	5.50	5.49	5.46	5.49	6.31d
CsF	4.75	5.44	6.08	8.15	7.30	7.68	5.97	6.33
BaO	1.69	2.12	2.52	3.49	3.08	3.42	2.31	2.63
PbS	0.11	0.52	0.98	1.60	0.93	1.13	0.47	0.49
ScN	0	0.10	0.35	0.70	0	0	0.09	0.23	1.58e
SrTiO3	1.77	1.95	2.14	2.33	1.40	1.49	2.01	2.15	4.20e
MnO	0.73	1.30	1.74	2.58	2.13	2.51	1.74	1.81	3.85f
FeO	0	0	0.25	0.84	0	0	0.32	0.44	1.66f
CoO	0	0.01	0.50	1.37	0	0.24	0.69	1.20	2.62f
NiO	0.52	1.15	1.57	2.07	0.58	1.28	1.76	2.50	4.10f, 3.54d
ZnO	0.53	0.94	1.38	2.06	2.58	2.77	1.41	1.53
Cu2O	0.50	0.66	0.72	0.85	1.40	1.40	0.73	0.73	1.44d
CeO2	2.03	2.05	2.11	2.16	0.95	0.94	2.02	2.13
ME	–0.88	–0.48	–0.06	1.39	–0.04	0.07	–0.18
MAE	0.88	0.48	0.32	1.42	0.77	0.73	0.21

The EXX-OEP results from previous works are also shown. The ME and MAE are with respect to BJS.

From ref (76).

From ref (77).

From ref (36).

From ref (78).

From ref (13) (LDA correlation potential[79] was added to EXX-OEP).

Table 5

Total Energies (in Ry/Cell) Obtained by Evaluating the HF Total-Energy Expression (i.e., eq for Ex and No Correlation) with Orbitals Generated from Various Exchange-Only Potentialsa

solid	LDA	PBE	EV93	AK13	BR(0.8)	S	BJBR(0.8)	EXX-OEPb	BJS
Ne	0.018	0.007	0.002	–0.004	0.123	0.069	–0.012		–257.359
Ar	0.016	0.021	0.009	0.043	0.265	0.225	0.007		–1057.361
Kr	0.058	0.024	0.000	0.037	0.656	0.578	0.004		–5577.773
Xe	0.065	0.023	–0.001	0.044	1.000	0.924	0.006		–14894.471
C	0.038	0.022	0.012	0.025	0.182	0.115	0.004	–0.004	–151.588
Si	0.064	0.026	0.001	0.014	0.457	0.346	0.004	–0.014	–1158.334
Ge	0.156	0.097	0.031	0.033	1.263	1.091	0.026		–8390.306
Se	0.207	0.108	0.020	0.053	1.959	1.705	0.030		–14571.924
BN	0.039	0.019	0.007	0.015	0.170	0.106	0.000	–0.008	–158.615
SiC	0.052	0.022	0.001	0.006	0.309	0.221	0.003		–655.021
GaAs	0.156	0.093	0.031	0.032	1.249	1.077	0.024		–8404.414
InP	0.121	0.058	0.009	0.027	1.180	1.050	0.022		–12444.585
CdS	0.129	0.056	0.001	0.017	1.120	0.998	0.016		–11984.420
LiH	0.019	0.005	0.006	0.034	0.024	0.022	0.000		–16.124
LiCl	0.046	0.013	0.001	0.045	0.262	0.215	0.001		–937.133
BeO	0.083	0.024	–0.009	0.000	0.292	0.172	–0.012		–359.099
MgO	0.048	0.006	–0.020	–0.015	0.264	0.163	–0.015	–0.031	–550.099
CsF	0.094	0.028	–0.010	0.024	1.111	0.977	–0.003		–15773.653
BaO	0.101	0.035	–0.001	0.019	1.159	1.014	0.005		–16422.281
PbS	0.163	0.091	0.042	0.063	2.058	1.826	0.041		–42637.183
ScN	0.062	0.028	0.003	0.018	0.556	0.448	0.016		–1636.115
SrTiO3	0.187	0.073	0.003	0.048	1.546	1.217	0.019		–8512.403
MnO	0.236	0.093	–0.024	–0.058	1.127	0.836	0.025		–4930.288
FeO	0.717	0.496	0.114	–0.099	1.844	1.431	0.115		–5386.094
CoO	1.121	0.372	0.119	–0.009	2.284	1.274	0.100		–5868.435
NiO	0.628	0.321	0.132	0.030	1.891	1.319	0.110	–0.298	–6377.441
ZnO	0.255	0.123	0.016	–0.012	1.246	0.962	–0.002		–7478.489
Cu2O	0.383	0.156	0.009	–0.045	2.062	1.612	–0.058	–0.290	–13527.450
CeO2	0.120	0.053	0.015	0.052	1.956	1.711	0.029		–18023.332
MRE	10.3	3.9	1.9	8.5	35.0	25.6	0.1
MARE	10.3	3.9	2.4	9.1	35.0	25.6	0.8

The results from ref (36) obtained with the EXX-OEP orbitals are also shown. For ease of comparison, we show the difference with respect to the values obtained with the BJS orbitals shown in the last column. A negative value indicates a more negative total energy than with BJS orbitals. The MRE and MARE (in pcm) are with respect to BJS.

From ref (36).

Electronic Structure

Starting with the comparison of the fundamental band gaps (Table ) obtained with the Slater and BR(0.8) potentials alone (i.e., not augmented with the additional BJ term), we can see that these two potentials lead to relatively similar values in many cases. Actually, the agreement can be considered to be excellent for the rare-gas solids, C, Se, BN, SiC, LiH, LiCl, BeO, MgO, SrTiO3, Cu2O, and CeO2. For these systems, the disagreement is below 0.1 eV (slightly more for Ne, which has a very large band gap). The largest discrepancies between Slater and BR(0.8) are found for CsF, MnO, NiO, and CdS (the Slater band gaps are larger by at least 0.4 eV), and for Si, the difference of 0.3 eV is also relatively important since it represents ∼30% of the band gap. InP and CoO are described as (semi)metal with BR(0.8), whereas the Slater potential leads to a nonzero band gap for these systems. The results obtained with the Slater and BR(0.8) potentials when they are combined with the second term in eq show the same trends. However, in some cases, like C, Si, or BN, slightly larger differences between BJS and BJBR(0.8) can be noticed. Also, Ge, GaAs, InP, and CoO are described as semiconductors by BJS and BJBR(0.8), and the discrepancies between the two methods are above 0.4 eV. By considering the BJS band gaps as reference, the best agreement is obtained with the BJBR(0.8) potential, which leads to the lowest mean absolute error (MAE) (Table ) among all tested potentials. Note, however, that the EV93, BR(0.8), and Slater potentials lead to small mean error (ME) of ∼0.05 eV, indicating that these potentials do not show a particular trend to under- or overestimate the band gap with respect to BJS.

From the results in Table , it is also interesting to note that when the second term in eq is added to the Slater or BR(0.8) potential the band gap is increased in some cases but decreased in others. GaAs, InP, LiH, CoO, NiO, and CeO2 are examples that show an increase of roughly 1 eV, whereas for Ne, CsF, BaO, MnO, Cu2O, and ZnO, there is a decrease that can also reach 1 eV. On the other hand, for Si, LiCl, and MgO, there is very little change in the band gap.

For the comparison with the EXX-OEP results (also shown in Table when available), it has already been shown[17,77,80,81] that the LDA and PBE band gaps are much smaller. In passing, we note that this underestimation by LDA/PBE is strongly reduced if correlation from the random-phase approximation (RPA-OEP) is added to the EXX-OEP[17,82] because, usually, total (exchange and correlation) semilocal approximations are more accurate than exchange or correlation alone. The EV93 and AK13 potentials systematically increase the band gap with respect to LDA and PBE, and in all cases, the magnitude of the band gap follows the order LDA < PBE < EV93 < AK13. However, from the comparison of the Slater/BR-based potentials with EV93 and AK13, no systematic trend can be observed, and their relative performances with respect to EXX-OEP depend on the solid. For instance, for the rare-gas solids, the Slater/BR-based potentials are the most accurate, whereas AK13 leads to strong overestimations of the band gap by several electron volts. For the transition-metal oxides, one or another of the Slater/BR-based methods is closer to EXX-OEP, whereas for Ge, GaAs, and MgO, AK13 leads to the best agreement with EXX-OEP. Thus, overall, there is no method that reproduces the EXX-OEP band gaps systematically better than the others. It is noteworthy that for Cu2O the agreement between the Slater/BR(0.8) and EXX-OEP potentials is excellent, whereas in ref (36), a qualitative agreement with EXX-OEP could be obtained only with the generalized BJ (gBJ) potential including the universal correction (UC).[44] In the present work, it is shown that adding the BJ term to the Slater or BR(0.8) potential destroys the agreement with EXX-OEP for Cu2O. For Si, Ge, GaAs, MnO, CoO, and NiO, it is clear that the use of the Slater/BJS potential instead of BR(0.8)/BJBR(0.8) leads to a much better agreement with EXX-OEP.

In order to provide some insight into the results discussed so far for the fundamental band gap, we show in Figures and 2 one-dimensional plots of the exchange potentials in Kr and BaO, respectively. The electron densities of the VBM and CBM are also shown. The trends in the band gap can be explained by the following features in the potentials, which are rather similar in various solids. With some exceptions, like transition-metal oxides with a d–d band gap (see, e.g., ref (47) for a detailed discussion on Cu2O), the VBM and CBM extend in different regions of space, and, typically, the VBM is mainly located around the atoms and in the bonding regions, whereas an important part of the CBM is in the interstitial region, where the VBM has essentially no contribution (see lower panels of Figures and 2). Therefore, in such situations, the band gap should depend strongly on the change in the magnitude of the potential when going from the valence region around the atoms to the interstitial. In this respect, we can see that LDA shows the smallest change and therefore the smallest band gaps. The GGA potentials are more positive than LDA in the interstitial region, which is the main source for the increase of the band gap with respect to LDA. The same mechanism explains the band gaps obtained with the BJS and BJBR(0.8) potentials, since they are also more repulsive than LDA and PBE in a rather large area of the interstitial as well as slightly more attractive closer to the nuclei. In the case of BaO, for instance (Figure ), we can see that in the interstitial region the BJS potential is more repulsive than BJBR(0.8), thus leading to an upward shift of the CBM and a band gap larger by 0.3 eV. In ref (36), it was noted that in the interstitial region, the EXX-OEP and BJ-based potentials are very similar and smooth, not showing the large peaks observed in the EV93 and AK13 potentials (and, to a lesser extent, also with PBE). Concerning the Slater and BR(0.8) potentials, also shown in Figures and 2, their shape looks rather similar to LDA since the intershell peaks are also absent. However, these two potentials are more negative by a factor of ∼3/2, which has a stronger effect in the high-density region close to the nuclei where the VBM is located, thus leading to a band gap that is enlarged with respect to LDA.

Figure 1

Upper panel: Exchange potentials vx in Kr plotted from the atom at (0,0,0) (d = 0) to the mid-distance to the atom at (1,0,0) (d = 5.598 Å). The maximum positive values of vx for PBE, EV93, and AK13 are 0.7, 1.5, and 15.7 Ry, respectively. Lower panel: Electron density of the VBM and CBM (normalized to one electron) plotted along the same path as that in the upper panel.

Figure 2

Upper panel: Exchange potentials vx in BaO plotted from the Ba atom at (0,0,0) (d = 0) to the O atom at (1/2,1/2,1/2) (d = 4.783 Å). The maximum positive values of vx for PBE, EV93, and AK13 are 0.2, 0.6, and 2.9 Ry, respectively. Lower panel: Electron density of the VBM and CBM (normalized to one electron) plotted along the same path as that in the upper panel.

Previously,[36] it was shown that only the gBJ potential with the universal correction (gBJUC) could lead to qualitative agreement with EXX-OEP for the fundamental band gap and EFG in Cu2O. This was explained by the similarities of the EXX-OEP and gBJUC potentials in the relevant regions of space (see Figure 8b of ref (36)). As noticed above, the results obtained with the Slater and BR(0.8) potentials agree closely with EXX-OEP for the band gap and reasonably well for the EFG (see below). Figure shows that, indeed, the Slater and BR(0.8) potentials (without the BJ response term) agree very well with EXX-OEP near the Cu atom.

Figure 3

Exchange potentials vx in Cu2O plotted starting at a distance of 1 Å from the Cu atom at site (1/2,1/2,0) (d = 0) in the direction of the O atom at site (3/4,3/4,3/4) (d = 3.538 Å). The EXX-OEP from ref (36) is also shown.

Figure shows the difference vx,↑ – vx,↓ between the spin-up and spin-down exchange potentials in NiO. The EXX-OEP taken from ref (36) is also included for comparison. The occupation of the Ni-3d states in NiO (spin-up full and t2g↓ occupied for the Ni atom at the left upper corner) is such that the shape of vx,↑ – vx,↓ reflects the empty eg↓ orbitals. Depending on the method, the band gap in NiO can be of d–d character, charge-transfer character, or a mixture of both. For instance (see, e.g., refs (36) and (83)), LDA/PBE, EXX-OEP, and the BJ-based potentials lead to a d–d band gap, whereas with LDA+U and the onsite-hybrid functionals, a mixed d–d/charge-transfer band gap is obtained. The HF method leads to a band gap of pure charge-transfer character.[36,84] As discussed in ref (36), the more the angular eg character is pronounced, the larger the d–d band gap should be. Indeed, from Figure , we can see that the shape of vx,↑ – vx,↓ on the Ni atom correlates rather well with the results for the band gap in Table . For instance, the potentials LDA, BJBR(0.8), BJS, and EXX-OEP (in this order) lead to band gaps that are larger and larger by step of ∼1 eV (similarly for the magnetic moment, see below), whereas the features of an eg orbital become more and more pronounced. Note also that the LDA and BR(0.8) potentials look rather similar and their band gaps differ by only 0.06 eV. In ref (36), an excellent agreement with the EXX-OEP for NiO could be obtained with the gBJ potential and with optimized parameters.

Figure 4

Two-dimensional plots of the difference between spin-up and spin-down exchange potentials (vx,↑ – vx,↓) in a (001) plane of antiferromagnetic NiO. The contour lines start at −2 Ry (blue color) and end at 2 Ry (red color), with an interval of 0.235 Ry. The Ni atom with a full spin-up 3d shell is at the left upper corner. The plot for EXX-OEP is taken from ref (36).

Concerning the core states, indicated in Table (for LiH, the Li-1s state was considered for the present analysis), we show in the middle panel of Figures S1–S29 the mean absolute relative error (MARE) (BJS is the reference) on the energy position of the cores states with respect to the Fermi energy (set at the VBM for solids with a nonzero band gap). In the vast majority of cases, BJBR(0.8) (considering, for the moment, only the results for γ = 0.8) leads to the lowest MARE (below ∼0.3–0.4%). The MARE for the Slater and BR(0.8) potentials (not shown) is larger than that for all other potentials by 1 order of magnitude (in the range 2–6%), which is a consequence of the too negative potential in the neighboring region of the nuclei, leading to orbitals that are too localized around the nuclei (i.e., too deep in energy). This is the same problem that was encountered with the gBJUC potential, as reported in ref (36).

Magnetic Moment and EFG

Turning now to magnetic properties, Table shows the spin magnetic moment μS in the antiferromagnetic transition-metal monoxides, and from the results, we can see that the trends among them are different. For instance, although BR(0.8) leads to one of the largest magnetic moment in MnO, it leads to the smallest value in CoO and NiO. This observation should be due mainly to the BJ response term that enhances magnetism more in CoO and NiO, which have a nonspherical 3d shell, than in MnO, which has a spherical 3d shell. Also, the agreement between the BR(0.8) and Slater potentials is rather good for MnO (the differences are of a few 0.01μB), but it is not for FeO and NiO, for which the discrepancies reach ∼0.2μB when the BJ response term is not added. As is well-known,[13,75,83,85−87] the LDA and standard GGAs like PBE strongly underestimate μS with respect to both experiment and EXX-OEP, which is a general problem of systems with localized 3d electrons. The EV93, AK13, and BJS/BJBR(0.8) potentials improve the results; however, the values are still clearly underestimated with respect to EXX-OEP as for the band gap. In ref (36), a magnetic moment of 1.86μB for NiO was obtained with the gBJ potential.

Table 3

Spin Magnetic Moment μS (in μB) of the Transition-Metal Atom in Antiferromagnetic MnO, FeO, CoO, and NiO Calculated from Different Exchange-Only Potentials

solid	LDA	PBE	EV93	AK13	BR(0.8)	S	BJBR(0.8)	BJS	EXX-OEP
MnO	4.18	4.23	4.30	4.39	4.39	4.43	4.25	4.28	4.81a
FeO	3.41	3.44	3.48	3.51	3.44	3.62	3.46	3.54	3.85a
CoO	2.44	2.46	2.53	2.59	2.36	2.50	2.53	2.61	2.88a
NiO	1.30	1.43	1.51	1.58	1.17	1.39	1.53	1.63	1.89a,1.91b

From ref (13) (LDA correlation potential[79] was added to EXX-OEP).

From ref (36).

The results for the EFG in Se and at the Cu site in Cu2O are shown in Table . For both systems, the Slater and BR(0.8) potentials lead to very similar values. However, for Cu2O, a (moderate) discrepancy of 1.5 × 1021 V/m2 between BJS and BJBR(0.8) is obtained. Similarly, as for the band gap in Cu2O discussed above, it was also shown in ref (36) that a good agreement with the EXX-OEP value (−17.7 × 1021 V/m2) could be obtained only with the gBJ potential including the universal correction (−15 × 1021 V/m2), whereas the magnitude of the EFG obtained with all other potentials tested in ref (36) was below 10 × 1021 V/m2. However, in the present work, a magnitude well above 10 × 1021 V/m2 is obtained with the Slater and BR(0.8) potentials, which is in line with the agreement with EXX-OEP for the band gap.

Table 4

EFG (in 1021 V/m2) of Se and Cu in Cu2O Calculated from Different Exchange-Only Potentials

solid	LDA	PBE	EV93	AK13	BR(0.8)	S	BJBR(0.8)	BJS	EXX-OEP
Se	–49.1	–53.2	–55.5	–60.7	–68.0	–67.0	–51.4	–50.8
Cu2O	–4.9	–5.7	–6.8	–8.0	–13.2	–13.3	–7.2	–5.7	–17.7a

From ref (36).

HF Total Energy

Next, we consider the total energy as a measure of the difference between the orbitals generated by the various exchange potentials, and Table shows the HF total energies calculated with the various sets of orbitals. The results show that, compared to the values obtained with the BJS orbitals, the EV93 and BJBR(0.8) orbitals lead to differences that are the smallest. The mean relative error (MRE) and mean absolute relative error (MARE) are the smallest for BJBR(0.8), which is followed rather closely by EV93. In some cases, like the rare-gas solids, AK13 leads to somewhat larger differences, but it improves overall upon LDA and PBE. As already mentioned above, the Slater and BR(0.8) potentials (without the BJ response term) localize the core orbitals too much; therefore, the large contribution to the total energy coming from the high-density region inside the atoms is very inaccurate, leading to total energies that differ from all other methods by 0.2–2 Ry/cell.

We recall that the EXX-OEP is, within the space of multiplicative potentials, the one providing the orbitals that minimize the HF total energy.[14,15] In our previous work,[36] we proposed a set of parameters for the gBJ potential such that the deviations from EXX-OEP are in the range 0.001–0.003 Ry/cell (for C, Si, BN, and MgO), representing an improvement with respect to the original BJBR(0.8), for which the disagreement with EXX-OEP was more on the order of 0.01–0.02 Ry/cell. In the present work, we can see that the BJS orbitals lead, in most cases (six exceptions), to total energies that are lower than those from BJBR(0.8). Similarly, the EV93 and AK13 orbitals give lower total energies than BJS orbitals in six or seven cases. Overall, in terms of total energy, the BJS potential seems to be the closest to the EXX-OEP since it leads to the lowest value for more than half of the solids. In particular, it is important to note that the replacement of the Slater potential in eq by BR(0.8) leads, in general, to a degradation of the results when compared with the EXX-OEP. This conclusion is in line with the results from Becke and Johnson,[22] who showed that the BJS orbitals usually lead to slightly lower HF total energies for isolated atoms. One of the exceptions was the Ne atom, which is also in agreement with our results for solid Ne.

Influence of the Parameter γ

In their work comparing the Slater and BR potentials for atoms and molecules, Heßelmann and Manby[23] optimized for each system the parameter γ in eq . More specifically, they adjusted the value of γ such that the exchange energy (calculated according to the last line of eq ) yields the exact HF value. As mentioned above, the gBJ potential proposed in ref (36) as an approximation to the EXX-OEP contains three parameters (see ref (36) for details), one of which is γ that was varied between 0.4 and 1.4. In a similar way, the effect of γ in BJBR has been investigated in the present work in order to see the extent to which the agreement with BJS could be improved with respect to the standard value γ = 0.8 used in this work so far. γ has been varied in steps of 0.1 between 0.4 and 1.4. With the BJBR(γ) potential, the results for the fundamental band gap, MARE on the core states, and HF total energy are displayed in Figures S1–S29 of the Supporting Information (because, in most cases, the results from BR(γ) alone are very bad, they are not shown).

A concise summary of the observed trends is the following. For most solids, an increase of γ in BJBR(γ) leads to an increase of the band gap, with the exceptions of Ne, BaO, MnO, and ZnO. In some cases, e.g., the rare gases, CsF or PbS, the parameter γ has little effect on the band gap. For about 10 of the solids studied, a value of γ larger than 1.2 leads to better agreement with BJS, whereas for most other solids, a value around γ ∼ 0.8 leads to a quite reasonable agreement with BJS. For C, Si, Ge, InP, and NiO, a value of γ above 1.4 would be required for a perfect agreement for the band gap with BJS. Concerning the total energy (right panel of Figures S1–S29), it is rather satisfying to see that for many of the solids (about 20 of them) the optimal value for γ lies in the range 0.7–0.9, enclosing the original value 0.8. Furthermore, for most of these solids, a value for γ in this range also leads to the most negative total energy. This finding is at variance with the results from Heßelmann and Manby,[23] since, for most molecules, they considered the optimized value of γ to be in the range 1.1–1.3. However, note that the molecules in their test set contain only light atoms; furthermore, they used a different procedure for optimizing γ, as mentioned above. FeO, CoO, and NiO are the worst cases, which would require a value larger than 1.4 in order to reach both the best agreement with BJS and the lowest total energy. The results for the MARE on the core states (middle panel of Figures S1–S29) show that for about half of the cases there is a very good correlation with the total energy such that the agreement with BJS is reached for both quantities with very similar values of γ. This seems to be particularly the case for solids containing heavy atoms, whose contribution to the total energy comes mainly from their core states. As a conclusion, a value of ∼0.8 for γ in the BJBR(γ) potential seems to be the most reasonable choice, overall.

Regarding the influence of the parameter γ in BJBR(γ) on the magnetic moment μS in MnO, FeO, CoO, and NiO, the results are shown in Figures S30–S33. For MnO and FeO, μS gets smaller when γ is increased, whereas the reverse occurs for NiO. For CoO, the smallest value of μS is obtained with γ = 0.9. For MnO and FeO, a better agreement with BJS is obtained for a small γ, whereas NiO requires a value for γ larger than 1.4. Figures S34 and S35 show the results for the EFG in Se and Cu2O, respectively, where we can see that for Se a perfect match with the BJS potential is obtained for γ = 0.9, whereas for Cu2O, the discrepancy with respect to BJS is always larger than 1 × 1021 V/m2.

Summary and Conclusions

The BR potential leads to calculations that are 1 or 2 orders of magnitude faster than those with the Slater potential; therefore, it is tempting to use the former instead of the latter, especially for applications on large systems. However, until now, no such investigation on the differences between these two potentials in solids has been done, and the present work attempts to fill the gap. For this purpose, we have compared the results obtained with the Slater and BR potentials for the electronic structure, EFG, magnetic moment, and total energy. The test set consists of semiconductors and insulators of various types.

The results indicate that, in many cases, the BR potential is a good approximation to the Slater potential. In particular, for the purpose of comparing the BJ potential with the EXX-OEP, it does not really matter which version of BJ (i.e., BJS or BJBR) is used. However, this is not systematically the case. For instance, in the strongly correlated systems, FeO, CoO, and NiO, the band gaps and/or the magnetic moments can differ significantly. In addition, rather large differences in the band gap were also observed for Si, Ge, and systems containing heavy atoms like CsF, BaO, and CdS. Interestingly, in such cases, the agreement with EXX-OEP is better when BJS is used. The influence of the parameter γ (in BJBR) on the results has also been investigated, and the conclusion is that γ = 0.8 is a rather good (but not universal) choice, as no other value of γ seems to lead to better results on average.