Fully numerical calculations on atoms with fractional occupations. Range-separated exchange functionals

TL;DR

This paper extends a finite element approach for fully numerical atomic calculations to include fractional occupations, enabling rapid, highly accurate computations of atomic energies and potentials, including for heavy atoms and range-separated functionals.

Contribution

It introduces a numerical method for atoms with fractional occupations, incorporating range-separated exchange functionals and efficient potential calculations, advancing atomic structure modeling.

Findings

Achieves microhartree accuracy with few basis functions.

Reproduces literature results for atoms with Z=1 to 86.

Validates erfc kernel implementation against Gaussian basis set results.

Abstract

A recently developed finite element approach for fully numerical atomic structure calculations [S. Lehtola, Int. J. Quantum Chem. 119, e25945 (2019)] is extended to the description of atoms with spherically symmetric densities via fractionally occupied orbitals. Specialized versions of Hartree-Fock as well as local density and generalized gradient approximation density functionals are developed, allowing extremely rapid calculations at the basis set limit on the ground and low-lying excited states even for heavy atoms. The implementation of range-separation based on the Yukawa or complementary error function (erfc) kernels is also described, allowing complete basis set benchmarks of modern range-separated hybrid functionals with either integer or fractional occupation numbers. Finally, computation of atomic effective potentials at the local density or generalized gradient approximation levels for the superposition of atomic potentials (SAP) approach [S. Lehtola, J. Chem. Theory Comput. 15, 1593 (2019)] that has been shown to be a simple and efficient way to initialize electronic structure calculations is described. The present numerical approach is shown to afford beyond microhartree accuracy with a small number of numerical basis functions, and to reproduce literature results for the ground states of atoms and their cations for $1 \leq Z \leq 86 $. Our results indicate that the literature values deviate by up to 10 {\mu}Eh from the complete basis set limit. The numerical scheme for the erfc kernel is shown to work by comparison to results from large Gaussian basis set calculations from the literature. Spin-restricted ground states are reported for Hartree-Fock and Hartree-Fock-Slater calculations with fractional occupations for $1 \leq Z \leq 118$.

Tables

Table 1. Table 1: Total energies of neutral atoms for calculations with the VWN functional and a PBE functional based on VWN correlation. The differences Δ E = E ( present work) − E ( Kraisler 2010 ) \Delta E=E(\text{present work)}-E(\text{Kraisler 2010}) are calculated relative to the fully numerical values from ref. 82 and are reported in microhartree.

$Z$	$E\left(\text{VWN}\right)$		$\Delta E\left(\text{VWN}\right)$	$E\left(\text{PBE}\right)$		$\Delta E\left(\text{PBE}\right)$	$Z$	$E\left(\text{VWN}\right)$		$\Delta E\left(\text{VWN}\right)$	$E\left(\text{PBE}\right)$		$\Delta E\left(\text{PBE}\right)$
1	-0.	478671	0	-0.	499963	-0	44	-4439.	044607	1	-4443.	255631	-6
2	-2.	834836	0	-2.	893287	-1	45	-4683.	360538	1	-4687.	685035	-6
3	-7.	343957	0	-7.	462726	-0	46	-4935.	368406	0	-4939.	811859	-5
4	-14.	447209	-0	-14.	630525	-0	47	-5195.	037351	1	-5199.	600560	-6
5	-24.	353614	0	-24.	606283	-0	48	-5462.	390982	1	-5467.	070748	-5
6	-37.	470031	0	-37.	795116	-0	49	-5737.	313809	-0	-5742.	114862	-6
7	-54.	136799	0	-54.	537743	-0	50	-6019.	972345	0	-6024.	895821	-6
8	-74.	527410	0	-75.	003219	-1	51	-6310.	419326	1	-6315.	467009	-6
9	-99.	114192	0	-99.	668342	-1	52	-6608.	650476	0	-6613.	811930	-6
10	-128.	233481	-0	-128.	869661	-1	53	-6914.	777857	1	-6920.	056261	-6
11	-161.	447625	-0	-162.	176267	-1	54	-7228.	856106	1	-7234.	254478	-7
12	-199.	139406	-0	-199.	958820	-1	55	-7550.	561866	0	-7556.	083016	-7
13	-241.	321156	-0	-242.	236076	-1	56	-7880.	111578	-0	-7885.	752916	-7
14	-288.	222945	0	-289.	236535	-1	57	-8217.	648931	91	-8223.	406822	220
15	-340.	005794	-0	-341.	120757	-2	58	-8563.	489711	880	-8569.	364450	663
16	-396.	743948	0	-397.	951200	-1	59	-8917.	715777	6966	-8923.	707847	6539
17	-458.	671463	-0	-459.	976078	-2	60	-9280.	405670	16518	-9286.	515623	16052
18	-525.	946195	0	-527.	352025	-2	61	-9651.	650420	13676	-9657.	878495	13601
19	-598.	206032	-0	-599.	716752	-2	62	-10031.	516930	6437	-10037.	864135	6383
20	-675.	742283	0	-677.	355243	-2	63	-10420.	023146	1	-10426.	490411	-8
21	-758.	685248	0	-760.	397795	-3	64	-10817.	148260	858	-10823.	727509	285
22	-847.	314902	-0	-849.	129808	-3	65	-11223.	108037	4860	-11229.	800083	3513
23	-941.	786662	-0	-943.	704413	-2	66	-11637.	977781	11030	-11644.	783485	9215
24	-1042.	218348	-0	-1044.	239902	-3	67	-12061.	832318	18521	-12068.	752563	16445
25	-1148.	644093	0	-1150.	765417	-3	68	-12494.	746152	26769	-12501.	781831	24564
26	-1261.	223291	6017	-1263.	441835	4291	69	-12936.	809752	19185	-12943.	957962	20751
27	-1380.	193787	716	-1382.	508399	831	70	-13388.	048594	1	-13395.	317842	-9
28	-1505.	672905	1	-1508.	087914	-4	71	-13848.	234767	1	-13855.	623680	-9
29	-1637.	793358	0	-1640.	310279	-4	72	-14317.	517965	2	-14325.	032671	-9
30	-1776.	573850	0	-1779.	194575	-4	73	-14795.	971453	1	-14803.	612704	-9
31	-1921.	851924	-0	-1924.	582672	-3	74	-15283.	610347	2	-15291.	380462	-9
32	-2073.	829860	0	-2076.	672928	-4	75	-15780.	381133	2	-15788.	268506	-8
33	-2232.	587154	0	-2235.	545023	-4	76	-16286.	434007	2	-16294.	440422	-9
34	-2398.	134930	0	-2401.	196896	-4	77	-16801.	850893	2	-16809.	976281	-8
35	-2570.	626651	-0	-2573.	796934	-4	78	-17326.	660985	3	-17334.	912620	-10
36	-2750.	147940	1	-2753.	430126	-4	79	-17860.	796573	2	-17869.	175326	-10
37	-2936.	342160	0	-2939.	739646	-5	80	-18404.	274220	1	-18412.	777007	-10
38	-3129.	453161	1	-3132.	963153	-5	81	-18956.	962102	1	-18965.	593468	-10
39	-3329.	525142	0	-3333.	148098	-5	82	-19519.	010773	2	-19527.	771776	-10
40	-3536.	771074	-0	-3540.	515940	-5	83	-20090.	453943	1	-20099.	346370	-9
41	-3751.	295618	0	-3755.	160742	-6	84	-20671.	273855	2	-20680.	287630	-9
42	-3973.	162595	-0	-3977.	149787	-5	85	-21261.	559507	2	-21270.	697436	-10
43	-4202.	348934	1	-4206.	446961	-5	86	-21861.	346869	3	-21870.	611766	-9

Table 2. Table 2: Total energies of atomic cations for calculations with the VWN functional and a PBE functional based on VWN correlation. The differences Δ E = E ( present work) − E ( Kraisler 2010 ) \Delta E=E(\text{present work)}-E(\text{Kraisler 2010}) are calculated relative to the fully numerical values from ref. 82 and are reported in microhartree. a An incorrect value was reported in ref. 82 for \ce La+; see main text.

