Initialization

Needs["MaTeX`"]
SetOptions[MaTeX,"Preamble"->{"\\usepackage{amssymb,amsmath,latexsym,amsfonts,amsthm,mathpazo,xcolor,bm,mhchem}"}];

HaToeV=27.21138602;

DFAs

LDA exchange

Dirac

{α=-(3/2)((((3/(4π))))^(1/3)),-(3/2)((((3/(4π))))^(1/3))//N}

ϵ[ρ_]=α (ρ^(1/3));	v[ρ_]=α (ρ^(4/3));

Simplify[(∂_ρ)v[ρ]==(4/3) α (ρ^(1/3))]

LDA correlation

VWN5

(Subscript[Υ, x])[ζ_]=(((((1+ζ))^(4/3))+(((1-ζ))^(4/3)))/2);	(Subscript[Υ, c])[ζ_]=(((Subscript[Υ, x])[ζ]-1)/((Subscript[Υ, x])[1]-1));	(((Subscript[Υ, x])[ζ]-1)/((Subscript[Υ, x])[1]-1))==((1/2)(((((1+ζ))^(4/3))+(((1-ζ))^(4/3))-2)/((2^(1/3))-1)))//Simplify
(Subscript[e, c]^VWN)[(Subscript[r, s]),ζ]=(Subscript[e, c]^VWN)[(Subscript[r, s]),0]+(Subscript[e, a])[(Subscript[r, s])](((Subscript[Υ, c])[ζ])/((Subscript[Υ, c]^″)[0]))(1-(ζ^4))+((Subscript[e, c]^VWN)[(Subscript[r, s]),1]-(Subscript[e, c]^VWN)[(Subscript[r, s]),0])(Subscript[Υ, c])[ζ](ζ^4);
(∂_(ζ,ζ))(Subscript[Υ, c])[ζ]==(4/(9 ((2^(1/3))-1)))/.ζ->0

D[(Subscript[Υ, c])[ζ],ζ]==(4/3)(((((1+ζ))^(1/3))-(((1-ζ))^(1/3)))/(2((2^(1/3))-1)))//Simplify

D[(Subscript[Υ, c])[ζ],{ζ,2}]==(1/((Subscript[Υ, x])[1]-1))D[(Subscript[Υ, x])[ζ],{ζ,2}]//Simplify
D[(Subscript[Υ, x])[ζ],{ζ,2}]==(4/9)(((((1+ζ))^(-2/3))+(((1-ζ))^(-2/3)))/2)//Simplify

D[(Subscript[e, a])[x],x]==A (+ (2/x)-((4 b)/( (((b+2 x))^2)+((Q[b,c])^2)))-(( Derivative[1,0,0][X][x,b,c])/(X[x,b,c]))-((b x0)/(X[x0,b,c])) ( (2/(x-x0))-((4 (b+2 x0))/( (((b+2 x))^2)+((Q[b,c])^2))) -(( Derivative[1,0,0][X][x,b,c])/(X[x,b,c]))))//Simplify

(Subscript[e, a])[x_]=A (Log[((x^2)/(X[x,b,c]))]+((2b)/(Q[b,c])) ArcTan[((Q[b,c])/(2x+b))]-((b x0)/(X[x0,b,c]))(Log[((((x-x0))^2)/(X[x,b,c]))]+2(((b+2x0))/(Q[b,c]))ArcTan[((Q[b,c])/(2x+b))]));

eDFAs

eLDA exchange

Ground state exchange functional for glomium

The ground state two-electron glomium system has an uniform electron density

ρ==(2/(2(π^2)(R^3)))	(Subscript[ρ, α])==(Subscript[ρ, β])==(1/(2(π^2)(R^3)))

The reduced HF energy of glomium is

(Subscript[e, HF]^((0)))==(4/(3 π R))

This energy can be decomposed as

(Subscript[t, s]^((0)))==0	(Subscript[e, H]^((0)))==(8/(3 π R))	(Subscript[e, x]^((0)))==-(4/(3 π R))

Knowing that the exchange functional has the following form

(Subscript[Ε, x,α]^((0)))==∫(Subscript[e, x]^((0)))(Subscript[ρ, α]) \[DifferentialD]r==(Subscript[C, x])∫(Subscript[ρ, α]^(4/3)) \[DifferentialD]r

This yields

(Subscript[Ε, x]^((0)))==-(8/(3π R))==(Subscript[C, x]^((0)))((((1/(2(π^2)(R^3)))))^(4/3))2(π^2)(R^3)	\[Implies]	(Subscript[C, x]^((0)))==-(4/3) ((((2/π)))^(1/3))

The ground-state exchange functional for glomium is

(Subscript[e, x]^((0)))[ρ]==(Subscript[C, x]^((0)))(ρ^(1/3))	(Subscript[C, x]^((0)))==-(4/3) ((((2/π)))^(1/3))

Doubly-excited state exchange functional for glomium

Any doubly-excited state two-electron glomium system has the same uniform electron density as the ground state

ρ==(2/(2(π^2)(R^3)))	(Subscript[ρ, α])==(Subscript[ρ, β])==(1/(2(π^2)(R^3)))

The reduced HF energy of the first doubly-excited state of glomium is

(Subscript[e, HF]^((1)))==(3/(2(R^2)))+(176/(105 π R))

(3/(2(R^2)))+(176/(105 π R))/.R->(1/((π^(2/3)) (ρ^(1/3))))

This energy can be decomposed as

(Subscript[t, s]^((1)))==(3/(2(R^2)))	(Subscript[e, H]^((1)))==(352/(105 π R))	(Subscript[e, x]^((1)))==-(176/(105 π R))

Knowing that the exchange functional has the following form

(Subscript[Ε, x,α]^((1)))==∫(Subscript[e, x]^((1)))(Subscript[ρ, α]) \[DifferentialD]r==(Subscript[C, x]^((1)))∫(Subscript[ρ, α]^(4/3)) \[DifferentialD]r

This yields

(Subscript[Ε, x]^((1)))==-(176/(105 π R))==(Subscript[C, x]^((1)))((((1/(2(π^2)(R^3)))))^(4/3))2(π^2)(R^3)	\[Implies]	(Subscript[C, x]^((1)))==-(176/105) ((((2/π)))^(1/3))

The double-excited state exchange functional for glomium is

(Subscript[e, x]^((1)))[ρ]==(Subscript[C, x]^((1)))(ρ^(1/3))	(Subscript[C, x]^((1)))==-(176/105) ((((2/π)))^(1/3))

Weight-dependent exchange functional for glomium

We can now combine these two exchange functionals to create a weight-dependent exchange functional

(Subscript[e, x]^w)[ρ]==(1-w)(Subscript[e, x]^((0)))[ρ]+w (Subscript[e, x]^((1)))[ρ]

(Subscript[e, x]^w)[ρ]==(1-w)(Subscript[C, x]^((0)))(ρ^(1/3))+w (Subscript[C, x]^((1)))[ρ](ρ^(1/3))
           == (Subscript[C, x]^((w)))(ρ^(1/3))

The weight dependent exchange functional for glomium is then

(Subscript[e, x]^w)[ρ]==(Subscript[C, x]^w)(ρ^(1/3))	(Subscript[C, x]^w)==(1-w)(Subscript[C, x]^((0)))+w (Subscript[C, x]^((1)))	(Subscript[C, x]^((0)))==-(4/3) ((((2/π)))^(1/3))	(Subscript[C, x]^((1)))==-(176/105) ((((2/π)))^(1/3))

Amazingly, the weight dependence of the exchange functional can be transfered to the (Subscript[C, x]) coefficient.
This is obvious but kind of nice.

LDA-centered weight-dependent exchange functional for molecules

In order to create a more “universal” functional (that does not depend on the number of electrons in particular), we can “shift” the weight-dependent exchange functional designed for glomium to make it LDA centered.
The corresponding definition of the so-called LDA-centered weight-dependent exchange functional is then

(Subscript[(OverBar[e]), x]^w)[ρ]==(1-w)((Subscript[e, x]^((0)))[ρ]+((Subscript[e, x]^LDA)[ρ]-(Subscript[e, x]^((0)))[ρ]))+w ((Subscript[e, x]^((1)))[ρ]+((Subscript[e, x]^LDA)[ρ]-(Subscript[e, x]^((0)))[ρ]))
	  ==(1-w)((Subscript[C, x]^((0)))+(Subscript[C, x]^LDA)-(Subscript[C, x]^((0))))(ρ^(1/3))+w ((Subscript[C, x]^((1)))+(Subscript[C, x]^LDA)-(Subscript[C, x]^((0))))(ρ^(1/3))
	  ==((Subscript[C, x]^LDA)+w ((Subscript[C, x]^((1)))-(Subscript[C, x]^((0)))))(ρ^(1/3))
	  ==(Subscript[(OverBar[C]), x]^w)(ρ^(1/3))

where we recall that the usual LDA exchange functional (D30) is given by

(Subscript[e, x]^LDA)[ρ]==(Subscript[C, x]^LDA)(ρ^(1/3))	(Subscript[C, x]^LDA)==-(3/2)((((3/(4π))))^(1/3))

The LDA-centered weight dependent exchange functional is then

(Subscript[(OverBar[e]), x]^w)[ρ]==(Subscript[(OverBar[C]), x]^w)(ρ^(1/3))	(Subscript[(OverBar[C]), x]^w)==(Subscript[C, x]^LDA)+w ((Subscript[C, x]^((1)))-(Subscript[C, x]^((0))))	(Subscript[C, x]^((0)))==-(4/3) ((((2/π)))^(1/3))	(Subscript[C, x]^((1)))==-(176/105) ((((2/π)))^(1/3))	(Subscript[C, x]^LDA)==-(3/2)((((3/(4π))))^(1/3))

This is the final form of our functional. The Cx coefficient is weight dependent and it appears as a correction over Dirac’s Cx coefficient.

eLDA correlation

Ground state correlation functional for glomium

The ground-state correlation energy of glomium can be very accurately computed with Hylleraas-type calculations.
We use a Pade-type fit to obtain the functional which reads

(Subscript[e, c]^((0)))[ρ]==((a^((0)))/(1+(b^((0)))(ρ^(-1/6))+(c^((0)))(ρ^(-1/3))))	(a^((0)))==-0.0238184	(b^((0)))==0.00575719	(c^((0)))==0.0830576

The first coefficient is actually not a fitting coefficient but is obtained via the exact value of the correlation energy is the high-density limit.

Doubly-excited state correlation functional for glomium

Similarly, the correlation functional for the first doubly-excited states is

(Subscript[e, c]^((1)))[ρ]==((a^((1)))/(1+(b^((1)))(ρ^(-1/6))+(c^((1)))(ρ^(-1/3))))	(a^((1)))==-0.0144633	(b^((1)))==-0.0504501	(c^((1)))==0.0331287

Weight-dependent exchange functional for glomium

We can now combine these two correlations functionals to create a weight-dependent correlation functional

(Subscript[e, c]^w)[ρ]==(1-w)(Subscript[e, c]^((0)))[ρ]+w (Subscript[e, c]^((1)))[ρ]

LDA-centered weight-dependent correlation functional for molecules

The corresponding LDA-centered weight-dependent correlation functional is then

(Subscript[(OverBar[e]), c]^w)[ρ]==(1-w)((Subscript[e, c]^((0)))[ρ]+((Subscript[e, c]^LDA)[ρ]-(Subscript[e, c]^((0)))[ρ]))+w ((Subscript[e, c]^((1)))[ρ]+((Subscript[e, c]^LDA)[ρ]-(Subscript[e, c]^((0)))[ρ]))

where we use the VWN5 functional as the LDA correlation functional (cf DFAs section)

(Subscript[e, c]^LDA)[ρ]≡(Subscript[e, c]^VWN5)[ρ]

Note that this functional is only valid for closed-shell system.
If one wants to extend this to open-shell system, one needs to take into account spin polarization.
There are several ways of doing this...
The extension to GGAs is also more tricky as the GGA is not usual a multiplicative factor compared to the LDA functional.

With[{M=10},
(Subscript[ψ, n_])[u_]=(u^(n-1));
H[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u]((((u^2)/(4(R^2)))-1)(Subscript[ψ, j]^″)[u]+(((5u)/(4(R^2)))-(2/u))  (Subscript[ψ, j]^′)[u]+(1/u)(Subscript[ψ, j])[u])(u^2) (Sqrt[1-((u^2)/(4 (R^2)))])\[DifferentialD]u,{i,M},{j,M}];
S[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u](Subscript[ψ, j])[u](u^2) (Sqrt[1-((u^2)/(4 (R^2)))])\[DifferentialD]u,{i,M},{j,M}];
];