$Z$	$E\left(\text{VWN}\right)$		$\Delta E\left(\text{VWN}\right)$	$E\left(\text{PBE}\right)$		$\Delta E\left(\text{PBE}\right)$	$Z$	$E\left(\text{VWN}\right)$		$\Delta E\left(\text{VWN}\right)$	$E\left(\text{PBE}\right)$		$\Delta E\left(\text{PBE}\right)$
1	0.	000000	0	0.	000000	0	44	-4438.	767375	1	-4442.	987586	-6
2	-1.	941703	0	-1.	993741	-1	45	-4683.	055688	-0	-4687.	391656	-6
3	-7.	142818	-0	-7.	257274	-0	46	-4935.	023835	-0	-4939.	477807	-6
4	-14.	115512	-0	-14.	299957	-0	47	-5194.	755725	0	-5199.	329663	-6
5	-24.	038275	-0	-24.	294185	-0	48	-5462.	065582	1	-5466.	758188	-6
6	-37.	037413	0	-37.	365761	-1	49	-5737.	101461	1	-5741.	909394	-6
7	-53.	585407	-0	-53.	989486	-1	50	-6019.	697098	0	-6024.	624904	-6
8	-74.	016721	-0	-74.	500466	-1	51	-6310.	085269	0	-6315.	135445	-6
9	-98.	450427	0	-99.	012379	-1	52	-6608.	317272	0	-6613.	492318	-6
10	-127.	418114	-0	-128.	061506	-1	53	-6914.	378893	1	-6919.	666895	-6
11	-161.	250340	0	-161.	979238	-2	54	-7228.	394173	0	-7233.	799075	-7
12	-198.	855669	0	-199.	679196	-1	55	-7550.	416737	0	-7555.	942187	-7
13	-241.	100595	0	-242.	016828	-2	56	-7879.	920522	984	-7885.	569221	813
14	-287.	918773	-0	-288.	932217	-1	57a	-8217.	455429	1545	-8223.	221428	2463
15	-339.	618451	-0	-340.	732787	-2	58	-8563.	294377	7617	-8569.	175665	6993
16	-396.	356540	0	-397.	574834	-1	59	-8917.	518819	20242	-8923.	517325	19853
17	-458.	184562	0	-459.	496045	-2	60	-9280.	230184	12512	-9286.	344462	14304
18	-525.	360439	-0	-526.	770845	-2	61	-9651.	482296	5894	-9657.	715083	7570
19	-598.	039506	-0	-599.	553203	-2	62	-10031.	348348	2185	-10037.	700262	3563
20	-675.	514035	-0	-677.	132344	-3	63	-10419.	820399	1	-10426.	294380	-8
21	-758.	442642	1812	-760.	161910	1575	64	-10816.	942313	782	-10823.	528176	-8
22	-847.	065015	2847	-848.	886280	2949	65	-11222.	899249	11493	-11229.	597850	8739
23	-941.	523838	0	-943.	448454	-2	66	-11637.	769436	20564	-11644.	578605	20710
24	-1041.	944126	0	-1043.	972993	-3	67	-12061.	641951	12969	-12068.	562285	16381
25	-1148.	368924	-0	-1150.	502293	-3	68	-12494.	571851	7644	-12501.	608100	10480
26	-1260.	927746	-0	-1263.	157817	-3	69	-12936.	633640	4018	-12943.	786682	5267
27	-1379.	896444	0	-1382.	218969	-3	70	-13387.	827982	1	-13395.	103732	-10
28	-1505.	370040	0	-1507.	793428	-4	71	-13847.	999554	2	-13855.	396289	-9
29	-1637.	485140	0	-1640.	010907	-4	72	-14317.	267886	1533	-14324.	788056	2091
30	-1776.	217890	0	-1778.	850041	-3	73	-14795.	705335	1	-14803.	354716	-8
31	-1921.	629140	-0	-1924.	365546	-4	74	-15283.	334783	2	-15291.	113757	-9
32	-2073.	533337	0	-2076.	379628	-4	75	-15780.	100605	2	-15788.	005290	-8
33	-2232.	220332	0	-2235.	179888	-4	76	-16286.	150952	2	-16294.	166211	-8
34	-2397.	770127	-0	-2400.	845945	-4	77	-16801.	535551	3	-16809.	673022	-8
35	-2570.	180737	-0	-2573.	360358	-4	78	-17326.	305178	1	-17334.	567568	-10
36	-2749.	623528	1	-2752.	911995	-4	79	-17860.	511437	0	-17868.	901185	-10
37	-2936.	183045	0	-2939.	584557	-5	80	-18403.	949940	1	-18412.	465847	-10
38	-3129.	240801	-0	-3132.	756842	-5	81	-18956.	753577	1	-18965.	392384	-10
39	-3329.	295616	-0	-3332.	926478	-5	82	-19518.	743995	2	-19527.	509994	-10
40	-3536.	524561	871	-3540.	277393	131	83	-20090.	133387	2	-20099.	028863	-9
41	-3751.	036938	0	-3754.	910732	-5	84	-20670.	954321	2	-20679.	981706	-9
42	-3972.	894262	-0	-3976.	890719	-5	85	-21261.	180969	2	-21270.	328813	-8
43	-4202.	074923	0	-4206.	181883	-5	86	-21860.	912366	3	-21870.	184141	-9

Table 3. Table 3: Comparison of the total energies of spherically symmetric atoms with the LC-BLYP functional with the range separation constant ω = 0.3 𝜔 0.3 \omega=0.3 reproduced with finite element calculations (present work) and a Gaussian basis set calculation (Anderson  et al. , 2017 ) .

atom	finite element		Gaussian basis
\ceH-	-0.	519949	-0.	51995
He	-2.	866811	-2.	86681
\ceLi-	-7.	435511	-7.	43551
Be	-14.	584723	-14.	58472
N	-54.	482223	-54.	48222
\ceF-	-99.	766050	-99.	76604
Ne	-128.	816627	-128.	81661
\ceNa-	-162.	136564	-162.	13655
Mg	-199.	907036	-199.	90702
P	-341.	069932	-341.	06992
\ceCl-	-460.	080588	-460.	08057
Ar	-527.	321257	-527.	32124

Table 4. Table 4: Non-relativistic spin-restricted spherical HF configurations for all elements in the periodic table. Entries in italic indicate a lower-lying configuration was identified but it failed to converge.

H	$1s^1$	-0.357710	Nb	$\text{[Kr]}5s^{2}4d^3$	-3753.033166	Tl	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^1$	-18961.740923
He	$1s^2$	-2.861680	Mo	$\text{[Kr]}4d^6$	-3974.815043	Pb	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^2$	-19523.831389
Li	$\text{[He]}2s^1$	-7.378133	Tc	$\text{[Kr]}4d^7$	-4204.141902	Bi	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^3$	-20095.328128
Be	$\text{[He]}2s^2$	-14.573023	Ru	$\text{[Kr]}4d^8$	-4441.038680	Po	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^4$	-20676.283998
B	$\text{[He]}2s^{2}2p^1$	-24.384693	Rh	$\text{[Kr]}4d^9$	-4685.600291	At	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^5$	-21266.749081
C	$\text{[He]}2s^{2}2p^2$	-37.344157	Pd	$\text{[Kr]}4d^{10}$	-4937.921024	Rn	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^6$	-21866.772241
N	$\text{[He]}2s^{2}2p^3$	-53.852155	Ag	$\text{[Kr]}4d^{10}5s^1$	-5197.639939	Fr	$\text{[Rn]}7s^1$	-22475.826522
O	$\text{[He]}2s^{2}2p^4$	-74.297532	Cd	$\text{[Kr]}4d^{10}5s^2$	-5465.133143	Ra	$\text{[Rn]}7s^2$	-23094.303666
F	$\text{[He]}2s^{2}2p^5$	-99.067145	In	$\text{[Kr]}4d^{10}5s^{2}5p^1$	-5740.082317	Ac	$\text{[Rn]}7s^{2}6d^1$	-23722.073196
Ne	$\text{[He]}2s^{2}2p^6$	-128.547098	Sn	$\text{[Kr]}4d^{10}5s^{2}5p^2$	-6022.746221	Th	$\text{[Rn]}7s^{2}6d^2$	-24359.362900
Na	$\text{[Ne]}3s^1$	-161.808533	Sb	$\text{[Kr]}4d^{10}5s^{2}5p^3$	-6313.211503	Pa	$\text{[Rn]}7s^{2}5f^3$	-25006.406325
Mg	$\text{[Ne]}3s^2$	-199.614636	Te	$\text{[Kr]}4d^{10}5s^{2}5p^4$	-6611.551696	U	$\text{[Rn]}7s^{2}5f^4$	-25663.398242
Al	$\text{[Ne]}3s^{2}3p^1$	-241.782323	I	$\text{[Kr]}4d^{10}5s^{2}5p^5$	-6917.837495	Np	$\text{[Rn]}7s^{2}5f^5$	-26330.321976
Si	$\text{[Ne]}3s^{2}3p^2$	-288.637472	Xe	$\text{[Kr]}4d^{10}5s^{2}5p^6$	-7232.138364	Pu	$\text{[Rn]}5f^{7}7s^1$	-27007.271797
P	$\text{[Ne]}3s^{2}3p^3$	-340.381142	Cs	$\text{[Xe]}6s^1$	-7553.899845	Am	$\text{[Rn]}5f^9$	-27694.356363
S	$\text{[Ne]}3s^{2}3p^4$	-397.202080	Ba	$\text{[Xe]}6s^2$	-7883.543827	Cm	$\text{[Rn]}5f^{10}$	-28391.573019
Cl	$\text{[Ne]}3s^{2}3p^5$	-459.286063	La	$\text{[Xe]}6s^{2}5d^1$	-8220.935691	Bk	$\text{[Rn]}5f^{11}$	-29098.977559
Ar	$\text{[Ne]}3s^{2}3p^6$	-526.817513	Ce	$\text{[Xe]}6s^{2}5d^2$	-8566.342481	Cf	$\text{[Rn]}5f^{12}$	-29816.624759
K	$\text{[Ar]}4s^1$	-599.124244	Pr	$\text{[Xe]}6s^{2}4f^{2}5d^1$	-8920.094872	Es	$\text{[Rn]}5f^{13}$	-30544.570349
Ca	$\text{[Ar]}4s^2$	-676.758186	Nd	$\text{[Xe]}6s^{2}4f^4$	-9282.434373	Fm	$\text{[Rn]}5f^{14}$	-31282.870930
Sc	$\text{[Ar]}4s^{2}4p^1$	-759.556762	Pm	$\text{[Xe]}6s^{2}4f^5$	-9653.359914	Md	$\text{[Rn]}5f^{14}7s^1$	-32031.135295
Ti	$\text{[Ar]}4s^{2}3d^2$	-847.933865	Sm	$\text{[Xe]}6s^{2}4f^6$	-10032.949725	No	$\text{[Rn]}5f^{14}7s^2$	-32789.512140
V	$\text{[Ar]}4s^{2}3d^3$	-942.147322	Eu	$\text{[Xe]}4f^{7}6s^2$	-10421.286649	Lr	$\text{[Rn]}5f^{14}7s^{2}6d^1$	-33557.812903
Cr	$\text{[Ar]}4s^{2}3d^4$	-1042.342957	Gd	$\text{[Xe]}6s^{1}4f^9$	-10818.487373	Rf	$\text{[Rn]}5f^{14}7s^{2}6d^2$	-34336.316816
Mn	$\text{[Ar]}3d^7$	-1148.803487	Tb	$\text{[Xe]}4f^{11}$	-11224.646666	Db	$\text{[Rn]}5f^{14}7s^{2}6d^3$	-35125.088022
Fe	$\text{[Ar]}3d^8$	-1261.579698	Dy	$\text{[Xe]}4f^{12}$	-11639.819030	Sg	$\text{[Rn]}5f^{14}6d^6$	-35924.293864
Co	$\text{[Ar]}3d^9$	-1380.817569	Ho	$\text{[Xe]}4f^{13}$	-12064.074984	Bh	$\text{[Rn]}5f^{14}6d^7$	-36733.871607
Ni	$\text{[Ar]}3d^{10}$	-1506.669759	Er	$\text{[Xe]}4f^{14}$	-12497.495312	Hs	$\text{[Rn]}5f^{14}6d^8$	-37553.863992
Cu	$\text{[Ar]}3d^{10}4s^1$	-1638.899667	Tm	$\text{[Xe]}4f^{14}6s^1$	-12939.976389	Mt	$\text{[Rn]}5f^{14}6d^9$	-38384.313942
Zn	$\text{[Ar]}3d^{10}4s^2$	-1777.848116	Yb	$\text{[Xe]}4f^{14}6s^2$	-13391.456193	Ds	$\text{[Rn]}5f^{14}6d^{10}$	-39225.264332
Ga	$\text{[Ar]}3d^{10}4s^{2}4p^1$	-1923.166449	Lu	$\text{[Xe]}4f^{14}6s^{2}6p^1$	-13851.687533	Rg	$\text{[Rn]}5f^{14}6d^{10}7s^1$	-40076.301440
Ge	$\text{[Ar]}3d^{10}4s^{2}4p^2$	-2075.150884	Hf	$\text{[Xe]}4f^{14}6s^{2}5d^2$	-14320.929628	Cn	$\text{[Rn]}5f^{14}6d^{10}7s^2$	-40937.797856
As	$\text{[Ar]}3d^{10}4s^{2}4p^3$	-2233.924574	Ta	$\text{[Xe]}4f^{14}6s^{2}5d^3$	-14799.321729	Nh	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^1$	-41809.456590
Se	$\text{[Ar]}3d^{10}4s^{2}4p^4$	-2399.595885	W	$\text{[Xe]}4f^{14}5d^6$	-15286.959470	Fl	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^2$	-42691.493680
Br	$\text{[Ar]}3d^{10}4s^{2}4p^5$	-2572.270918	Re	$\text{[Xe]}4f^{14}5d^7$	-15783.943765	Mc	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^3$	-43583.961779
Kr	$\text{[Ar]}3d^{10}4s^{2}4p^6$	-2752.054977	Os	$\text{[Xe]}4f^{14}5d^8$	-16290.259414	Lv	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^4$	-44486.902550
Rb	$\text{[Kr]}5s^1$	-2938.319660	Ir	$\text{[Xe]}4f^{14}5d^9$	-16805.965623	Ts	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^5$	-45400.354767
Sr	$\text{[Kr]}5s^2$	-3131.545686	Pt	$\text{[Xe]}4f^{14}5d^{10}$	-17331.121868	Og	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^6$	-46324.355815
Y	$\text{[Kr]}5s^{2}5p^1$	-3331.559557	Au	$\text{[Xe]}4f^{14}5d^{10}6s^1$	-17865.342083
Zr	$\text{[Kr]}5s^{2}4d^2$	-3538.662298	Hg	$\text{[Xe]}4f^{14}5d^{10}6s^2$	-18408.991495

Table 5. Table 5: Non-relativistic spin-restricted spherical HFS configurations for all elements in the periodic table. Entries in italic indicate a lower-lying configuration was identified but it failed to converge.