Eig0=Table[{R,Sort[Eigenvalues[{N[H[R],100],N[S[R],100]}]][[1]]-((8/(3 π R)))},{R,{(1/10),(1/5),(1/2),1,2,5,10,20,50,100,200}}];
Eig2=Table[{R,Sort[Eigenvalues[{N[H[R],100],N[S[R],100]}]][[2]]-((3/(R^2))+(352/(105 π R)))},{R,{(1/10),(1/5),(1/2),1,2,5,10,20,50,100,200}}];

Table[{R,Sort[Eigenvalues[{N[H[R],100],N[S[R],100]}]][[1]]-((8/(3 π R)))},{R,Table[(10^(-n)),{n,8,4,-1}]}]
Table[{R,Sort[Eigenvalues[{N[H[R],100],N[S[R],100]}]][[2]]-((3/(R^2))+(352/(105 π R)))},{R,Table[(10^(-n)),{n,8,4,-1}]}]

FindFit[Eig0,(a/(1+b (R^(1/2))+c R))/.a->-0.04763687334741862`,{b,c},R,MaxIterations->10000]
(*FindFit[Eig1,(a/(1+b (R^(1/2))+c R))/.a->-0.05656280,{b,c},R,MaxIterations->10000]*)
FindFit[Eig2,(a/(1+b (R^(1/2))+c R))/.a->-0.02892648,{b,c},R,MaxIterations->10000]

ϵc0[ρ_]:=Module[{R=(1/((π^(2/3)) (ρ^(1/3)))),a=-0.0238184,b=0.00843194,c=0.178161},(a/(1+b (R^(1/2))+c R))]
ϵc2[ρ_]:=Module[{R=(1/((π^(2/3)) (ρ^(1/3)))),a=-0.01446325,b=-0.0738888,c=0.0710620},(a/(1+b (R^(1/2))+c R))]

Show[{Plot[{2ϵc0[(2/(2(π^2)(R^3)))],2ϵc1[(2/(2(π^2)(R^3)))],2ϵc2[(2/(2(π^2)(R^3)))]},{R,0,200},PlotRange->All],ListPlot[{Eig0,Eig1,Eig2}]}]

Uniform scaling

The (Subscript[C, x])coefficient which gives the exact exchange energy for the hydrogen atom is

-((20 )/27)((((π/2)))^(1/3))//N

((-((20 )/27)((((π/2)))^(1/3)))/(-(3/4)((((3/π)))^(1/3))))//N

The ratio of this coefficient with the LDA (Subscript[C, x]) coefficient does not violate the Lieb-Oxford bound (which is equal to 1.174).

The (Subscript[C, x]^((0))) coefficient does violate the Lieb-Oxford bound by a significant amount

((-(4/3)((((1/π)))^(1/3)))/(-(3/4)((((3/π)))^(1/3))))//N

The (Subscript[C, x]^((0))) and (Subscript[C, x]^((1))) coefficients can be scaled in such a way that it doesn’t does not violate the Lieb-Oxford bound

Solve[α((-(4/3)((((1/π)))^(1/3)))/(-(3/4)((((3/π)))^(1/3))))==1.174]

We can also scale the (Subscript[C, x]^((0))) and (Subscript[C, x]^((1))) coefficients in such a way that it does reproduce the hydrogen atom value

Solve[α((-(4/3)((((1/π)))^(1/3)))/(-(3/4)((((3/π)))^(1/3))))==((-((20 )/27)((((π/2)))^(1/3)))/(-(3/4)((((3/π)))^(1/3))))]