H	$1s^1$	-0.406534	Nb	$\text{[Kr]}4d^{3}5s^2$	-3747.428127	Tl	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^1$	-18948.496862
He	$1s^2$	-2.723640	Mo	$\text{[Kr]}5s^{1}4d^5$	-3969.125868	Pb	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^2$	-19510.422489
Li	$\text{[He]}2s^1$	-7.174881	Tc	$\text{[Kr]}4d^{6}5s^1$	-4198.246878	Bi	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^3$	-20081.732046
Be	$\text{[He]}2s^2$	-14.223291	Ru	$\text{[Kr]}4d^8$	-4434.888516	Po	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^4$	-20662.460965
B	$\text{[He]}2s^{2}2p^1$	-24.050406	Rh	$\text{[Kr]}4d^9$	-4679.115070	At	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^5$	-21252.645251
C	$\text{[He]}2s^{2}2p^2$	-37.053605	Pd	$\text{[Kr]}4d^{10}$	-4931.010033	Rn	$\text{[Xe]}4f^{14}5d^{10}6s^{2}6p^6$	-21852.321426
N	$\text{[He]}2s^{2}2p^3$	-53.567903	Ag	$\text{[Kr]}4d^{10}5s^1$	-5190.567420	Fr	$\text{[Rn]}7s^1$	-22461.201212
O	$\text{[He]}2s^{2}2p^4$	-73.925425	Cd	$\text{[Kr]}4d^{10}5s^2$	-5457.821825	Ra	$\text{[Rn]}7s^2$	-23079.470637
F	$\text{[He]}2s^{2}2p^5$	-98.456607	In	$\text{[Kr]}4d^{10}5s^{2}5p^1$	-5732.640932	Ac	$\text{[Rn]}7s^{2}5f^1$	-23707.189388
Ne	$\text{[He]}2s^{2}2p^6$	-127.490741	Sn	$\text{[Kr]}4d^{10}5s^{2}5p^2$	-6015.182678	Th	$\text{[Rn]}5f^{2}7s^2$	-24344.622650
Na	$\text{[Ne]}3s^1$	-160.628228	Sb	$\text{[Kr]}4d^{10}5s^{2}5p^3$	-6305.500906	Pa	$\text{[Rn]}5f^{3}7s^2$	-24991.833379
Mg	$\text{[Ne]}3s^2$	-198.248792	Te	$\text{[Kr]}4d^{10}5s^{2}5p^4$	-6603.649656	U	$\text{[Rn]}7s^{1}5f^5$	-25648.893676
Al	$\text{[Ne]}3s^{2}3p^1$	-240.346857	I	$\text{[Kr]}4d^{10}5s^{2}5p^5$	-6909.683446	Np	$\text{[Rn]}7s^{1}5f^6$	-26315.863733
Si	$\text{[Ne]}3s^{2}3p^2$	-287.145287	Xe	$\text{[Kr]}4d^{10}5s^{2}5p^6$	-7223.657213	Pu	$\text{[Rn]}7s^{1}5f^7$	-26992.780160
P	$\text{[Ne]}3s^{2}3p^3$	-338.804261	Cs	$\text{[Xe]}6s^1$	-7545.272707	Am	$\text{[Rn]}5f^{8}7s^1$	-27679.697021
S	$\text{[Ne]}3s^{2}3p^4$	-395.481609	Ba	$\text{[Xe]}6s^2$	-7874.734118	Cm	$\text{[Rn]}5f^{9}7s^1$	-28376.667807
Cl	$\text{[Ne]}3s^{2}3p^5$	-457.333996	La	$\text{[Xe]}6s^{2}4f^1$	-8212.148603	Bk	$\text{[Rn]}5f^{10}7s^1$	-29083.745568
Ar	$\text{[Ne]}3s^{2}3p^6$	-524.517426	Ce	$\text{[Xe]}6s^{2}4f^2$	-8557.852692	Cf	$\text{[Rn]}5f^{11}7s^1$	-29800.983007
K	$\text{[Ar]}4s^1$	-596.699051	Pr	$\text{[Xe]}4f^{3}6s^2$	-8911.927706	Es	$\text{[Rn]}5f^{12}7s^1$	-30528.432552
Ca	$\text{[Ar]}4s^2$	-674.160118	Nd	$\text{[Xe]}4f^{4}6s^2$	-9274.451612	Fm	$\text{[Rn]}5f^{13}7s^1$	-31266.146407
Sc	$\text{[Ar]}4s^{2}3d^1$	-757.000629	Pm	$\text{[Xe]}4f^{5}6s^2$	-9645.500832	Md	$\text{[Rn]}5f^{14}7s^1$	-32014.176598
Ti	$\text{[Ar]}4s^{2}3d^2$	-845.497930	Sm	$\text{[Xe]}4f^{6}6s^2$	-10025.150892	No	$\text{[Rn]}5f^{14}7s^2$	-32772.269829
V	$\text{[Ar]}3d^{3}4s^2$	-939.796100	Eu	$\text{[Xe]}4f^{7}6s^2$	-10413.476735	Lr	$\text{[Rn]}5f^{14}7s^{2}6d^1$	-33540.454380
Cr	$\text{[Ar]}3d^{4}4s^2$	-1040.034946	Gd	$\text{[Xe]}4f^{8}6s^2$	-10810.552897	Rf	$\text{[Rn]}5f^{14}6d^{2}7s^2$	-34318.854809
Mn	$\text{[Ar]}4s^{1}3d^6$	-1146.366756	Tb	$\text{[Xe]}4f^{9}6s^2$	-11216.453617	Db	$\text{[Rn]}5f^{14}6d^{4}7s^1$	-35107.525943
Fe	$\text{[Ar]}4s^{1}3d^7$	-1258.917212	Dy	$\text{[Xe]}4f^{10}6s^2$	-11631.252911	Sg	$\text{[Rn]}5f^{14}6d^6$	-35906.506548
Co	$\text{[Ar]}4s^{1}3d^8$	-1377.819755	Ho	$\text{[Xe]}4f^{11}6s^2$	-12055.024619	Bh	$\text{[Rn]}5f^{14}6d^7$	-36715.824635
Ni	$\text{[Ar]}3d^{9}4s^1$	-1503.210775	Er	$\text{[Xe]}4f^{12}6s^2$	-12487.842443	Hs	$\text{[Rn]}5f^{14}6d^8$	-37535.505151
Cu	$\text{[Ar]}3d^{10}4s^1$	-1635.226377	Tm	$\text{[Xe]}4f^{13}6s^2$	-12929.779972	Mt	$\text{[Rn]}5f^{14}6d^9$	-38365.584348
Zn	$\text{[Ar]}3d^{10}4s^2$	-1773.909886	Yb	$\text{[Xe]}4f^{14}6s^2$	-13380.910702	Ds	$\text{[Rn]}5f^{14}6d^{10}$	-39206.098757
Ga	$\text{[Ar]}3d^{10}4s^{2}4p^1$	-1919.085911	Lu	$\text{[Xe]}4f^{14}6s^{2}5d^1$	-13840.976253	Rg	$\text{[Rn]}5f^{14}6d^{10}7s^1$	-40056.951158
Ge	$\text{[Ar]}3d^{10}4s^{2}4p^2$	-2070.946515	Hf	$\text{[Xe]}4f^{14}6s^{2}5d^2$	-14310.121254	Cn	$\text{[Rn]}5f^{14}6d^{10}7s^2$	-40918.195130
As	$\text{[Ar]}3d^{10}4s^{2}4p^3$	-2229.571620	Ta	$\text{[Xe]}4f^{14}6s^{2}5d^3$	-14788.392156	Nh	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^1$	-41789.700671
Se	$\text{[Ar]}3d^{10}4s^{2}4p^4$	-2395.043625	W	$\text{[Xe]}4f^{14}6s^{1}5d^5$	-15275.846800	Fl	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^2$	-42671.589032
Br	$\text{[Ar]}3d^{10}4s^{2}4p^5$	-2567.446685	Re	$\text{[Xe]}4f^{14}5d^{6}6s^1$	-15772.541265	Mc	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^3$	-43563.886976
Kr	$\text{[Ar]}3d^{10}4s^{2}4p^6$	-2746.866101	Os	$\text{[Xe]}4f^{14}5d^8$	-16278.531177	Lv	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^4$	-44466.621119
Rb	$\text{[Kr]}5s^1$	-2932.972209	Ir	$\text{[Xe]}4f^{14}5d^9$	-16793.845129	Ts	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^5$	-45379.818244
Sr	$\text{[Kr]}5s^2$	-3125.998090	Pt	$\text{[Xe]}4f^{14}5d^{10}$	-17318.533845	Og	$\text{[Rn]}5f^{14}6d^{10}7s^{2}7p^6$	-46303.505356
Y	$\text{[Kr]}5s^{2}4d^1$	-3325.964742	Au	$\text{[Xe]}4f^{14}5d^{10}6s^1$	-17852.550237
Zr	$\text{[Kr]}5s^{2}4d^2$	-3533.076869	Hg	$\text{[Xe]}4f^{14}5d^{10}6s^2$	-18395.920112

Equations

$${\chi }_{nlm}=r^{-1}{B}_n(r)Y_l^m(\theta ,\varphi )$$

$$\frac{1}{r_{12}}=$$

$$(ij|kl)=$$

$$\frac{1}{r}=\frac{{\varphi }_{\text{sr}}(r)}{r}+\frac{1-{\varphi }_{\text{sr}}(r)}{r}=\frac{{\varphi }_{\text{sr}}(r)}{r}+\frac{{\varphi }_{\text{lr}}(r)}{r}=\frac{1-{\varphi }_{\text{lr}}(r)}{r}+\frac{{\varphi }_{\text{lr}}(r)}{r},$$

$$\frac{{\varphi }_{\mathrm{sr}}(r)}{r_{12}}$$

$$\times \sum _{M=-L}^L{\left(Y_L^M({\Omega }_1)\right)}^{\ast }{Y}_L^M({\Omega }_2)$$

$${\mathcal{G}}_L^{\text{Coulomb}}(r_{>},r_{<})=\frac{r_{<}^L}{r_{>}^{L+1}}.$$

$$I_{ijkl}^L$$

$$\times B_k(r_2)B_l(r_2){\mathcal{G}}_L(r_{>},r_{<},\mu ),$$

$$\frac{e^{-\lambda |\mathbf{r}-{\mathbf{r}}^{\prime }|}}{|\mathbf{r}-{\mathbf{r}}^{\prime }|}$$

$$\times \sum _{M=-L}^L{Y}_L^M({\Omega }_1){\left(Y_L^M({\Omega }_2)\right)}^{\ast }$$

$${\mathcal{G}}_L^{\text{Yukawa}}(r_{>},r_{<},\lambda )=(2L+1)\lambda i_L(\lambda r_{<})k_L(\lambda r_{>}).$$

$${\mathcal{G}}_L(R,r;\mu )=\mu {\Phi }_L(\Xi ,\xi )$$

$${\Phi }_n(\Xi ,\xi )$$

$$F_n(\Xi ,\xi )$$

$$H_n(\Xi ,\xi )$$

$${\Phi }_n(\Xi ,\xi )$$

$$D_{n,0}(\Xi )$$

$$D_{n,k}(\Xi )$$

$${\psi }_{nlm}(\mathbf{r})=R_{nl}(r)Y_l^m(\hat{\mathbf{r}}).$$

$$n(\mathbf{r})$$

$$=\sum _{nl}{f}_{nl}\frac{\left(2l+1\right)R_{nl}^2(r)}{4\pi }=n(r)$$

$$P_{\mu \nu }^{\sigma }={\delta }_{l_{\mu },l_{\nu }}{\delta }_{m_{\mu },m_{\nu }}{P}_{\mu \nu }^{l_{\mu };\sigma }$$

$$\nabla n={\partial }_{r}n{\hat{\mathbf{e}}}_r$$

$$K_{\mu \nu }^{xc;\sigma }=$$

$$K_{\mu \nu }=\sum _{\sigma \tau }(\mu \sigma |\nu \tau )P_{\sigma \tau }.$$

$$K_{\mu \nu }^{l_{\text{out}}}$$

$$C=-\frac{1}{2Z}\frac{{\rm d}\log n(r)}{{\rm d}r}\Bigg{|}_{r=0}=-\frac{1}{2Z}\frac{n^{\prime}(0)}{n(0)}$$

$$n(0)$$

$$=P_{\mu \nu }\underset{r\to 0}{\lim }\frac{\frac{{\mathrm{d}}^2}{{\mathrm{dr}}^2}{B}_{\mu }(r)B_{\nu }(r)}{\frac{{\mathrm{d}}^2}{{\mathrm{dr}}^2}{r}^2}$$

$$=P_{\mu \nu }{B}_{\mu }^{\prime }(r)B_{\nu }^{\prime }(r).$$

$$n^{\prime }(0)$$

$$=P_{\mu \nu }\underset{r\to 0}{\lim }\left[\frac{\frac{{\mathrm{d}}^3}{{\mathrm{dr}}^3}{B}_{\mu }(r)B_{\nu }(r)}{\frac{{\mathrm{d}}^2}{{\mathrm{dr}}^2}{r}^2}-2\frac{\frac{{\mathrm{d}}^3}{{\mathrm{dr}}^3}{B}_{\mu }(r)B_{\nu }(r)}{\frac{{\mathrm{d}}^3}{{\mathrm{dr}}^3}{r}^3}\right]=P_{\mu \nu }{B}_{\mu }^{\prime \prime }(0)B_{\nu }^{\prime }(0).$$

$$V(\mathbf{r})$$

$$V_{ij}(r)$$

Full text

The implementation of range-separation based on the Yukawa or complementary error function (erfc) kernels is also described, allowing complete basis set benchmarks of modern range-separated hybrid functionals with either integer or fractional occupation numbers. Finally, computation of atomic effective potentials at the local density or generalized gradient approximation levels for the superposition of atomic potentials (SAP) approach [S. Lehtola, J. Chem. Theory Comput. 15, 1593 (2019)] that has been shown to be a simple and efficient way to initialize electronic structure calculations is described.

The present numerical approach is shown to afford beyond microhartree accuracy with a small number of numerical basis functions, and to reproduce literature results for the ground states of atoms and their cations for $1\leq Z\leq 86$. Our results indicate that the literature values deviate by up to 10 $\mu E_{h}$ from the complete basis set limit. The numerical scheme for the erfc kernel is shown to work by comparison to results from large Gaussian basis set calculations from the literature. Spin-restricted ground states are reported for Hartree–Fock and Hartree–Fock–Slater calculations with fractional occupations for $1\leq Z\leq 118$.

I Introduction

Atoms are the simplest possible unit in chemistry, which is why electronic structure studies on atoms have a long and venerated history. Thanks to the high amount of symmetry that may be used to reduce the number of degrees of freedom in the atomic problem, fully numerical electronic structure approaches on atoms have been possible for a very long time (Lehtola, 2019a); for instance, a fully numerical configuration interaction calculation on the oxygen atom was reported by Hartree and coworkers over 80 years ago (Hartree et al., 1939).

As the atomic hamiltonian is spherically symmetric, the exact wave function should be rotationally invariant as well. Although the necessary symmetry requirements can straightforwardly be enforced in wave function approaches, the application of density functional theory (Hohenberg and Kohn, 1964; Kohn and Sham, 1965) (DFT) on atoms is surprisingly tricky. In the usual DFT approach, a single Slater determinant is employed, with all orbitals below the Fermi level being fully occupied. Non-relativistically, all $2l+1$ atomic orbitals sharing the principal quantum number $n$ and angular quantum number $l$ should be completely degenerate; however, this behavior is broken by conventional DFT as well as Hartree–Fock (HF) already on the first row. Different choices for the occupied orbitals on the $2p$ shell yield different final energies for e.g. B and F, which may lead to several kcal/mol differences in the total energy – with a symmetrized density yielding yet another result (Baerends et al., 1997). One possibility to obtain comparable results is to employ a standard set of electronic occupations (Hay, 1977), but such an approach does not yield the lowest possible energy.

Pursuing the lowest energy is not unproblematic, either. While HF is infamous for possessing variational solutions that break symmetry in systems with a high degree of symmetry (Prat, 1972), symmetry breaking is a problem in DFT as well (Goursot et al., 1995). In atoms, broken symmetries often arise for open shells, and the effect of non-spherical densities is known to be more pronounced with functionals at the generalized gradient approximation (GGA) and especially the meta-GGA (mGGA) level than at the local density approximation (LDA) level (Kutzler and Painter, 1987; Philipsen and Baerends, 1996; Tao and Perdew, 2005; Johnson et al., 2007); even optimized effective potential exact-exchange calculations are subject to spurious energy splittings (Pittalis et al., 2006). Inclusion of current density dependence leads to improvement of GGA and mGGA results (Tao and Perdew, 2005; Johnson et al., 2007), but the proper orbital degeneracy is still not fully restored.

Symmetry breaking effects in atoms can be seen already at the simplest possible level of DFT, that is, the exchange-only LDA, which is also commonly known as Hartree–Fock–Slater (HFS) theory. For example, HFS calculations on the F atom reveal millihartree decreases of the total energy upon addition of $d$ as well as $f$ functions, which is at variance to the generally accepted electronic configuration of fluorine as $1s^{2}2s^{2}2p^{5}$. Interestingly, this kind of symmetry breaking sometimes happens even in the case of closed-shell atoms; see, for instance, our recent finite element reproduction (Lehtola, 2019b) of calculations on atomic anions (Anderson et al., 2017) where symmetry breaking was observed for \ceH-, Be, \ceLi-, and \ceNa-.

In addition to being degenerate due to symmetry (as often in atoms), orbitals may also be degenerate by accident. Since the aufbau rule implies populating the orbitals in increasing energy, it tempting to divide the occupations evenly in the case of degeneracies. This paves the way to the use of fractional occupations, which in the case of atoms naturally yield a spherically symmetric density thanks to Unsöld’s theorem (Unsöld, 1927); the use of fractional occupations can be formally justified within the theory of ensemble representable densities (Englisch and Englisch, 1984a, b).

Fractionally occupied orbitals should especially be used in the case where there is a negative gap between the highest occupied and lowest unoccupied orbital no matter which way the orbitals are occupied; this happens when the highest occupied and lowest unoccupied orbital switch places during the orbital optimization. In this case, the total energy can be lowered by moving a fraction of an electron from the highest occupied orbital to the lowest unoccupied orbital, and at some point the two levels should cross.

Fractional occupations have been shown to yield better results for strongly correlated systems (Dunlap and Mei, 1983; Wang and Schwarz, 1996; Schipper et al., 1998; Takeda et al., 2003; Nygaard and Olsen, 2013a). However, fractional occupations can only be justified at the Fermi level (Valiev and Fernando, 1995), and more recently it has been shown that energy minimization naturally leads to integer occupations below the Fermi level, and possible fractional occupations at the Fermi level for independent particle models like HF and DFT (Giesbertz and Baerends, 2010).

While in some systems it is clear a priori from symmetry arguments or the orbital energies how many orbitals should be fractionally occupied, this is generally not the case. However, fractional occupations can be obtained as (Kraisler et al., 2009) the zero-temperature limit of finite-temperature DFT (FT-DFT) (Mermin, 1965; Stoitsov and Petkov, 1988). In a finite-temperature approach, the fractional orbital occupation numbers are determined by the orbital energies according to some smearing scheme that is typically controlled by a single parameter, an electronic temperature. Because of the simplicity and favorable computational scaling of FT-DFT, it has become a powerful tool for approximate modeling of systems exhibiting strong correlation; such approaches have been used to obtain promising results for a variety of systems (Chai, 2012, 2014; Wu and Chai, 2015; Yeh and Chai, 2016; Seenithurai and Chai, 2016; Wu et al., 2016; Chai, 2017; Seenithurai and Chai, 2017; Lin et al., 2017; Yeh et al., 2018; Seenithurai and Chai, 2018; Deng and Chai, 2019; Chung and Chai, 2019; Grimme, 2013; Grimme and Hansen, 2015).

Finite electronic temperatures may also be used to aid the convergence of self-consistent field calculations of molecules (Rabuck and Scuseria, 1999); in the solid state, the use of fractional occupation numbers is often mandatory in order to attain convergence (Kratzer and Neugebauer, 2019). Although finite temperature approaches are more attractive for DFT where all electrons experience the same potentials, finite temperature approaches can also be used in the context of HF calculations where they may offer good active spaces for post-HF calculations on strongly correlated systems (Slavíček and Martínez, 2010).

Although several types of smearing schemes have been suggested, including Fermi–Dirac (Mermin, 1965), Gaussian smearing (Fu and Ho, 1983), Methfessel–Paxton smearing (Methfessel and Paxton, 1989), cold smearing (Marzari et al., 1999), and others (Holender et al., 1995), they have been shown to yield similar results if the parameters are adjusted properly (Springborg et al., 1998; Grotheer and Fähnle, 1998; Slavíček and Martínez, 2010); however, the behavior with respect to temperature needs to be carefully checked in each case to ensure convergence (Mehl, 2000). Note that the evaluation of forces in finite-temperature calculations require the consideration of an additional entropic term that arises from the non-integer occupations and that depends on the smearing function (Weinert and Davenport, 1992; Warren and Dunlap, 1996).

Regardless of the used temperature, calculations with fractional occupations are more involved than those with integer occupations. Convergence acceleration techniques such as direct inversion in the iterative subspace (Pulay, 1980, 1982) (DIIS) become invalid when the orbital occupation pattern changes, even though the self-consistent field problem itself may become easier with fractional occupation numbers (Rabuck and Scuseria, 1999). Determining the correct occupations is hard, since the orbital occupations depend on the orbital energies, which in turn depend on the orbital occupations. The changes in the occupations may also cause changes in the shapes of the orbitals, meaning that the orbitals, their energies and their occupations need to be solved self-consistently. Several approaches have been proposed for solving this problem both for zero (Averill and Painter, 1992; Cancès et al., 2003; Kraisler et al., 2009; Nygaard and Olsen, 2013b) and finite electronic temperatures (Gillan, 1989; Fernando et al., 1989; Grumbach et al., 1994; Chetty et al., 1995; Holender et al., 1995; Marzari et al., 1997).

In systems with a high degree of symmetry such as atoms, the fractional occupations can be defined by symmetry block. Fractional occupations for atoms are typically defined in terms of atomic shells, over which the electrons are equally divided. For instance, the $1s^{2}2s^{2}2p^{5}$ configuration for F implies that the hole in the $2p$ shell be equally divided, resulting in the minority spin occupations $2p_{x}^{2/3}2p_{y}^{2/3}2p_{z}^{2/3}$; a spin-restricted variant would employ occupations of $2p_{x}^{5/6}2p_{y}^{5/6}2p_{z}^{5/6}$ in both spin channels. Indeed, this is the method of choice for fully numerical density functional calculations on atoms (Lehtola, 2019a), and it has been used *e.g. *in ref. 66 for local density calculations on $1\leq Z\leq 92$ at the ground state electronic configuration from experiment, and in ref. 67 for Perdew–Burke–Ernzerhof (PBE) (Perdew et al., 1996, 1997) calculations on $1\leq Z\leq 20$ and $31\leq Z\leq 36$.

Atomic calculations with fractional occupation numbers are also typically used to generate pseudopotentials (Fuchs and Scheffler, 1999; Oliveira and Nogueira, 2008), numerical atomic orbital basis functions (Ozaki and Kino, 2004; Blum et al., 2009), and Gaussian basis sets (Andzelm et al., 1985; Godbout et al., 1992; Porezag and Pederson, 1999). Spin-restricted spherically symmetric atoms may also be used for setting up frozen core calculations within all-electron approaches, and to determine approximate binding energies (te Velde et al., 2001). We have also recently shown that the radial potential from atomic calculations with fractional occupation numbers can be used to formulate efficient initial guesses for electronic structure calculations on polyatomic systems via the superposition of atomic potentials (SAP) approach (Lehtola, 2019c).

In the typical case, electrons are divided evenly among the $2l+1$ orbitals that are degenerate by symmetry. However, the fractional occupations can be generalized beyond integer occupations per shell, in case accidental degeneracy is also present. Early multiconfigurational HF calculations on atoms found that the $3d$ orbitals become occupied before the $4s$ orbitals in transition metals (Slater et al., 1969; Abdulnur et al., 1972), which was solved by moving fractions of an electron between the shells. One example of this approach is the iron atom, where the [Ar]$3d_{1}^{5}4s_{1}^{1}$ and [Ar]$3d_{2}^{5}4s_{0}^{1}$ configurations both turn out to have a negative gap in the local-density approximation (Janak, 1978), the upper and lower indices denoting spin-up and spin-down electrons, respectively. With the Vosko–Wilk–Nusair (VWN) local density functional, the lowest-energy configuration is found to be [Ar]$3d_{1.3984}^{5}s_{0.602}^{1}$ (Kraisler et al., 2009).

A systematic, non-relativistic study for spherical atoms $1\leq Z\leq 86$ has recently been presented by Kraisler, Makov and Kelson for the local density and PBE functionals based on three local density functionals, employing 16 000 point grids and wave functions converged to 2 $\mu E_{h}$ (Kraisler et al., 2010). It was found in ref. 82 that the ground state of most atoms does not involve fractional splitting of electrons between shells, indicating that a fully numerical program for modeling atoms with spherical densities would go a long way towards the final solution.

While several programs exist for either wave function or density functional based fully numerical calculations on atoms (Lehtola, 2019a), we are not aware of any publicly available software that supports hybrid functionals, except the recently published HelFEM program (Lehtola, 2019b, 2018), which also includes a fully numerical approach for diatomic molecules that similarly supports hybrid functionals (Lehtola, 2019d). Most publicly available programs for fully numerical density functional calculations on atoms target the generation of projector-augmented wave (PAW) setups (Blöchl, 1994) or the generation of pseudopotentials (Schwerdtfeger, 2011). Although Hartree–Fock pseudopotential generators have been available for some time (Trail and Needs, 2005; Al-Saidi et al., 2008), which allowed the use of non-self-consistent pseudopotentials for hybrid functionals (Wu et al., 2009), surprisingly, the self-consistent generation of pseudopotentials for hybrid functionals has only been described last year (Yang et al., 2018), explaining the scarcity of such programs.

Interestingly, the work of Yang et al in ref. 90 did not employ fractionally occupied Hartree–Fock calculations, but rather followed Slater’s multiconfigurational approach, which is at odds with the density functional description used in the work, as the exact exchange and density functional parts experience different electron densities. In contrast, when fractional occupations are employed as in the present work, the exchange exact operator becomes independent of the magnetic quantum number $m$ as will be shown in II.2, and both the density functional and exact exchange operators are evaluated with the same density matrix.

Although a general-use atomic program like the one in HelFEM can be straightforwardly adapted to calculations on spherically symmetric densities by employing fractional occupation numbers in the construction of the density matrix, a more efficient approach is afforded by taking the assumption of the spherical symmetry of the density matrix deeper in the algorithms. As a result, some or even all of the angular integrals can be eliminated from the calculations, reducing the problem to a small number of dimensions; indeed, this is* *exactly what is done in the multiconfigurational HF approach Slater proposed 90 years ago (Slater, 1929).

In the present work, we describe the extension of the atomic program in HelFEM to the description of atoms with spherical symmetric density via fractional occupation numbers. Alike the other programs in HelFEM, the spherically symmetric atomic program is interfaced to the libxc library of density functionals (Lehtola et al., 2018) and can be used with all supported density functionals therein. Specialized implementations for atomic calculations with fractional occupations are developed for local density (LDA) and generalized gradient (GGA) functionals as well as HF exchange, yielding significant reductions in the dimensionality of the problem, whereas meta-GGA functionals can be used via an interface to the algorithms previously developed in ref. 14.

Importantly, we also describe the implementation of Yukawa and complementary error function (erfc) range-separated exchange for atomic calculations in HelFEM with either fractional or integer occupations, allowing complete basis set benchmarks of recently developed exchange-correlation functionals such as the CAM-QTP family by Bartlett and coworkers (Verma and Bartlett, 2014; Jin and Bartlett, 2016; Haiduke and Bartlett, 2018), the N12-SX and revM11 functionals by Truhlar and coworkers (Peverati and Truhlar, 2012; Verma et al., 2019), and the $\omega$B97X-V and $\omega$B97M-V functionals by Mardirossian and Head-Gordon (without the non-local correlation part) (Mardirossian and Head-Gordon, 2014, 2016). While the spherical harmonics decomposition for the Yukawa kernel is well known, the decomposition for the erfc kernel has only been derived some time ago (Ángyán et al., 2006) and does not appear to have been implemented within a generally applicable fully numerical approach for atoms. Results for H and He with relatively low-order B-spline basis sets have, however, been published almost simultaneously to our work (Zapata et al., 2019). Finally, we also describe the analytic calculation of the radial potentials necessary for the SAP orbital guess (Lehtola, 2019c).

In the next section, we derive the equations for fractionally occupied HF and DFT at the LDA and GGA levels within the used finite element approach. Then, in the Results section, we present applications of the program to reproducing ground states for the neutral atoms and cations $1\leq Z\leq 86$ and compare with ref. 82; we reproduce the long-range corrected density functional calculations on closed-shell atoms of ref. 15 to show that the range-separation scheme works; and finally we report the non-relativistic ground states of all atoms in the periodic table at HF and HFS levels of theory. The article concludes with a brief summary and discussion section.

II Method

A basis set of the form

[TABLE]

is adopted as in the integer-occupation program described in ref. 14. Here, $B_{n}(r)$ are the piecewise polynomial shape functions of the finite element method, which have been discussed extensively in refs. 1 and 14 to which we refer for further details.

II.1 Range-separated exchange

As discussed in refs. 1 and 14, the key to fully numerical electronic structure calculations on atoms is the Laplace expansion

[TABLE]

that factorizes the two-electron integrals

[TABLE]

into a radial and an angular part.

In range-separated density functional theory (Gill et al., 1996; Leininger et al., 1997), the Coulomb interaction is split into a short-range (sr) and a long-range (lr) part as

[TABLE]

where $\phi_{{\rm\text{sr}}}(r)=1-\phi_{\text{lr}}(r)$ is a splitting function. Typically, the short-range part is described using density functional theory, and the long-range part with HF theory, but in practice many functionals employ more flexibility: for instance, the CAM-B3LYP functional (Yanai et al., 2004) contains 19% short-range and 65% long-range exact exchange.

The evaluation of the range-separated exchange functionals is simple if one has access to the Green’s function expansion of the range-separated kernel as

[TABLE]

where $\mathcal{G}_{L}(r_{>},r_{<},\mu)$ is the Green’s function, $r_{>}$ and $r_{<}$ are the greater and smaller of $r_{1}$ and $r_{2}$, respectively, and $\mu$ is the range separation parameter. The Green’s function for the (unscreened) classical Coulomb interaction can be identified from (LABEL:laplace) as

[TABLE]

The implementation of the integrals in HelFEM is based on the primitive integrals defined in ref. 14 as

[TABLE]

where $B_{i}(r)$ are the piecewise polynomial basis functions of (LABEL:basis).

II.1.1 Yukawa kernel

The Yukawa-screened (Yukawa, 1935) potential, $\phi^{\text{sr}}(r_{12})=\exp(-\lambda r_{12})$ has a relatively well-known simple expansion

[TABLE]

where $i_{L}$ and $k_{L}$ are regular and irregular modified Bessel functions that are regular at zero and infinity, respectively. Due to its separability, Yukawa-screened functionals are easy to handle in fully numerical approaches. Indeed, the Yukawa Green’s function is employed in several recently developed linear scaling approaches for solving the HF or Kohn–Sham equations for bound orbitals in molecular systems via the Helmholtz kernel (Harrison et al., 2004; Frediani et al., 2013; Jensen, 2014; Solala et al., 2017; Parkkinen et al., 2017). The Yukawa interaction is also straightforward to implement in calculations with Slater-type orbitals (Seth and Ziegler, 2012; Seth et al., 2013; Rico et al., 2013). It turns out that Yukawa screening can also be implemented with Gaussian-type orbitals in a rather straightforward manner (Akinaga and Ten-no, 2008), as analogous integrals also arise within $r_{12}$ wave function theory (Ten-no, 2004, 2007). Such implementations are, however, rare at the moment, even though it has been claimed that Yukawa screening yields more accurate atomization and charge transfer excitation energies than erfc screening (Akinaga and Ten-No, 2009). The Green’s function for the Yukawa interaction can be read from (LABEL:sr-yukawa) as

[TABLE]

As the Yukawa interaction factorizes in $r_{>}$ and $r_{<}$, it can be implemented in a similar fashion to the full Coulomb interaction, (LABEL:Gl-coulomb), along the lines of ref. 14.

II.1.2 erfc kernel

Most range-separated functionals, however, are based on the complementary error function (erfc) kernel $\phi_{\text{sr}}(r)=\text{erfc }(\mu r)$. Such functionals are easy to implement in Gaussian-basis programs, requiring but simple modifications to the two-electron integrals (Heyd et al., 2003; Ahlrichs, 2006), as well as plane wave programs since the kernel has a simple Fourier transform which is strongly attenuated at large momentum. In contrast, the implementation of the erfc kernel is more complicated in real-space approaches. Fortunately, spherical harmonic expansions for the erfc Green’s functions are available in the literature (Marshall, 2002; Ángyán et al., 2006), but their form is more involved than that of the Yukawa function in (LABEL:sr-yukawa). The main complication is that the Green’s function does not factorize in $r_{<}$ and $r_{>}$, which means that two-dimensional quadrature is always required. In the approach of ref. 100, new variables are introduced as $\Xi=\mu R$ and $\xi=\mu r$ and

[TABLE]

where $\Phi_{L}$ is a scaled radial function given by

[TABLE]

(Note that the lower limit of the sum in (LABEL:Fn) is incorrect in ref. 100, where it reads $p=1$ instead of $p=0$.) \EqrangerefphinHn are numerically unstable in the short range, which is why when either $\xi<0.4$, or $\Xi<0.5$ and $0<\xi<2\Xi$ (Ángyán et al., 2006), the Green’s function is evaluated with a Taylor expansion

[TABLE]

Despite the lack of factorization of the erfc Green’s function, its evaluation can be carried out analogously to the Coulomb and Yukawa kernels. The primitive integrals, (LABEL:primint), can be divided into two cases thanks to the finite support of the piecewise polynomial basis functions, as discussed in ref. 14. In an intraelement integral, both $ij$ and $kl$ are within the same element, whereas in an interelement integral $ij$ are in one element and $kl$ are in another. In analogy to the scheme for Coulomb integrals discussed in ref. 14, the interelement integrals are evaluated with $N_{\text{quad}}$ quadrature points in both $ij$ and $kl$, whereas the intraelement integrals employ $N_{\text{quad}}$ points in $ij$, whereas the $kl$ quadrature is split into $N_{\text{quad}}$ intervals, all of which employ a fresh set of $N_{\text{quad}}$ quadrature points.

II.2 Self-consistent field calculations with fractional occupations

It is well known that atomic orbitals can be written in the form

[TABLE]

Employing smeared occupations as

[TABLE]

where $f_{nl}$ is the occupation number of all the $2l+1$ orbitals on the $(n,l)$ shell, one immediately sees that the density matrix is diagonal in $l$ and $m$

[TABLE]

and that the elements of the density matrix only depend on the value of $l$.

The spherical averaging yields huge simplifications for density functional calculations. As now the density is only a function of the radial coordinate, also its gradient

[TABLE]

only depends on the radial coordinate. Following the usual projective approach (Pople et al., 1992; Lehtola, 2019b), the LDA and GGA matrix elements

[TABLE]

become greatly simplified as only the radial terms are picked up, and as the same radial basis is used for all $l,m$; see (LABEL:basis). Note, however, that meta-GGAs that depend on the kinetic energy density cannot be handled in the same fashion, as the kinetic energy density is not manifestly dependent only on the radial coordinate as discussed e.g. in ref. 122. Alike the exact exchange discussed below, the meta-GGA potential turns out to depend on the $l$ channel. Meta-GGA functionals can be used in the present program via a fractional-occupation interface to the full atomic routines discussed in ref. 14.

The Coulomb matrix arising from (LABEL:laplace) trivially reduces to a single term as the spherically symmetric density only consists of a single $L=0$, $M=0$ component. Exact exchange – either with the full Coulomb form of (LABEL:Gl-coulomb) or the range-separated versions in (LABEL:Gl-Yukawa,_Gl-erfc) – is a bit more complicated, as both the integrals and the density matrix carry a dependence on the orbital angular momenta in the well-known equation

[TABLE]

Employing the blocking of the density matrix given in (LABEL:densmat), the exchange matrix can be written as

[TABLE]

where $L$ is a coupled angular momentum with $z$ projection $M=m_{\text{out}}-m_{\text{in}}$. Rearranging the contractions, it is then seen that

[TABLE]

where the evaluation is done from the insidemost bracket out.

II.3 Cusp condition

One way to diagnose atomic wave functions is the Kato–Steiner cusp condition (Kato, 1957; Steiner, 1963)

[TABLE]

which yields the value $C=1$ for the exact HF or density functional solution (Nagy and Sen, 2001). The electron density $n(r)$ at the nucleus was obtained in ref. 14 via l’Hôpital’s rule as

[TABLE]

Its derivative at the nucleus also turns out to have a simple expression:

[TABLE]

The value of the cusp is printed out at the end of all atomic calculations in HelFEM.

II.4 Effective radial potential for SAP

In the SAP approach discussed in ref. 78, approximate orbitals for a molecule are obtained by diagonalizing an effective one-body Hamiltonian in an external potential obtained as a superposition of radial atomic potentials. Once the atomic ground state has been found with any supported method in HelFEM, including HF and hybrid and meta-GGA functionals, the radial effective potential for the SAP approach can be calculated based on any LDA or GGA functional. Extensions to the exact exchange, as in the optimized effective potential method (Sharp and Horton, 1953), as well as generalized Kohn–Sham methods for the radial potentials from meta-GGA functionals are left for future work.

If the radial potential is self-consistent, i.e. the same functional was used for both the atomic orbitals and the potential, the SAP guess will reproduce the atomic orbitals exactly (Lehtola, 2019c). The atomic potential comprises Coulomb and exchange-correlation contributions, the calculation of which is presented in the following.

II.4.1 Coulomb potential

Employing the Laplace expansion, (LABEL:laplace), the Coulomb potential at a point** $\boldsymbol{r}$** for a spherically symmetric charge distribution is

[TABLE]

Expressing the orbitals as in (LABEL:atorb) yields potential matrix elements of the form

[TABLE]

and one gets three cases depending on whether $r$ is inside the element where $i$ and $j$ reside, or not. Let the element begin at $r_{b}$ and end at $r_{e}$. Now

[TABLE]

Like the two-electron integrals discussed above, the in-element potential $r_{b}<r<r_{e}$ has to be evaluated by slices at every radial quadrature point $(r_{0},r_{1},\dots,r_{n-1})$

[TABLE]

II.4.2 Exchange-correlation potential

The functional derivative satisfies

[TABLE]

and so

[TABLE]

Integrating by parts one gets

[TABLE]

from which one can identify

[TABLE]

Expressing the functional in terms of

[TABLE]

one has

[TABLE]

and so

[TABLE]

or for an open shell system

[TABLE]

where $\sigma\neq\sigma^{\prime}$.

To guarantee accuracy, the gradient terms have to be evaluated analytically. Fortunately, there’s only radial dependence, so the gradient

[TABLE]

can be replaced by a radial derivative, and the divergence with

[TABLE]

Now,

[TABLE]

where

[TABLE]

where we have defined $g^{\tau}=\partial_{r}n^{\tau}$ and $l^{\tau}=\partial_{r}^{2}n^{\tau}$, and the extra $2A_{r}/r$ term from the divergence, (LABEL:divA), yielding

[TABLE]

Thus, altogether, the radial exchange(-correlation) potential is given by

[TABLE]

where the various derivatives of the exchange-correlation functional are available in libxc (Lehtola et al., 2018).

III Results

To demonstrate the new routines, we reproduce literature values for the ground states of the neutral and cationic atoms $1\leq Z\leq 86$ with the VWN functional, as well as a Perdew–Burke–Ernzerhof (PBE) (Perdew et al., 1996, 1997) functional based on VWN correlation (Kraisler et al., 2010), accurate values for which were obtained in ref. 82 with 10 000 radial grid points. We found that by using a radial basis consisting of ten 15-node elements and a practical infinity $r_{\infty}=40a_{0}$, the energy has converged beyond microhartree accuracy, even though this basis contains just 139 radial basis functions – almost two orders of magnitude fewer degrees of freedom than used in ref. 82. These results, with differences to the reference data from ref. 82 are shown in 1 for neutral atoms and 2 for their cations.

The agreement to most part is excellent: large positive differences that indicate the value of ref. 82 is lower are seen for the species for which the calculations in ref. 82 transferred fractional charge across shells. Otherwise, the differences become noticeable for heavier atoms, nearing $-10\,\mu E_{h}$ when $Z\to 86$, indicating that the data in ref. 82 is not fully converged to the basis set limit. In our VWN calculation on the \ceLa+ cation, it was discovered that the energy for the $4f_{0}^{0.945}5d_{0}^{0.575}6s_{0}^{0.480}$ state reported in ref. 82 was incorrect; the correct energy is $-8217.456974$ \bibnoteEli Kraisler, private communication, 2019..

Next, we demonstrate that the erfc range-separation scheme works by reproducing literature values for the total energies of the spherically symmetric atoms on the first two rows using the long-range corrected BLYP functional (Iikura et al., 2001; Lee et al., 1988; Anderson et al., 2017). In ref. 14 we investigated the accuracy of the aug-pc-$\infty$ Gaussian basis set that was used in ref. 15. The study was restricted to ${}^{n}S$ states to avoid symmetry breaking effects, which were still observed for \ceH-, He, \ceLi-, and \ceNa-, as was discussed in the Introduction. Reproducing symmetry preserving data with Erkale (Lehtola et al., 2012), we found that the truncation error of the aug-pc-$\infty$ basis set is less than 1 $\mu E_{h}$ for light atoms and tens of $\mu E_{h}$ for heavier atoms in Hartree–Fock and BHHLYP (Becke, 1993) calculations.

Because the screening is evaluated analytically in Gaussian-basis calculations (Heyd et al., 2003) and the accuracy of the aug-pc-$\infty$ basis set has been established (Lehtola, 2019b), the values reported in ref. 15 offer an ideal reference for the present work. The comparison of results obtained in the present work with \eqrangerefphinDnk and a numerical basis set with five 15-node radial elements and a practical infinity $r_{\infty}=40a_{0}$ is shown in 3, demonstrating excellent agreement between the fully numerical and Gaussian basis calculations. The values are in full agreement after rounding to the same accuracy for the light atoms, while the fully numerical values are slightly below the Gaussian-basis values for the heavier atoms, as expected based on the basis set truncation errors observed in ref. 14.

Finally, the spin-restricted ground states for all atoms in the periodic table at HF and HFS levels of theory are shown in LABEL:tab:HF-gs,_HFS-gs, respectively; these calculations also used ten 15-node radial elements. The data reveal that in some cases a lower-lying configuration has been seen in the brute force search, but that it failed to converge. In the HF calculations, the $6s^{2}4f^{2}5d^{1}$ state of Pr converges without problems; however, the $6s^{2}4f^{3}$ configuration has a lower energy but its wave function failed to converge. Similar issues are observed in the HFS calculations for Cf, Es, and Fm, where the $5f^{n-1}7s^{1}$ state converges without problems, but a lower energy is observed for a $5f^{n}$ configuration the wave function of which fails to converge.

The HF results can be compared to the high-accuracy data for multiconfigurational HF of Saito (Saito, 2009). Because the present calculations are spin-restricted with fractional occupations, the energies are higher than those reported in ref. 132. However, the agreement for the noble gases is perfect, underlining the high accuracy of the computational approach used in the present work, which was outlined in ref. 14, even though only 139 radial basis functions were employed.

IV Summary and discussion

We have described new efficient implementations of range-separated functionals as well as fractional occupations for atomic electronic structure calculations with HelFEM, and demonstrated that beyond microhartree accucacy can be achieved with just 139 numerical radial basis functions. We have tested the program by reproducing local density approximation (LDA) and generalized gradient approximation (GGA) total energies for $1\leq Z\leq 86$ at the basis set limit, and shown that the literature values deviate from the complete basis set limit by up to $10\ \mu E_{h}$. The approaches developed in the present work could be straightforwardly extended to fractional occupations per shell in future work, requiring the addition of a logic to formulate the fractional occupation numbers.

The capabilities added to HelFEM in the present work allow for self-consistent benchmarking of density functionals at the basis set limit, which is useful for development and implementation purposes. For instance, Clementi–Roetti wave functions (Clementi and Roetti, 1974) are often used for non-self-consistent benchmarks of density functionals, but the availability of a program for self-consistent calculations is certain to help future developments as numerical instabilities in the functional may not be detected in non-self-consistent calculations.

Furthermore, we have reported the non-relativistic spin-restricted ground state configurations of all atoms in the periodic table at HF and HFS levels of theory. Such knowledge is useful for implementations of the superposition of atomic densities guess (Almlöf et al., 1982; Van Lenthe et al., 2006), which is often implemented based on spin-restricted fractionally occupied calculations. The present approach is also useful for implementations of the SAP guess (Lehtola, 2019c). For instance, the implementation of SAP now available in the development version of the Psi4 program (Parrish et al., 2017) is based on HFS potentials tabulated during the present work. Instead of the 4000 point tabulation used in ref. 78 with unknown error, the ten-element calculations of the present work yield 751-point tabulations that reproduce the sub-microhartree-level accuracy of the original calculation.

The atomic orbitals obtained from the present approach may also be useful for initializing fully numerical molecular electronic structure calculations via either a superposition of atomic densities, or in combination to the extended Hückel rule developed in ref. 78.

Funding information

This work has been supported by the Academy of Finland through project number 311149.

Acknowledgments

I thank Dirk Andrae, Volker Blum, Eli Kraisler, Jacek Kobus, Leeor Kronik, Micael Oliveira, Dage Sundholm, Edward Valeev, and Lucas Visscher for discussions. Computational resources provided by CSC – It Center for Science Ltd (Espoo, Finland) and the Finnish Grid and Cloud Infrastructure (persistent identifier urn:nbn:fi:research-infras-2016072533) are gratefully acknowledged.