Initialization

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

HaToeV=27.21138602;

(Subscript[δ, a_,b_]):=KroneckerDelta[a,b];

(Subscript[F, 0])[t_]:=If[t==0,1,(1/2)(Sqrt[(π/t)])Erf[(Sqrt[t])]];

Γ[x_]:=Gamma[x];

$Assumptions->R>0;

Γ[x_]=Gamma[x];

Γ[x_,t_]=Gamma[x,t];

(Subscript[δ, a_,b_])=KroneckerDelta[a,b];

(Subscript[δ, x_])=KroneckerDelta[x];

(Subscript[H, n_])[x_]=HermiteH[n,x];

SortEigensystem[eigsys_]:=Transpose[Sort[Transpose[eigsys]]]

SortEigensystemGreater[eigsys_]:=Transpose[Sort[Transpose[eigsys],Greater]]

TakeMOs[c_,n_]:=Take[c\[Transpose],n]\[Transpose];

Needs["MaTeX`"]

Needs["NumericalDifferentialEquationAnalysis`"]

Id=(({{1, 0, 0}, {0, 1, 0}, {0, 0, 1}}));

AngQ[Shell_]:=Which[
Shell=="S",{{0,0,0}}
,
Shell=="P",{{1,0,0},{0,1,0},{0,0,1}}
,
Shell=="SP",{{0,0,0},{1,0,0},{0,1,0},{0,0,1}}
,
Shell=="D",{{2,0,0},{0,2,0},{0,0,2},{1,1,0},{1,0,1},{0,1,1}}
,
Shell=="F",{{3,0,0},{0,3,0},{0,0,3},{2,1,0},{2,0,1},{1,2,0},{0,2,1},{1,0,2},{0,1,2},{1,1,1}}
];

$Assumptions=α>0&&r>0&&x∈Reals&&y∈Reals&&z∈Reals&&Z>0&&u>0&&AB>0;

Cofactor[m_List?MatrixQ,{i_Integer,j_Integer}]:=If[Length[m]==1,1,(-1)^(i+j) Det[Drop[Transpose[Drop[Transpose[m],{j}]],{i}]]]

MinorMatrix[m_List?MatrixQ]:=Map[Reverse,Minors[m],{0,1}]

CofactorMatrix[m_List?MatrixQ]:=MapIndexed[#1 (-1)^(Plus@@#2)&,MinorMatrix[m],{2}]

NormCoeff[α_,{ax_,ay_,az_}]:=(((((2α)/π)))^(3/4))((((4α))^((ax+ay+az)/2))/((((2ax-1)!!(2ay-1)!!(2az-1)!!))^(1/2)))

(Subscript[F, m_])[t_]:=If[t==0,(1/(2m+1)),((Γ[m+(1/2)]-Γ[m+(1/2),t])/(2(t^(m+1/2))))];

Matrix elements

𝒮[{α_,A_},{β_,B_}]:=(((((2α)/π)))^(3/4))(((((2β)/π)))^(3/4))((((π/(α+β))))^(3/2))Exp[-((α β)/(α+β)) ((Norm[A-B])^2)]

𝒯[{α_,A_},{β_,B_}]:=((α β)/(α+β))(3-((2α β)/(α+β))((Norm[A-B])^2))𝒮[{α,A},{β,B}]

𝒱[{α_,A_},{β_,B_},{ZC_,C_}]:=-((2π ZC)/(α+β))(((((2α)/π)))^(3/4))(((((2β)/π)))^(3/4))(Subscript[F, 0])[(α+β)((Norm[((α A+β B)/(α+β))-C])^2)]Exp[-((α β)/(α+β)) ((Norm[A-B])^2)]

𝒰[{α_,A_},{β_,B_},{γ_,C_},{δ_,D_}]:=((2(π^(5/2)) )/((α+β)(γ+δ)(Sqrt[α+β+γ+δ])))(((((2α)/π)))^(3/4))(((((2β)/π)))^(3/4))(((((2γ)/π)))^(3/4))(((((2δ)/π)))^(3/4))(Subscript[F, 0])[(((α+β)(γ+δ))/(α+β+γ+δ))((Norm[((α A+β B)/(α+β))-((γ C+δ D)/(γ+δ))])^2)]Exp[-((α β)/(α+β)) ((Norm[A-B])^2)]Exp[-((γ δ)/(γ+δ)) ((Norm[C-D])^2)]

DFAs

eDFAs

eLDA exchange

eLDA correlation

Explicitly-correlated calculations

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}];
];

With[{M=10},
(Subscript[ψ, n_])[u_]=(u^(n-1));
H1[R_]=Table[(∫_0^(2R))(Subscript[ψ, i])[u]((((u^2)/(4(R^2)))-1)(Subscript[ψ, j]^″)[u]+(((7u)/(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}];
];

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

FindFit[Eig0,(a/(1+b (R^(1/2))+c R))/.a->-0.0238184,{b,c},R,MaxIterations->10000]
FindFit[Eig1,(a/(1+b (R^(1/2))+c R))/.a->-0.0282814,{b,c},R,MaxIterations->10000]
FindFit[Eig2,(a/(1+b (R^(1/2))+c R))/.a->-0.0144633,{b,c},R,MaxIterations->10000]

ϵc0[ρ_]:=Module[{R=(1/((π^(2/3)) (ρ^(1/3)))),a=-0.0238184,b=0.007923352705637004`,c=0.17820175850290984`},(a/(1+b (R^(1/2))+c R))]
ϵc1[ρ_]:=Module[{R=(1/((π^(2/3)) (ρ^(1/3)))),a=-0.0282814,b=0.004011883522661856`,c=0.14262596914269524`},(a/(1+b (R^(1/2))+c R))]
ϵc2[ρ_]:=Module[{R=(1/((π^(2/3)) (ρ^(1/3)))),a=-0.0144633,b=-0.07411135123824349`,c=0.07108999554813503`},(a/(1+b (R^(1/2))+c R))]

SizeLegend=1.5;
SizeEpilog=3;
SizeLabel=2;
OptionsEc={
AspectRatio->1,
PlotMarkers->{"OpenMarkers",Large},
Frame->True,
FrameStyle->Directive[Thick,20],
ImageSize->Large,
PlotTheme->"Scientific",
PlotRange->All,FrameLabel->{
MaTeX["\\text{Radius of the glome $R$ (a.u.)}",Magnification->SizeLabel],
MaTeX["\\text{Reduced correlation energy $\\epsilon_c^{(I)}$ (a.u.)}",Magnification->SizeLabel]
}
};
Show[{ListPlot[{
Table[{Eig0[[k,1]],Eig0[[k,2]]},{k,Length[Eig0]}],
(*Table[{Eig1[[k,1]],Eig1[[k,2]]},{k,Length[Eig1]}],*)
Table[{Eig2[[k,1]],Eig2[[k,2]]},{k,Length[Eig2]}]
},OptionsEc,PlotLegends->Placed[{
MaTeX["\\text{Ground state $I=0$}",Magnification->SizeLabel],
(*MaTeX["\\text{Singly-excited state $I=1$}",Magnification->SizeLabel],*)
MaTeX["\\text{Doubly-excited state $I=1$}",Magnification->SizeLabel]},{Right,Bottom}]],Plot[{ϵc0[(1/((π^2) (R^3)))],(*ϵc1[(1/((π^2) (R^3)))],*)ϵc2[(1/((π^2) (R^3)))]},{R,0,150},PlotRange->All,PlotStyle->Thick,PlotTheme->"Scientific"]}]

Ground state correlation functional for glomium

Singly-excited state correlation functional for glomium

Doubly-excited state correlation functional for glomium

Weight-dependent exchange functional for glomium

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.

Uniform scaling

Hydrogen molecule

Molecular coordinates

ZA=1;	ZB=1;	rA[R_]={0,0,-(R/2)};	rB[R_]={0,0,+(R/2)};

Basis functions

With[{nG=3,d={0.1543289673,0.5353281423,0.4446345422},α={3.425250914,0.6239137298,0.1688554040}},
(Subscript[χ, A])[{x_,y_,z_},R_]=(∑_(k=1)^nG)d[[k]](((((2α[[k]])/π)))^(3/4))Exp[-α[[k]]((Norm[{x,y,z}-rA[R]])^2)];
(Subscript[χ, B])[{x_,y_,z_},R_]=(∑_(k=1)^nG)d[[k]](((((2α[[k]])/π)))^(3/4))Exp[-α[[k]]((Norm[{x,y,z}-rB[R]])^2)]
];

Integrals over basis functions

(*With[{nG=3,d={0.1543289673,0.5353281423,0.4446345422},α={3.425250914,0.6239137298,0.1688554040}},
(Subscript[S, AA])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]𝒮[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]}]];	(Subscript[S, BB])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]𝒮[{α[[(Subscript[k, 1])]],rB[R]},{α[[(Subscript[k, 2])]],rB[R]}]];	(Subscript[S, AB])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]𝒮[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rB[R]}]];
(Subscript[t, AA])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]𝒯[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]}]];	(Subscript[t, BB])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]𝒯[{α[[(Subscript[k, 1])]],rB[R]},{α[[(Subscript[k, 2])]],rB[R]}]];	(Subscript[t, AB])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]𝒯[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rB[R]}]];
(Subscript[v, AA])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]](𝒱[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]},{ZA,rA[R]}]+𝒱[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]},{ZB,rB[R]}])];(Subscript[v, BB])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]](𝒱[{α[[(Subscript[k, 1])]],rB[R]},{α[[(Subscript[k, 2])]],rB[R]},{ZA,rA[R]}]+𝒱[{α[[(Subscript[k, 1])]],rB[R]},{α[[(Subscript[k, 2])]],rB[R]},{ZB,rB[R]}])];(Subscript[v, AB])[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]](𝒱[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rB[R]},{ZA,rA[R]}]+𝒱[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rB[R]},{ZB,rB[R]}])];
(Subscript[h, AA])[R_]=(Subscript[t, AA])[R]+(Subscript[v, AA])[R];	(Subscript[h, BB])[R_]=(Subscript[t, BB])[R]+(Subscript[v, BB])[R];	(Subscript[h, AB])[R_]=(Subscript[t, AB])[R]+(Subscript[v, AB])[R];
AAAA[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)(∑_((Subscript[k, 3])=1)^nG)(∑_((Subscript[k, 4])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]d[[(Subscript[k, 3])]]d[[(Subscript[k, 4])]]𝒰[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]},{α[[(Subscript[k, 3])]],rA[R]},{α[[(Subscript[k, 4])]],rA[R]}]];
BBBB[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)(∑_((Subscript[k, 3])=1)^nG)(∑_((Subscript[k, 4])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]d[[(Subscript[k, 3])]]d[[(Subscript[k, 4])]]𝒰[{α[[(Subscript[k, 1])]],rB[R]},{α[[(Subscript[k, 2])]],rB[R]},{α[[(Subscript[k, 3])]],rB[R]},{α[[(Subscript[k, 4])]],rB[R]}]];
ABAB[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)(∑_((Subscript[k, 3])=1)^nG)(∑_((Subscript[k, 4])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]d[[(Subscript[k, 3])]]d[[(Subscript[k, 4])]]𝒰[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rB[R]},{α[[(Subscript[k, 3])]],rA[R]},{α[[(Subscript[k, 4])]],rB[R]}]];
AABB[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)(∑_((Subscript[k, 3])=1)^nG)(∑_((Subscript[k, 4])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]d[[(Subscript[k, 3])]]d[[(Subscript[k, 4])]]𝒰[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]},{α[[(Subscript[k, 3])]],rB[R]},{α[[(Subscript[k, 4])]],rB[R]}]];
AAAB[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)(∑_((Subscript[k, 3])=1)^nG)(∑_((Subscript[k, 4])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]d[[(Subscript[k, 3])]]d[[(Subscript[k, 4])]]𝒰[{α[[(Subscript[k, 1])]],rA[R]},{α[[(Subscript[k, 2])]],rA[R]},{α[[(Subscript[k, 3])]],rA[R]},{α[[(Subscript[k, 4])]],rB[R]}]];
BBAB[R_]=Simplify[(∑_((Subscript[k, 1])=1)^nG)(∑_((Subscript[k, 2])=1)^nG)(∑_((Subscript[k, 3])=1)^nG)(∑_((Subscript[k, 4])=1)^nG)d[[(Subscript[k, 1])]]d[[(Subscript[k, 2])]]d[[(Subscript[k, 3])]]d[[(Subscript[k, 4])]]𝒰[{α[[(Subscript[k, 1])]],rB[R]},{α[[(Subscript[k, 2])]],rB[R]},{α[[(Subscript[k, 3])]],rA[R]},{α[[(Subscript[k, 4])]],rB[R]}]];
]*)

Spatial orbitals in minimal basis:

(*(Subscript[ϕ, 1])[{x_,y_,z_},R_]=(((Subscript[χ, A])[{x,y,z},R]+(Subscript[χ, B])[{x,y,z},R])/(Sqrt[2(1+(Subscript[S, AB])[R])]));	(Subscript[ϕ, 2])[{x_,y_,z_},R_]=(((Subscript[χ, A])[{x,y,z},R]-(Subscript[χ, B])[{x,y,z},R])/(Sqrt[2(1-(Subscript[S, AB])[R])]));*)

One-electron densities

(*(Subscript[n, 0])[{x_,y_,z_},R_]=2(((Subscript[ϕ, 1])[{x,y,z},R])^2);	(Subscript[n, 1])[{x_,y_,z_},R_]=(((Subscript[ϕ, 1])[{x,y,z},R])^2)+(((Subscript[ϕ, 2])[{x,y,z},R])^2);	(Subscript[n, 2])[{x_,y_,z_},R_]=2(((Subscript[ϕ, 2])[{x,y,z},R])^2);*)

(*nw[{x_,y_,z_},R_,w_]=(1-w[[1]]-w[[2]])(Subscript[n, 0])[{x,y,z},R]+w[[1]] (Subscript[n, 1])[{x,y,z},R]+w[[2]] (Subscript[n, 2])[{x,y,z},R];*)

Integrals over orbitals

(*(Subscript[h, 1])[R_]=Simplify[(1/(2(1+(Subscript[S, AB])[R])))((Subscript[h, AA])[R]+(Subscript[h, BB])[R]+2(Subscript[h, AB])[R])];	(Subscript[h, 2])[R_]=Simplify[(1/(2(1-(Subscript[S, AB])[R])))((Subscript[h, AA])[R]+(Subscript[h, BB])[R]-2(Subscript[h, AB])[R])];*)

(*(Subscript[J, 11])[R_]=Simplify[((((1/(2(1+(Subscript[S, AB])[R])))))^2)(AAAA[R]+BBBB[R]+2AABB[R]+4ABAB[R]+4AAAB[R]+4BBAB[R])];
(Subscript[J, 22])[R_]=Simplify[((((1/(2(1-(Subscript[S, AB])[R])))))^2)(AAAA[R]+BBBB[R]+2AABB[R]+4ABAB[R]-4AAAB[R]-4BBAB[R])];
(Subscript[J, 12])[R_]=Simplify[((1/(2(1+(Subscript[S, AB])[R]))))((1/(2(1-(Subscript[S, AB])[R]))))(AAAA[R]+BBBB[R]+2AABB[R]-4ABAB[R])];*)

(*(Subscript[K, 11])[R_]=(Subscript[J, 11])[R];
(Subscript[K, 22])[R_]=(Subscript[J, 22])[R];
(Subscript[K, 12])[R_]=Simplify[((1/(2(1+(Subscript[S, AB])[R]))))((1/(2(1-(Subscript[S, AB])[R]))))(AAAA[R]+BBBB[R]-2AABB[R])];*)

Nuclear energy

(*(Subscript[Ε, NN])[R_]=((ZA ZB)/(Norm[rA[R]-rB[R]]));*)

Hartree-Fock energies

(Subscript[Ε, HF])[R_]={
2(Subscript[h, 1])[R]+2(Subscript[J, 11])[R]-(Subscript[K, 11])[R],
(Subscript[h, 1])[R]+(Subscript[h, 2])[R]+(1/2)((Subscript[J, 11])[R]+(Subscript[J, 22])[R]+2(Subscript[J, 12])[R])-(1/2)((Subscript[K, 11])[R]+(Subscript[K, 22])[R]),
2(Subscript[h, 2])[R]+2(Subscript[J, 22])[R]-(Subscript[K, 22])[R]
}+(Subscript[Ε, NN])[R];

(*(Subscript[Εw, HF])[R_,w_]=(1-w[[1]]-w[[2]])2(Subscript[h, 1])[R]+w[[1]]((Subscript[h, 1])[R]+(Subscript[h, 2])[R])+w[[2]]2(Subscript[h, 2])[R]+(((1-w[[1]]-w[[2]]))^2)(2(Subscript[J, 11])[R]-(Subscript[K, 11])[R])+((w[[1]])^2)((Subscript[J, 12])[R]+(Subscript[K, 12])[R])+((w[[2]])^2)(2(Subscript[J, 22])[R]-(Subscript[K, 22])[R])+2(1-w[[1]]-w[[2]])w[[2]](2(Subscript[J, 12])[R]-(Subscript[K, 12])[R])+(Subscript[Ε, NN])[R];*)

(*(Subscript[ΔΕw, HF])[R_,w_]=(1-w[[1]]-w[[2]])(Subscript[Ε, HF])[R][[1]]+w[[1]](Subscript[Ε, HF])[R][[2]]+w[[2]](Subscript[Ε, HF])[R][[3]]-(Subscript[Εw, HF])[R,w];*)

(Subscript[ϵw, HF])[R_,w_]={(Subscript[h, 1])[R]+(1-w)(2(Subscript[J, 11])[R]-(Subscript[K, 11])[R])+w(2(Subscript[J, 12])[R]-(Subscript[K, 12])[R]),(Subscript[h, 2])[R]+(1-w)(2(Subscript[J, 12])[R]-(Subscript[K, 12])[R])+w(2(Subscript[J, 22])[R]-(Subscript[K, 22])[R])};

LDA energies

(*(Subscript[Ε, LDA])[R_]:={
2(Subscript[h, 1])[R]+2(Subscript[J, 11])[R]+NIntegrate[ϵxcLDA[(Subscript[n, 0])[{x,y,z},R]](Subscript[n, 0])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}],
(Subscript[h, 1])[R]+(Subscript[h, 2])[R]+(1/2)((Subscript[J, 11])[R]-NIntegrate[ϵxcLDA[(Subscript[n, 0])[{x,y,z},R]](Subscript[n, 0])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}])+(1/2)((Subscript[J, 22])[R]-NIntegrate[ϵxcLDA[(Subscript[n, 2])[{x,y,z},R]](Subscript[n, 2])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}])+(Subscript[J, 12])[R],
2(Subscript[h, 2])[R]+2(Subscript[J, 22])[R]+NIntegrate[ϵxcLDA[(Subscript[n, 2])[{x,y,z},R]](Subscript[n, 2])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]
}+(Subscript[Ε, NN])[R];*)

(*(Subscript[Εw, LDA])[R_,w_]:=(1-w[[1]]-w[[2]])2(Subscript[h, 1])[R]+w[[1]]((Subscript[h, 1])[R]+(Subscript[h, 2])[R])+w[[2]] 2(Subscript[h, 2])[R]+2(((1-w[[1]]-w[[2]]))^2)(Subscript[J, 11])[R]+((w[[1]])^2)(Subscript[J, 12])[R]+2((w[[2]])^2)(Subscript[J, 22])[R]+2(1-w[[1]]-w[[2]])w[[2]] (Subscript[J, 12])[R]+NIntegrate[ϵxcLDA[nw[{x,y,z},R,w]]nw[{x,y,z},R,w],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]+(Subscript[Ε, NN])[R];*)

(*(Subscript[ΔΕw, LDA])[R_,w_]:=(1-w[[1]]-w[[2]])(Subscript[Ε, LDA])[R][[1]]+w[[1]](Subscript[Ε, LDA])[R][[2]]+w[[2]](Subscript[Ε, LDA])[R][[3]]-(Subscript[Εw, LDA])[R,w];*)

(*(Subscript[ϵw, LDA])[R_,w_]:={
(Subscript[h, 1])[R]+2(1-w)(Subscript[J, 11])[R]+2w (Subscript[J, 12])[R]+(1/2)NIntegrate[(dϵxcLDAdn[nw[{x,y,z},R,w]]nw[{x,y,z},R,w]+ϵxcLDA[nw[{x,y,z},R,w]])(Subscript[n, 0])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}],
(Subscript[h, 1])[R]+(Subscript[h, 2])[R]+(1-w[[1]]-w[[2]])(Subscript[J, 11])[R]+2w (Subscript[J, 12])[R]+(1/2)NIntegrate[(dϵxcLDAdn[nw[{x,y,z},R,w]]nw[{x,y,z},R,w]+ϵxcLDA[nw[{x,y,z},R,w]])(Subscript[n, 0])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}],
(Subscript[h, 2])[R]+2(1-w)(Subscript[J, 12])[R]+2w (Subscript[J, 22])[R]+(1/2)NIntegrate[(dϵxcLDAdn[nw[{x,y,z},R,w]]nw[{x,y,z},R,w]+ϵxcLDA[nw[{x,y,z},R,w]])(Subscript[n, 2])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]
};*)

eLDA energies

(*(Subscript[Ε, eLDA])[R_]:={
2(Subscript[h, 1])[R]+2(Subscript[J, 11])[R]+NIntegrate[ϵxceLDA[(Subscript[n, 0])[{x,y,z},R],0](Subscript[n, 0])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}],
(Subscript[h, 1])[R]+(Subscript[h, 2])[R]+(Subscript[J, 12])[R]+NIntegrate[ϵxceLDA[(Subscript[n, 1])[{x,y,z},R],0](Subscript[n, 1])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}],
2(Subscript[h, 2])[R]+2(Subscript[J, 22])[R]+NIntegrate[ϵxceLDA[(Subscript[n, 2])[{x,y,z},R],1](Subscript[n, 2])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]
}+(Subscript[Ε, NN])[R];*)

(*(Subscript[Εw, eLDA])[R_,w_]:=(1-w[[1]]-w[[2]])2(Subscript[h, 1])[R]+w[[1]]((Subscript[h, 1])[R]+(Subscript[h, 2])[R])+w[[2]] 2(Subscript[h, 2])[R]+2(((1-w[[1]]-w[[2]]))^2)(Subscript[J, 11])[R]+((w[[1]])^2)(Subscript[J, 12])[R]+2((w[[2]])^2)(Subscript[J, 22])[R]+2(1-w[[1]]-w[[2]])w[[2]] (Subscript[J, 12])[R]+NIntegrate[ϵxceLDA[nw[{x,y,z},R,w],w]nw[{x,y,z},R,w],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]*)

(*(Subscript[ΔΕw, eLDA])[R_,w_]:=(1-w[[1]]-w[[2]])(Subscript[Ε, eLDA])[R][[1]]+w[[1]](Subscript[Ε, eLDA])[R][[2]]+w[[2]](Subscript[Ε, eLDA])[R][[3]]-(Subscript[Εw, eLDA])[R,w];*)

(*(Subscript[ϵw, eLDA])[R_,w_]:={
(Subscript[h, 1])[R]+2(1-w)(Subscript[J, 11])[R]+2w (Subscript[J, 12])[R]+(1/2)NIntegrate[(dϵxceLDAdn[nw[{x,y,z},R,w],w]nw[{x,y,z},R,w]+ϵxceLDA[nw[{x,y,z},R,w],w])(Subscript[n, 1])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}],
(Subscript[h, 2])[R]+2(1-w)(Subscript[J, 12])[R]+2w (Subscript[J, 22])[R]+(1/2)NIntegrate[(dϵxceLDAdn[nw[{x,y,z},R,w],w]nw[{x,y,z},R,w]+ϵxceLDA[nw[{x,y,z},R,w],w])(Subscript[n, 2])[{x,y,z},R],{x,-∞,∞},{y,-∞,∞},{z,-∞,∞}]
};*)

CID energies

(*(Subscript[Ε, CISD])[R_]={(((Subscript[Ε, HF])[R][[1]]+(Subscript[Ε, HF])[R][[3]])/2)-(1/2)(Sqrt[((((Subscript[Ε, HF])[R][[1]]-(Subscript[Ε, HF])[R][[3]]))^2)+4(((Subscript[K, 12])[R])^2)]),(Subscript[h, 1])[R]+(Subscript[h, 2])[R]+(Subscript[J, 12])[R]+(Subscript[K, 12])[R],(((Subscript[Ε, HF])[R][[1]]+(Subscript[Ε, HF])[R][[3]])/2)+(1/2)(Sqrt[((((Subscript[Ε, HF])[R][[1]]-(Subscript[Ε, HF])[R][[3]]))^2)+4(((Subscript[K, 12])[R])^2)])};*)

Data

(*With[{Rmin=1.3,Rmax=1.5,ΔR=0.1},
ScanEHF=Table[Flatten[{R,(Subscript[Ε, HF])[R]}],{R,Rmin,Rmax,ΔR}];
ScanEwHF=Table[Flatten[{R,(Subscript[Εw, HF])[R,{(1/3),(1/3)}]}],{R,Rmin,Rmax,ΔR}];
ScanECID=Table[Flatten[{R,(Subscript[Ε, CID])[R]}],{R,Rmin,Rmax,ΔR}];
ScanELDA=Table[Flatten[{R,(Subscript[Ε, LDA])[R]}],{R,Rmin,Rmax,ΔR}];
ScanEwLDA=Table[Flatten[{R,(Subscript[Εw, LDA])[R,{(1/3),(1/3)}]}],{R,Rmin,Rmax,ΔR}];
ScanEeLDA=Table[Flatten[{R,(Subscript[Ε, eLDA])[R]}],{R,Rmin,Rmax,ΔR}];
ScanEweLDA=Table[Flatten[{R,(Subscript[Εw, eLDA])[R,{(1/3),(1/3)}]}],{R,Rmin,Rmax,ΔR}];
]*)

(*With[{R=1.4},(Subscript[Ε, HF])[R][[2;;3]]-(Subscript[Ε, HF])[R][[1]]]
With[{R=1.4},(Subscript[Ε, CISD])[R][[2;;3]]-(Subscript[Ε, CISD])[R][[1]]]
With[{R=1.4},(Subscript[Ε, LDA])[R][[2;;3]]-(Subscript[Ε, LDA])[R][[1]]]*)

Self-consistent ensemble KS calculation

Functionals

LDA exchange functionals

ϵxLDA[n_]=Module[{CxLDA=-(3/4)((((3/π)))^(1/3))},CxLDA (n^(1/3))];
ϵx0[n_]=Module[{Cx0=-(4/3) ((((1/π)))^(1/3))},Cx0 (n^(1/3))];
ϵx1[n_]=Module[{Cx1=-(16/10) ((((1/π)))^(1/3))},Cx1 (n^(1/3))];
ϵx2[n_]=Module[{Cx1=-(176/105) ((((1/π)))^(1/3))},Cx1 (n^(1/3))];
ϵxeLDA[n_,{w1_,w2_}]=ϵxLDA[n]+w1(ϵx1[n]-ϵx0[n])+w2(ϵx2[n]-ϵx0[n]);

dϵxLDAdn[n_]=Evaluate[D[ϵxLDA[n],n]];
dϵx0dn[n_]=Evaluate[D[ϵx0[n],n]];
dϵx1dn[n_]=Evaluate[D[ϵx1[n],n]];
dϵx2dn[n_]=Evaluate[D[ϵx2[n],n]];
dϵxeLDAdn[n_,{w1_,w2_}]=dϵxLDAdn[n]+w1(dϵx1dn[n]-dϵx0dn[n])+w2(dϵx2dn[n]-dϵx0dn[n]);

LDA correlation functionals

(* Intermediate quantities *)
X[x_,b_,c_]=(x^2)+b x+c;
x[d_]=(Sqrt[d]);
Q[b_,c_]=(Sqrt[4c-(b^2)]);
f[ζ_]=(((((1+ζ))^(4/3))+(((1-ζ))^(4/3))-2)/(2((2^(1/3))-1)));
(*interpolation function*)
vwn[A_,a_,b_,c_,d_]=A(Log[(((x[d])^2)/(X[x[d],b,c]))]+((2b)/(Q[b,c]))ArcTan[((Q[b,c])/(2x[d]+b))]-((b*a)/(X[a,b,c]))(Log[((((x[d]-a))^2)/(X[x[d],b,c]))]+((2(b+2a))/(Q[b,c]))ArcTan[((Q[b,c])/(2x[d]+b))]));

(*VWN5 correlation functional*)

ϵcVWN5[rs_,ζ_]=With[{f2=(4/(9((2^(1/3))-1)))},vwn[(0.0621814/2),-0.10498,3.72744,12.9352,rs]+vwn[-(1/(6(π^2))),-0.00475840,1.13107,13.0045,rs](((f[ζ])/f2))(((1-ζ))^4)+(vwn[(0.0310907/2),-0.32500,7.06042,18.0578,rs]-vwn[(0.0621814/2),-0.10498,3.72744,12.9352,rs])f[ζ](ζ^4)];

ϵcLDA[n_]=ϵcVWN5[((((3/(4π n))))^(1/3)),0];
ϵc0[n_]=Module[{a=-0.0238184`,b=0.00540994,c=0.0830766},(a/(1+b (n^(-1/6))+c (n^(-1/3))))];
ϵc1[n_]=Module[{a=-0.0282814`,b=0.002739250132910031`,c=0.06649138204575991`},(a/(1+b (n^(-1/6))+c (n^(-1/3))))];
ϵc2[n_]=Module[{a=-0.0144633`,b=-0.0506019,c=0.0331417},(a/(1+b (n^(-1/6))+c (n^(-1/3))))];
ϵceLDA[n_,{w1_,w2_}]=ϵcLDA[n]+w1(ϵc1[n]-ϵc0[n])+w2(ϵc2[n]-ϵc0[n]);

dϵcLDAdn[n_]=Evaluate[D[ϵcLDA[n],n]];
dϵc0dn[n_]=Evaluate[D[ϵc0[n],n]];
dϵc1dn[n_]=Evaluate[D[ϵc1[n],n]];
dϵc2dn[n_]=Evaluate[D[ϵc2[n],n]];
dϵceLDAdn[n_,{w1_,w2_}]=dϵcLDAdn[n]+w1(dϵc1dn[n]-dϵc0dn[n])+w2(dϵc2dn[n]-dϵc0dn[n]);

LDA exchange-correlation functionals

ϵxcLDA[n_]=ϵxLDA[n]+ϵcLDA[n];
ϵxc0[n_]=ϵx0[n]+ϵc0[n];
ϵxc1[n_]=ϵx1[n]+ϵc1[n];
ϵxc2[n_]=ϵx2[n]+ϵc2[n];
ϵxceLDA[n_,{w1_,w2_}]=ϵxeLDA[n,{w1,w2}]+ϵceLDA[n,{w1,w2}];

dϵxcLDAdn[n_]=dϵxLDAdn[n]+dϵcLDAdn[n];
dϵxc0dn[n_]=dϵx0dn[n]+dϵc0dn[n];
dϵxc1dn[n_]=dϵx1dn[n]+dϵc1dn[n];
dϵxc2dn[n_]=dϵx2dn[n]+dϵc2dn[n];
dϵxceLDAdn[n_,{w1_,w2_}]=dϵxeLDAdn[n,{w1,w2}]+dϵceLDAdn[n,{w1,w2}];

dϵxeLDAdw1[n_]=ϵx1[n]-ϵx0[n];
dϵceLDAdw1[n_]=ϵc1[n]-ϵc0[n];
dϵxceLDAdw1[n_]=ϵxc1[n]-ϵxc0[n];
dϵxeLDAdw2[n_]=ϵx2[n]-ϵx0[n];
dϵceLDAdw2[n_]=ϵc2[n]-ϵc0[n];
dϵxceLDAdw2[n_]=ϵxc2[n]-ϵxc0[n];

D[((a n)/(1+b (n^(-1/6))+c (n^(-1/3)))),n]==(a/(1+b (n^(-1/6))+c (n^(-1/3))))+n(1/(6 n))(a/(1+b (n^(-1/6))+c (n^(-1/3))))((b (n^(-1/6))+2c (n^(-1/3)))/(1+b (n^(-1/6))+c (n^(-1/3))))//FullSimplify

D[((a 2n)/(1+b (((2n))^(-1/6))+c (((2n))^(-1/3)))),n]==((a/(1+b (((2n))^(-1/6))+c (((2n))^(-1/3))))+2n(1/(6 2n))(a/(1+b (((2n))^(-1/6))+c (((2n))^(-1/3))))((b (((2n))^(-1/6))+2c (((2n))^(-1/3)))/(1+b (((2n))^(-1/6))+c (((2n))^(-1/3)))))2//FullSimplify

One- and two-electron integrals

Overlap

𝒮[{α_,A_,LA_},{β_,B_,LB_}]:=𝒮[{α,A,LA},{β,B,LB}]=Module[{EAB,Overlap},
EAB=Exp[-((α β)/(α+β))(A-B).(A-B)];

(* Lazy evaluation *)
sRR[i_,a_,b_]:=sRR[i,a,b];

Do[
(* Initial Conditions *)
sRR[i_,0,0]:=1;
sRR[i_,1,0]:=-(A[[i]]-((α A[[i]]+β B[[i]])/(α+β)));
(* Recurrence Index *)
sRR[i_,a_,0]:=-(A[[i]]- ((α A[[i]]+β B[[i]])/(α+β)) )sRR[i,a-1,0]+((a-1)/(2 (α+β)))sRR[i,a-2,0];
(* Transfer Recurrence Relation *)
sRR[i_,a_,b_]:= sRR[i,a+1,b-1]+(A[[i]]-B[[i]])sRR[i,a,b-1];

,{i,3}];
Return[EAB((((π/(α+β))))^(3/2))(∏_(i=1)^3)sRR[i,LA[[i]],LB[[i]]]];
];

𝒮c[{α_,A_,LA_,DA_},{β_,B_,LB_,DB_}]:=Module[{KA,KB,Int},
(* Degrees of contraction *)
KA=Length[α];	KB=Length[β];
(* Contracted integral *)
Int=(∑_(kA=1)^KA)(∑_(kB=1)^KB)NormCoeff[α[[kA]],LA]NormCoeff[β[[kB]],LB]DA[[kA]]DB[[kB]]𝒮[{α[[kA]],A,LA},{β[[kB]],B,LB}];
Return[Int];
];

ERI

G[m_,ζ_,η_,NPQ_]:=G[m,ζ,η,NPQ]=
Which[
NPQ==0,
((((-((ζ η)/(ζ+η))))^m)/((1+2 m) ((ζ η)/(ζ+η)))),
(NPQ^2) ((ζ η)/(ζ+η))<0.1,(1/(2 ((ζ η)/(ζ+η))))(((-((ζ η)/(ζ+η))))^m) ((1/((1/2)+m))-(((ζ) (η) (NPQ^2))/(((3/2)+m) (ζ+η)))+(((((ζ))^2) (((η))^2) (NPQ^4))/((5+2 m) (((ζ+η))^2)))-(((((ζ))^3) (((η))^3) (NPQ^6))/(6 ((7/2)+m) (((ζ+η))^3)))+(((((ζ))^4) (((η))^4) (NPQ^8))/(24 ((9/2)+m) (((ζ+η))^4)))-(((((ζ))^5) (((η))^5) (NPQ^10))/(120 ((11/2)+m) (((ζ+η))^5)))+(((((ζ))^6) (((η))^6) (NPQ^12))/(720 ((13/2)+m) (((ζ+η))^6)))),
True,
(1/(2 ((ζ η)/(ζ+η))))(((-((ζ η)/(ζ+η))))^m) ((((NPQ^2) ((ζ η)/(ζ+η))))^(-(1/2)-m)) (Γ[(1/2)+m]-Γ[(1/2)+m,(NPQ^2) ((ζ η)/(ζ+η))])
];

𝒰[{α_,A_,LA_},{β_,B_,LB_},{γ_,C_,LC_},{δ_,D_,LD_}]:=𝒰[{α,A,LA},{β,B,LB},{γ,C,LC},{δ,D,LD}]=Module[{AB,CD,P,Q,PQ,NAB,NCD,NPQ,ζ,η,Gp,i0,i1,i2,i3,int},
AB=A-B;	CD=C-D;	NAB=Norm[AB];	NCD=Norm[CD];
ζ=α+β;	η=γ+δ;	P=((α A+β B)/ζ);	Q=((γ C+δ D)/η);	PQ=P-Q;	NPQ=Norm[PQ];
Gp=((2(π^(5/2)))/(((ζ+η))^(3/2)))Exp[-((α β)/ζ)(NAB^2)-((γ δ)/η)(NCD^2)];
i0={0,0,0};	i1={1,0,0};	i2={0,1,0};	i3={0,0,1};
(* Initial Conditions *)
Do[G[m,ζ,η,NPQ],{m,0,(∑_(μ=1)^3)(LA[[μ]]+LB[[μ]]+LC[[μ]]+LD[[μ]])}];
(* (([0])^m) *)
iRR[m_,a_,b_,c_,d_]:=iRR[m,a,b,c,d];
iRR[m_,a_,b_,c_,d_]:=0 /; a[[1]]<0||b[[1]] <0||c[[1]] <0||d[[1]] <0||a[[2]]<0|| b[[2]] <0|| c[[2]]<0|| d[[2]]<0||a[[3]]<0||  b[[3]]<0|| c[[3]]<0|| d[[3]]<0;
iRR[m_,i0,i0,i0,i0]:= G[m,ζ,η,NPQ];
(* VRR [00|0d] *)
iRR[m_,i0,i0,i0,d_]:=Which[
d[[1]]>0,
((γ CD[[1]])/η)iRR[m,i0,i0,i0,d-i1]-((PQ[[1]])/η)iRR[m+1,i0,i0,i0,d-i1]+((d[[1]]-1)/(2η))(iRR[m,i0,i0,i0,d-2i1]+(1/η)iRR[m+1,i0,i0,i0,d-2 i1]),
d[[2]]>0,
((γ CD[[2]])/η)iRR[m,i0,i0,i0,d-i2]-((PQ[[2]])/η)iRR[m+1,i0,i0,i0,d-i2]+((d[[2]]-1)/(2η))(iRR[m,i0,i0,i0,d-2i2]+(1/η)iRR[m+1,i0,i0,i0,d-2 i2]),
d[[3]]>0,
((γ CD[[3]])/η)iRR[m,i0,i0,i0,d-i3]-((PQ[[3]])/η)iRR[m+1,i0,i0,i0,d-i3]+((d[[3]]-1)/(2η))(iRR[m,i0,i0,i0,d-2i3]+(1/η)iRR[m+1,i0,i0,i0,d-2 i3])];
(* VRR [0b|0d] *)
iRR[m_,i0,b_,i0,d_]:=Which[
b[[1]]>0,
((α AB[[1]])/ζ)iRR[m,i0,b-i1,i0,d]+((PQ[[1]])/ζ)iRR[m+1,i0,b-i1,i0,d]+((b[[1]]-1)/(2ζ))(iRR[m,i0,b-2i1,i0,d]+(1/ζ)iRR[m+1,i0,b-2 i1,i0,d])-((d[[1]])/(2ζ η))iRR[m+1,i0,b-i1,i0,d-i1],
b[[2]]>0,
((α AB[[2]])/ζ)iRR[m,i0,b-i2,i0,d]+((PQ[[2]])/ζ)iRR[m+1,i0,b-i2,i0,d]+((b[[2]]-1)/(2ζ))(iRR[m,i0,b-2i2,i0,d]+(1/ζ)iRR[m+1,i0,b-2 i2,i0,d])-((d[[2]])/(2ζ η))iRR[m+1,i0,b-i2,i0,d-i2],
b[[3]]>0,
((α AB[[3]])/ζ)iRR[m,i0,b-i3,i0,d]+((PQ[[3]])/ζ)iRR[m+1,i0,b-i3,i0,d]+((b[[3]]-1)/(2ζ))(iRR[m,i0,b-2i3,i0,d]+(1/ζ)iRR[m+1,i0,b-2 i3,i0,d])-((d[[3]])/(2ζ η))iRR[m+1,i0,b-i3,i0,d-i3]];
(* TRR [0b|cd] *)
iRR[m_,i0,b_,c_,d_]:=Which[
c[[1]]>0,
iRR[m,i0,b,c-i1,d+i1]-CD[[1]]iRR[m,i0,b,c-i1,d],
c[[2]]>0,
iRR[m,i0,b,c-i2,d+i2]-CD[[2]]iRR[m,i0,b,c-i2,d],
c[[3]]>0,
iRR[m,i0,b,c-i3,d+i3]-CD[[3]]iRR[m,i0,b,c-i3,d]];
(* TRR [ab|cd] *)
iRR[m_,a_,b_,c_,d_]:=Which[
a[[1]]>0,
iRR[m,a-i1,b+i1,c,d]-AB[[1]]iRR[m,a-i1,b,c,d],
a[[2]]>0,
iRR[m,a-i2,b+i2,c,d]-AB[[2]]iRR[m,a-i2,b,c,d],
a[[3]]>0,
iRR[m,a-i3,b+i3,c,d]-AB[[3]]iRR[m,a-i3,b,c,d]];
Return[N[Gp iRR[0,LA,LB,LC,LD]]];
];

𝒰c[{α_,A_,LA_,DA_},{β_,B_,LB_,DB_},{γ_,C_,LC_,DC_},{δ_,D_,LD_,DD_}]:=Module[{KA,KB,KC,KD,Int},
(* Degrees of contraction *)
KA=Length[α];	KB=Length[β];	KC=Length[γ];	KD=Length[δ];
(* Contracted integral *)
Int=(∑_(kA=1)^KA)(∑_(kB=1)^KB)(∑_(kC=1)^KC)(∑_(kD=1)^KD)NormCoeff[α[[kA]],LA]NormCoeff[β[[kB]],LB]NormCoeff[γ[[kC]],LC]NormCoeff[δ[[kD]],LD]DA[[kA]]DB[[kB]]DC[[kC]]DD[[kD]]𝒰[{α[[kA]],A,LA},{β[[kB]],B,LB},{γ[[kC]],C,LC},{δ[[kD]],D,LD}];
Return[Int];
];

F12

pf12G[{α_,A_,LA_},{β_,B_,LB_},{γ_,C_,LC_},{δ_,D_,LD_},λ_]:=pf12G[{α,A,LA},{β,B,LB},{γ,C,LC},{δ,D,LD},λ]=Module[{ζ,η,P,Q,GAB,GCD,AB,CD,PQ,fζ,fλ,fη,O,x,y,z,gζ,gλ,gη,RA,RC},

(* Pre-computed quantities *)
ζ=α+β;	P=((α A+β B)/ζ);	AB=A-B;	GAB=((((π/ζ)))^(3/2))Exp[-(((Norm[AB])^2)/((α^(-1))+(β^(-1))))];
η=γ+δ;	Q=((γ C+δ D)/η);	CD=C-D;	GCD=((((π/η)))^(3/2))Exp[-(((Norm[CD])^2)/((γ^(-1))+(δ^(-1))))];
PQ=P-Q;	fζ=((ζ^(-1))/((ζ^(-1))+(λ^(-1))+(η^(-1))));	fλ=((λ^(-1))/((ζ^(-1))+(λ^(-1))+(η^(-1))));	fη=((η^(-1))/((ζ^(-1))+(λ^(-1))+(η^(-1))));
gζ=((1-fζ)/(2ζ));	gλ=((fζ)/(2η));	gη=((1-fη)/(2η));	RA=(P-A)+fζ (Q-P);	RC=(Q-C)+fη (P-Q);

(* Lazy evaluation *)
Gabcd[i_,a_,b_,c_,d_]:=Gabcd[i,a,b,c,d];

(* Termination condition *)
Gabcd[i_,a_,b_,c_,d_]:=0 /; a<0||b <0||c <0||d <0;

(* Loop over cartesian directions *)
Do[

(* [00|00] *)
Gabcd[i,0,0,0,0]:=((fλ)^(1/2))Exp[-((((P[[i]]-Q[[i]]))^2)/((ζ^(-1))+(λ^(-1))+(η^(-1))))];

(* [(a+1)0|00] *)
Gabcd[i_,a_,0,0,0]:=RA[[i]]Gabcd[i,a-1,0,0,0]+(a-1)gζ Gabcd[i,a-2,0,0,0];

(* [a0|(c+1)0] *)
Gabcd[i_,a_,0,c_,0]:=RC[[i]]Gabcd[i,a,0,c-1,0]+a gλ Gabcd[i,a-1,0,c-1,0]+(c-1)gη Gabcd[i,a,0,c-2,0];

(* [a(b+1)|c0] *)
Gabcd[i_,a_,b_,c_,0]:=Gabcd[i,a+1,b-1,c,0]+AB[[i]]Gabcd[i,a,b-1,c,0];

(* [ab|c(d+1)] *)
Gabcd[i,a_,b_,c_,d_]:=Gabcd[i,a,b,c+1,d-1]+CD[[i]]Gabcd[i,a,b,c,d-1];

,{i,3}];

(* Final result *)
Return[N[GAB GCD (∏_(i=1)^3)Gabcd[i,LA[[i]],LB[[i]],LC[[i]],LD[[i]]],15]];
];

pf12S[{α_,A_,LA_},{β_,B_,LB_},{γ_,C_,LC_},{δ_,D_,LD_},λ_]:=pf12S[{α,A,LA},{β,B,LB},{γ,C,LC},{δ,D,LD},λ]=Module[{αSG,cSG,Int},

(* Exponents and coefficients for the Slater function *)
αSG={0.2209,1.004,3.622,12.16,45.87,254.4};
cSG={0.3144,0.3037,0.1681,0.09811,0.06024,0.03726};

(* Contraction *)
Int=(∑_(k=1)^6)cSG[[k]]pf12G[{α,A,LA},{β,B,LB},{γ,C,LC},{δ,D,LD},(λ^2)αSG[[k]]];

(* Final result *)
Return[N[Int,15]];
];

cf12S[{α_,A_,LA_,DA_},{β_,B_,LB_,DB_},{γ_,C_,LC_,DC_},{δ_,D_,LD_,DD_},λ_]:=Module[{KA,KB,KC,KD,Int},
(* Degrees of contraction *)
KA=Length[α];	KB=Length[β];	KC=Length[γ];	KD=Length[δ];
(* Contracted integral *)
Int=(∑_(kA=1)^KA)(∑_(kB=1)^KB)(∑_(kC=1)^KC)(∑_(kD=1)^KD)NormCoeff[α[[kA]],LA]NormCoeff[β[[kB]],LB]NormCoeff[γ[[kC]],LC]NormCoeff[δ[[kD]],LD]DA[[kA]]DB[[kB]]DC[[kC]]DD[[kD]]pf12S[{α[[kA]],A,LA},{β[[kB]],B,LB},{γ[[kC]],C,LC},{δ[[kD]],D,LD},λ];
Return[Int];
];

With[{α=1/2.,β=2.,γ=3.,δ=1.,A={0.,0.,0.},B={1.,0.,0.},C={0.,0.,1/2.},D={1/2.,-3/2.,0.},LA={1,1,1},LB={0,1,1},LC={2,1,1},LD={2,1,2},λ=1.4,r1={x1,y1,z1},r2={x2,y2,z2}},
{pf12S[{α,A,LA},{β,B,LB},{γ,C,LC},{δ,D,LD},λ],
(∏_(k=1)^3)NIntegrate[(((x1-A[[k]]))^(LA[[k]]))Exp[-α ((Norm[x1-A[[k]]])^2)](((x1-B[[k]]))^(LB[[k]]))Exp[-β ((Norm[x1-B[[k]]])^2)]Exp[-λ ((Norm[x1-x2])^2)](((x2-C[[k]]))^(LC[[k]]))Exp[-γ ((Norm[x2-C[[k]]])^2)](((x2-D[[k]]))^(LD[[k]]))Exp[-δ ((Norm[x2-D[[k]]])^2)],{x1,-∞,∞},{x2,-∞,∞}]}]

Potential

GEXPV[α_,β_,A_,B_]:=((((π/(α+β))))^(3/2))(2/(Sqrt[π]))Exp[-((α β)/(α+β))((Norm[A-B])^2)];

GmV[m_,ζ_,A_,B_]:=If[((Norm[A-B])^2)==0,((((ζ))^((1/2)+m))/(1+2 m)),If[((Norm[A-B])^2)(ζ)<0.1,(1/2)(((ζ))^((1/2)+m))((1/((1/2)+m))-(((ζ) ((Norm[A-B])^2))/((3/2)+m))+(((((ζ))^2) ((Norm[A-B])^4))/(5+2 m))-(((((ζ))^3) ((Norm[A-B])^6))/(6 ((7/2)+m)))+(((((ζ))^4) ((Norm[A-B])^8))/(24 ((9/2)+m)))-(((((ζ))^5) ((Norm[A-B])^10))/(120 ((11/2)+m)))+(((((ζ))^6) ((Norm[A-B])^12))/(720 ((13/2)+m)))-(((((ζ))^7) ((Norm[A-B])^14))/(5040 ((15/2)+m)))),(1/2)(((ζ))^((1/2)+m))(((((Norm[A-B])^2)(ζ)))^(-(1/2)-m))(Gamma[(1/2)+m]-Gamma[(1/2)+m,((Norm[A-B])^2)(ζ)])]]

VR[m_,a_,b_,α_,β_,A_,B_,R_]:=
Which[ Norm[b]>0,Which[(Subscript[b, [[1]]])>0,VR[m,a+(Subscript[Id, [[1;;3,1]]]),b-(Subscript[Id, [[1;;3,1]]]),α,β,A,B,R]+((Subscript[A, [[1]]])-(Subscript[B, [[1]]]))VR[m,a,b-(Subscript[Id, [[1;;3,1]]]),α,β,A,B,R],
(Subscript[b, [[2]]])>0,VR[m,a+(Subscript[Id, [[1;;3,2]]]),b-(Subscript[Id, [[1;;3,2]]]),α,β,A,B,R]+((Subscript[A, [[2]]])-(Subscript[B, [[2]]]))VR[m,a,b-(Subscript[Id, [[1;;3,2]]]),α,β,A,B,R],
(Subscript[b, [[3]]])>0,VR[m,a+(Subscript[Id, [[1;;3,3]]]),b-(Subscript[Id, [[1;;3,3]]]),α,β,A,B,R]+((Subscript[A, [[3]]])-(Subscript[B, [[3]]]))VR[m,a,b-(Subscript[Id, [[1;;3,3]]]),α,β,A,B,R]],
 Norm[a]>0,
Which[(Subscript[a, [[1]]])>0,(((α (Subscript[A, [[1]]])+β (Subscript[B, [[1]]]))/(α+β))-(Subscript[A, [[1]]]))VR[m,a-(Subscript[Id, [[1;;3,1]]]),b,α,β,A,B,R]+(((Subscript[a, [[1]]])-1)/(2 (α+β)))VR[m,a-2(Subscript[Id, [[1;;3,1]]]),b,α,β,A,B,R]-(1/(α+β))(((α (Subscript[A, [[1]]])+β (Subscript[B, [[1]]]))/(α+β))-(Subscript[R, [[1]]]))VR[m+1,a-(Subscript[Id, [[1;;3,1]]]),b,α,β,A,B,R]-(((Subscript[a, [[1]]])-1)/(2 (((α+β))^2)))VR[m+1,a-2(Subscript[Id, [[1;;3,1]]]),b,α,β,A,B,R],
(Subscript[a, [[2]]])>0,(((α (Subscript[A, [[2]]])+β (Subscript[B, [[2]]]))/(α+β))-(Subscript[A, [[2]]]))VR[m,a-(Subscript[Id, [[1;;3,2]]]),b,α,β,A,B,R]+(((Subscript[a, [[2]]])-1)/(2 (α+β)))VR[m,a-2(Subscript[Id, [[1;;3,2]]]),b,α,β,A,B,R]-(1/(α+β))(((α (Subscript[A, [[2]]])+β (Subscript[B, [[2]]]))/(α+β))-(Subscript[R, [[2]]]))VR[m+1,a-(Subscript[Id, [[1;;3,2]]]),b,α,β,A,B,R]-(((Subscript[a, [[2]]])-1)/(2 (((α+β))^2)))VR[m+1,a-2(Subscript[Id, [[1;;3,2]]]),b,α,β,A,B,R],
(Subscript[a, [[3]]])>0,(((α (Subscript[A, [[3]]])+β (Subscript[B, [[3]]]))/(α+β))-(Subscript[A, [[3]]]))VR[m,a-(Subscript[Id, [[1;;3,3]]]),b,α,β,A,B,R]+(((Subscript[a, [[3]]])-1)/(2 (α+β)))VR[m,a-2(Subscript[Id, [[1;;3,3]]]),b,α,β,A,B,R]-(1/(α+β))(((α (Subscript[A, [[3]]])+β (Subscript[B, [[3]]]))/(α+β))-(Subscript[R, [[3]]]))VR[m+1,a-(Subscript[Id, [[1;;3,3]]]),b,α,β,A,B,R]-(((Subscript[a, [[3]]])-1)/(2 (((α+β))^2)))VR[m+1,a-2(Subscript[Id, [[1;;3,3]]]),b,α,β,A,B,R]],
(a=={0,0,0})&&(b=={0,0,0}),
GEXPV[α,β,A,B]GmV[m,α+β,((α A+β B)/(α+β)),R]];
VR[m_,a_,b_,α_,β_,A_,B_,R_]:=0/;Norm[a]<0||Norm[b]<0

𝒱[{α_,A_,a_},{β_,B_,b_},{ZC_,C_}]:=𝒱[{α,A,a},{β,B,b},{ZC,C}]=-ZC VR[0,a,b,α,β,A,B,C];

𝒱c[{α_,A_,LA_,DA_},{β_,B_,LB_,DB_},{ZC_,C_}]:=Module[{KA,KB,Int},
(* Degrees of contraction *)
KA=Length[α];	KB=Length[β];
(* Contracted integral *)
Int=(∑_(kA=1)^KA)(∑_(kB=1)^KB)NormCoeff[α[[kA]],LA]NormCoeff[β[[kB]],LB]DA[[kA]]DB[[kB]]𝒱[{α[[kA]],A,LA},{β[[kB]],B,LB},{ZC,C}];
Return[Int];
];

Kinetic

DD[dx_,dy_,dz_,a_,b_,α_,β_,A_,B_]:=Which[dx>0,2β DD[dx-1,dy,dz,a,b+ Id[[1 ;; 3,1]],α,β,A,B]-(Subscript[b, [[1]]])DD[dx-1,dy,dz,a,b- Id[[1 ;; 3,1]],α,β,A,B],dy>0,2β DD[dx,dy-1,dz,a,b+ Id[[1 ;; 3,2]],α,β,A,B]-(Subscript[b, [[2]]])DD[dx,dy-1,dz,a,b- Id[[1 ;; 3,2]],α,β,A,B],dz>0,2β DD[dx,dy,dz-1,a,b+ Id[[1 ;; 3,3]],α,β,A,B]-(Subscript[b, [[3]]])DD[dx,dy,dz-1,a,b- Id[[1 ;; 3,3]],α,β,A,B]];
DD[dx_,dy_,dz_,a_,b_,α_,β_,A_,B_]:=0/; dx<0||dy<0||dz<0||(Subscript[a, [[1]]])<0||(Subscript[a, [[2]]])<0||(Subscript[a, [[3]]])<0||(Subscript[b, [[1]]])<0||(Subscript[b, [[2]]])<0||(Subscript[b, [[3]]])<0;
DD[dx_,dy_,dz_,a_,b_,α_,β_,A_,B_]:=-(1/2)𝒮[{α,A,a},{β,B,b}]/; dx==0&&dy==0&&dz==0;
𝒯[{α_,A_,a_},{β_,B_,b_}]:=𝒯[{α,A,a},{β,B,b}]=DD[2,0,0,a,b,α,β,A,B]+DD[0,2,0,a,b,α,β,A,B]+DD[0,0,2,a,b,α,β,A,B]

𝒯c[{α_,A_,LA_,DA_},{β_,B_,LB_,DB_}]:=Module[{KA,KB,Int},
(* Degrees of contraction *)
KA=Length[α];	KB=Length[β];
(* Contracted integral *)
Int=(∑_(kA=1)^KA)(∑_(kB=1)^KB)NormCoeff[α[[kA]],LA]NormCoeff[β[[kB]],LB]DA[[kA]]DB[[kB]]𝒯[{α[[kA]],A,LA},{β[[kB]],B,LB}];
Return[Int];
];

HF-related modules

Contracted integrals module

Int[NAt_,NB_,Atoms_,BF_]:=Module[{S,s,T,V,U,X,WP,u,v,w,ϵ,Cp,Ortho,nLD,cutoff},
WP=15;	cutoff=(10^(-7));
Ortho="SVD";
S=N[Table[𝒮c[BF[[a]],BF[[b]]],{a,NB},{b,NB}],WP];
T=N[Table[𝒯c[BF[[a]],BF[[b]]],{a,NB},{b,NB}],WP];
V=N[Table[(∑_(c=1)^NAt)𝒱c[BF[[a]],BF[[b]],Atoms[[c]]],{a,NB},{b,NB}],WP];
U=N[Table[𝒰c[BF[[a]],BF[[b]],BF[[c]],BF[[d]]],{a,NB},{b,NB},{c,NB},{d,NB}],WP];
Which[
Ortho=="Lowdin",X=MatrixPower[S,-1/2],
Ortho=="Canonical",X=Eigenvectors[S]\[Transpose].DiagonalMatrix[((Eigenvalues[S])^(-1/2))],
Ortho=="SVD",(
	s=Sort[Eigenvalues[S]];
	Print["Overlap eigenvalues in OBS = ",s];
	nLD=0;
	Do[If[s[[i]]<cutoff ,nLD++],{i,NB}];
	Print["Number of linear dependencies in OBS = ",nLD];
	{u,w,v}=SingularValueDecomposition[S,NB-nLD];	
	X=u.DiagonalMatrix[((Diagonal[w])^(-1/2))].v\[Transpose];
)
];
Return[{S,T,V,U,X}];
]

Restricted  ensemble Hartree-Fock module (RHF)

RHF[n_,w_,ENuc_,N_,S_,T_,V_,U_,X_]:=Module[{maxSCF,thresh,Fw,H,UJ,UK,C,Cp,ϵ,Γ0,Γ1,Γw,nSCF,Conv,EwHF,EwT,EwV,EwJ,EwK,Gap,EwGICHF,EHF,ΩHF,TabRes},

(* Parameters *)
maxSCF=128;	thresh=(10^(-7));

(* Check weight conditions *)
If[!0≤w≤1/2,Print["Weights are inconsistent!!"]];

(* Generate initial guess *)
H=T+V;
UJ=Flatten[U,{{1},{2},{3,4}}];
UK=Flatten[U,{{1},{4},{3,2}}];
{ϵ,Cp}=Chop[SortEigensystem[Eigensystem[X\[Transpose].H.X]]];	Cp=Cp\[Transpose];
C=X.Cp;

(* Main SCF loop *)
Conv=1;	nSCF=0;
While[Conv>thresh∧nSCF<maxSCF,
nSCF++;

(* Density matrices *)
Γ0=2C\[Transpose][[1;;(n/2)]]\[Transpose].C\[Transpose][[1;;(n/2)]];	
Γ1=If[(n/2)>1,2C\[Transpose][[1;;(n/2)-1]]\[Transpose].C\[Transpose][[1;;(n/2)-1]],0]+2C\[Transpose][[(n/2)+1;;(n/2)+1]]\[Transpose].C\[Transpose][[(n/2)+1;;(n/2)+1]];
Γw=(1-w)Γ0+w Γ1;

(* Weight-dependent Fock matrix *)
Fw=H+(UJ-(1/2)UK).Flatten[Γw];

(* Weight-dependent MOs *)
{ϵ,Cp}=Chop[SortEigensystem[Eigensystem[X\[Transpose].Fw.X]]];	Cp=Cp\[Transpose];
C=X.Cp;

(* Convergence criterion *)
Conv=Max[Abs[Fw.Γw.S-S.Γw.Fw]];

];

(* Print HF information *)
Print["MO energies: ",Chop[{ϵ}]//TableForm];
Print["MO coefficients: ",Chop[C]//TableForm];
Gap=If[N>(n/2),ϵ[[(n/2)+1]]-ϵ[[(n/2)]],0];
EwT= Tr[Γw.T];	EwV= Tr[Γw.V];	EwJ= ((Tr[Γw.UJ.Flatten[Γw]])/2);	EwK= -((Tr[Γw.UK.Flatten[Γw]])/4);
EwHF=EwT+EwV+EwJ+EwK;

(* GIC-HF energies *)
EwGICHF=EwT+EwV+(1/2)(1-w)Tr[Γ0.(UJ-(1/2)UK).Flatten[Γ0]]+(1/2)w Tr[Γ1.(UJ-(1/2)UK).Flatten[Γ1]];
TabRes=TableForm[{{n,N,w,NumberForm[ENuc+EwHF,{9,6}],NumberForm[ENuc,{9,6}],NumberForm[EwHF,{9,6}],NumberForm[EwGICHF,{9,6}],NumberForm[EwT,{9,6}],NumberForm[EwV,{9,6}],NumberForm[EwJ,{9,6}],NumberForm[EwK,{9,6}],nSCF,NumberForm[Conv,{6,3}],NumberForm[Gap,{6,3}]}},
TableHeadings->{None,{"n","(Subscript[N, Basis])","w","(Subscript[E, 0])","(Subscript[E, Nuc])","(Subscript[E, HF]^w)","(Subscript[E, GIC - HF]^w)","(Subscript[T, HF]^w)","(Subscript[V, HF]^w)","(Subscript[J, HF]^w)","(Subscript[K, HF]^w)","(Subscript[N, SCF])","Conv","HL Gap"}},TableAlignments->Center];
Print[TabRes];

(* Compute individual energies for ground and excited states *) 
EHF={Tr[Γ0.H],Tr[Γ1.H]}
+(1/2){Tr[Γ0.(UJ-(1/2)UK).Flatten[Γ0]],Tr[Γ1.(UJ-(1/2)UK).Flatten[Γ1]]};

(* Excitation energies *)
ΩHF=EHF[[2]]-EHF[[1]];

(* Print results *)
TabRes=TableForm[{{
NumberForm[ENuc+EHF[[1]],{9,6}],NumberForm[ENuc+EHF[[2]],{9,6}],NumberForm[ΩHF,{9,6}],NumberForm[ΩHF HaToeV,{9,6}]
}},
TableHeadings->{None,{"(Subscript[E, HF]^((0)))","(Subscript[E, HF]^((1)))","(Subscript[Ω, HF]^((1)))","(Subscript[Ω, HF]^((1)))(eV)"}},TableAlignments->Center];
Print[TabRes];

(* Done *)
Return[{Γ0,Γ1}];

];

GOK-DFT module

GOKDFT[w_,ENuc_,T_,V_,U_,Γ0_,Γ1_,χμ_,Quad_]:=Module[{H,UJ,UK,Γw,n0,n1,nw,ELDA,EeLDA,ΩLDA,ΩeLDA,ΔDD,ΔΔDD,TabRes},

(* Generate useful quantities *)
H=T+V;
UJ=Flatten[U,{{1},{2},{3,4}}];
UK=Flatten[U,{{1},{4},{3,2}}];
Γw=(1-w)Γ0+w Γ1;

(* One-electron density *)
n0=χμ.Γ0.χμ;	n1=χμ.Γ1.χμ;	nw=(1-w)n0+w n1;
(* Check electron number *)
Print["Integrating one-electron density..."];
Print[NInt[nw,Quad]];

(* Compute individual energies for ground and excited states *)
ELDA={Tr[Γ0.H],Tr[Γ1.H]}
+(1/2){Tr[Γ0.UJ.Flatten[Γ0]],Tr[Γ1.UJ.Flatten[Γ1]]}
+Chop[NInt[(ϵxcLDA[nw]+dϵxcLDAdn[nw]nw){n0,n1}-dϵxcLDAdn[nw](nw^2),Quad]];

(* Excitation energies *)
ΩLDA=ELDA[[2]]-ELDA[[1]];

(* Print results *)
TabRes=TableForm[{{
NumberForm[ENuc+ELDA[[1]],{9,6}],NumberForm[ENuc+ELDA[[2]],{9,6}],NumberForm[ΩLDA,{9,6}],NumberForm[ΩLDA HaToeV,{9,6}]
}},
TableHeadings->{None,{"(Subscript[E, LDA]^((0)))","(Subscript[E, LDA]^((1)))","(Subscript[Ω, LDA]^((1)))","(Subscript[Ω, LDA]^((1))) (eV)"}},TableAlignments->Center];
Print[TabRes];

(* Compute individual energies for ground and excited states *)
EeLDA={Tr[Γ0.H],Tr[Γ1.H]}
+(1/2){Tr[Γ0.UJ.Flatten[Γ0]],Tr[Γ1.UJ.Flatten[Γ1]]}
+Chop[NInt[(ϵxceLDA[nw,{0,w}]+dϵxceLDAdn[nw,{0,w}]nw){n0,n1}-dϵxceLDAdn[nw,{0,w}](nw^2),Quad]];
Print[{Tr[Γ0.T],Tr[Γ1.T]}];
Print[{Tr[Γ0.V],Tr[Γ1.V]}];
Print[(1/2){Tr[Γ0.UJ.Flatten[Γ0]],Tr[Γ1.UJ.Flatten[Γ1]]}];
Print[Chop[NInt[(ϵxeLDA[nw,{0,w}]+dϵxeLDAdn[nw,{0,w}]nw){n0,n1}-dϵxeLDAdn[nw,{0,w}](nw^2),Quad]]];
Print[Chop[NInt[(ϵceLDA[nw,{0,w}]+dϵceLDAdn[nw,{0,w}]nw){n0,n1}-dϵceLDAdn[nw,{0,w}](nw^2),Quad]]];
(* Derivative discontinuities *)
(*ΔDD=Chop[NInt[ dϵxeLDAdw2[nw]nw ,Quad]];
ΔDD={-w ΔDD,(1-w)ΔDD};
ΔΔDD=ΔDD[[2]]-ΔDD[[1]];
Print["x",ΔDD];

ΔDD=Chop[NInt[ dϵceLDAdw2[nw]nw ,Quad]];
ΔDD={-w ΔDD,(1-w)ΔDD};
ΔΔDD=ΔDD[[2]]-ΔDD[[1]];
Print["c",ΔDD];*)

ΔDD=Chop[NInt[ dϵxceLDAdw2[nw]nw ,Quad]];
ΔDD={-w ΔDD,(1-w)ΔDD};
ΔΔDD=ΔDD[[2]]-ΔDD[[1]];

(* Individual energies *)
EeLDA+=ΔDD;

(* Excitation energies *)
ΩeLDA=EeLDA[[2]]-EeLDA[[1]];

(* Print results *)
TabRes=TableForm[{{
NumberForm[ENuc+EeLDA[[1]],{9,6}],NumberForm[ΔDD[[1]],{9,6}],NumberForm[ENuc+EeLDA[[2]],{9,6}],NumberForm[ΔDD[[2]],{9,6}],
NumberForm[ΩeLDA,{9,6}],NumberForm[ΔΔDD,{9,6}],NumberForm[ΩeLDA HaToeV,{9,6}],NumberForm[ΔΔDD HaToeV,{9,6}]
}},
TableHeadings->{None,{"(Subscript[E, eLDA]^((0)))","(Subscript[Δ, DD]^((0)))","(Subscript[E, eLDA]^((1)))","(Subscript[Δ, DD]^((1)))","(Subscript[Ω, eLDA]^((1)))","(Subscript[ΔΔ, DD]^((1)))","(Subscript[Ω, eLDA]^((1))) (eV)","(Subscript[ΔΔ, DD]^((1)))(eV)"}},TableAlignments->Center];
Print[TabRes];

(* Done *)
Return[];

];

Workhorse module

NInt[f_,Quad_]:=Which[
Quad=="Cartesian",NIntegrate[f,{x,-∞,∞},{y,-∞,∞},{z,-∞,∞} ],
Quad=="Cylindrical",NIntegrate[f r,{r,0,∞},{ϕ,0,2π},{z,-∞,∞}],
Quad=="Spherical",NIntegrate[f (r^2)Sin[θ],{r,0,∞},{ϕ,0,2π},{θ,0,π}]
];

Workhorse[n_,w_,XYZ_,ZNuc_,Quad_,FileInt_,Basis_]:=Module[
{xyz,NB,NAt,NShell,S,T,V,U,X,ENuc,BF,Atoms,ϵ,ZA,KA,αμ,Aμ,aμ,Dμ,χμ,Γ0,Γ1,nInt,TabInt,tCPU},

(* Start *)
Print[Style["Start of calculation",Red,Bold]];

SetDirectory["~/Dropbox/Manuscripts/FarDFT/Data/int"];

(* Numerical integration *)
xyz=Which[
Quad=="Cartesian",{x,y,z},
Quad=="Cylindrical",{r Cos[ϕ],r Sin[ϕ],z},
Quad=="Spherical",{r Cos[ϕ]Sin[θ],r Sin[ϕ]Sin[θ],r Cos[θ]}
];

(* General information *)
NAt=Length[XYZ];	NShell=Length[Basis];
BF=Flatten[Table[{Basis[[i,3]],XYZ[[Basis[[i,1]]]],AngQ[Basis[[i,2]]][[l]],Basis[[i,4]]},{i,NShell},{l,Length[AngQ[Basis[[i,2]]]]}],1];
NB=Length[BF];
Atoms=Table[{ZNuc[[i]],XYZ[[i]]},{i,NAt}];
ENuc=(∑_(c=1)^NAt)(∑_(d=c+1)^NAt)((ZNuc[[c]]ZNuc[[d]])/(Norm[XYZ[[c]]-XYZ[[d]]]));

(* Read one-electron integrals *)
Print["Reading one-electron integrals..."];
tCPU=Timing[
TabInt=ReadList[StringJoin["Ov.",FileInt,".dat"],{Number,Number,Real}];	nInt=Length[TabInt];
S = Normal[SparseArray[
Join[
Table[{TabInt[[k,1]],TabInt[[k,2]]}->TabInt[[k,3]],{k,nInt}],
Table[{TabInt[[k,2]],TabInt[[k,1]]}->TabInt[[k,3]],{k,nInt}]
]
]
];

TabInt=ReadList[StringJoin["Kin.",FileInt,".dat"],{Number,Number,Real}];	nInt=Length[TabInt];
T=SparseArray[
Join[
Table[{TabInt[[k,1]],TabInt[[k,2]]}->TabInt[[k,3]],{k,nInt}],
Table[{TabInt[[k,2]],TabInt[[k,1]]}->TabInt[[k,3]],{k,nInt}]
]
];

TabInt=ReadList[StringJoin["Nuc.",FileInt,".dat"],{Number,Number,Real}];	nInt=Length[TabInt];
V=SparseArray[
Join[
Table[{TabInt[[k,1]],TabInt[[k,2]]}->TabInt[[k,3]],{k,nInt}],
Table[{TabInt[[k,2]],TabInt[[k,1]]}->TabInt[[k,3]],{k,nInt}]
]
];
];
Print["Done!"];
Print["CPU time = ",tCPU[[1]]];

(* Compute orthognalization matrix *)
X=MatrixPower[S,-1/2];

(* Read two-electron integrals *)
Print["Reading two-electron integrals..."];
tCPU=Timing[
TabInt=ReadList[StringJoin["ERI.",FileInt,".dat"],{Number,Number,Number,Number,Real}];	nInt=Length[TabInt];
U= SparseArray[
Join[
Table[{TabInt[[k,1]],TabInt[[k,3]],TabInt[[k,2]],TabInt[[k,4]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,3]],TabInt[[k,1]],TabInt[[k,2]],TabInt[[k,4]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,1]],TabInt[[k,3]],TabInt[[k,4]],TabInt[[k,2]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,3]],TabInt[[k,1]],TabInt[[k,4]],TabInt[[k,2]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,2]],TabInt[[k,4]],TabInt[[k,1]],TabInt[[k,3]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,4]],TabInt[[k,2]],TabInt[[k,1]],TabInt[[k,3]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,2]],TabInt[[k,4]],TabInt[[k,3]],TabInt[[k,1]]}->TabInt[[k,5]],{k,nInt}],
Table[{TabInt[[k,4]],TabInt[[k,2]],TabInt[[k,3]],TabInt[[k,1]]}->TabInt[[k,5]],{k,nInt}]
]
];
];
Print["Done!"];
Print["CPU time = ",tCPU[[1]]];

(* Compute one-electron density *)
ZA=ZNuc;	
αμ=Table[BF[[μ,1]],{μ,NB}];	Aμ=Table[BF[[μ,2]],{μ,NB}];	aμ=Table[BF[[μ,3]],{μ,NB}];	Dμ=Table[BF[[μ,4]],{μ,NB}];
χμ=Table[(KA=Length[αμ[[μ]]];	Simplify[(∑_(kA=1)^KA)NormCoeff[αμ[[μ,kA]],aμ[[μ]]]Dμ[[μ,kA]]Apply[Times,(((xyz-Aμ[[μ]]))^(aμ[[μ]]))]Exp[-αμ[[μ,kA]]((Norm[xyz-Aμ[[μ]]])^2)],(x|y|z|r|θ|ϕ)∈Reals]),{μ,NB}];

(* Perform RHF calculation *)
Print["Ensemble RHF calculation..."];
{tCPU,{Γ0,Γ1}}=Timing[RHF[n,w,ENuc,NB,S,T,V,U,X]];
Print["Done!"];
Print["CPU time = ",tCPU];

(* Compute eDFT energies *)
Print["Compute GOK-DFT energies..."];
tCPU=Timing[GOKDFT[w,ENuc,T,V,U,Γ0,Γ1,χμ,Quad]];
Print["Done!"];
Print["CPU time = ",tCPU[[1]]];

(* End *)
Return[Style["End of calculation",Red,Bold]];
];

Examples

(Subscript[H, 2])

STO-3G

Workhorse[2,1/2,{{0.,0.,-0.7},{0.,0.,0.7}},{1,1},
"Cylindrical","H2_sto-3g",
{
{1,"S",{3.42525091,0.62391373,0.16885540},{0.15432897,0.53532814,0.44463454}},
{2,"S",{3.42525091,0.62391373,0.16885540},{0.15432897,0.53532814,0.44463454}}
}
]

4-31G

Workhorse[2,{{0.,0.,-0.7},{0.,0.,0.7}},{1,1},
{
{1,"S",{18.7311370,2.8253937,0.6401217},{0.03349460,0.23472695,0.81375733}},
{1,"S",{0.1612778},{1.}},
{2,"S",{18.7311370,2.8253937,0.6401217},{0.03349460,0.23472695,0.81375733}},
{2,"S",{0.1612778},{1.}}
}
]

6-31G**

Workhorse[2,{{0.,0.,-0.7},{0.,0.,0.7}},{1,1},
{
{1,"S",{18.7311370,2.8253937,0.6401217},{0.03349460,0.23472695,0.81375733}},
{1,"S",{0.1612778},{1.}},
{1,"P",{1.1000000},{1.}},
{2,"S",{18.7311370,2.8253937,0.6401217},{0.03349460,0.23472695,0.81375733}},
{2,"S",{0.1612778},{1.}},
{2,"P",{1.1000000},{1.}}
}]

cc-pVDZ

Workhorse[2,0,{{0.,0.,-0.7},{0.,0.,0.7}},{1.,1.},
"Cylindrical","H2_cc-pvdz",
{
{1,"S",{13.0100000,1.9620000,0.4446000},{0.0196850,0.1379770,0.4781480}},
{1,"S",{0.1220000},{1.}},
{1,"P",{0.7270000},{1.}},
{2,"S",{13.0100000,1.9620000,0.4446000},{0.0196850,0.1379770,0.4781480}},
{2,"S",{0.1220000},{1.}},
{2,"P",{0.7270000},{1.}}
}]

Workhorse[2,0,{{0.,0.,-0.7},{0.,0.,0.7}},{1.,1.},
"Cylindrical","H2_cc-pvdz",
{
{1,"S",{13.0100000,1.9620000,0.4446000},{0.0196850,0.1379770,0.4781480}},
{1,"S",{0.1220000},{1.}},
{1,"P",{0.7270000},{1.}},
{2,"S",{13.0100000,1.9620000,0.4446000},{0.0196850,0.1379770,0.4781480}},
{2,"S",{0.1220000},{1.}},
{2,"P",{0.7270000},{1.}}
}]

cc-pVTZ

Workhorse[2,1/2,{{0.,0.,-0.7},{0.,0.,0.7}},{1.,1.},
"Cylindrical","H2_cc-pvtz",
{
{1,"S",{3.387000 (10^1),5.095000 (10^0),1.159000 (10^1),3.258000 (10^(-1)),1.027000 (10^(-1))},{6.068000 (10^(-3)),4.530800 (10^(-1)),2.028220 (10^(-1)),5.039030 (10^(-1)),3.834210 (10^(-1))}},
{1,"S",{3.258000 (10^(-1))},{1.}},
{1,"S",{1.027000 (10^(-1))},{1.}},
{1,"P",{1.407000 (10^0)},{1.}},
{1,"P",{3.880000 (10^(-1))},{1.}},
{1,"D",{1.057000 (10^0)},{1.}},
{2,"S",{3.387000 (10^1),5.095000 (10^0),1.159000 (10^1),3.258000 (10^(-1)),1.027000 (10^(-1))},{6.068000 (10^(-3)),4.530800 (10^(-1)),2.028220 (10^(-1)),5.039030 (10^(-1)),3.834210 (10^(-1))}},
{2,"S",{3.258000 (10^(-1))},{1.}},
{2,"S",{1.027000 (10^(-1))},{1.}},
{2,"P",{1.407000 (10^0)},{1.}},
{2,"P",{3.880000 (10^(-1))},{1.}},
{2,"D",{1.057000 (10^0)},{1.}}
}]

Be

STO-3G

Workhorse[4,0,{{0.,0.,0.}},{4},
"Spherical","Be_sto-3g",
{
{1,"S",{30.1678710,5.4951153,1.4871927},{0.15432897,0.53532814,0.44463454}},
{1,"S",{1.3148331,0.3055389,0.0993707},{-0.0999672,0.3995128,0.7001155}},
{1,"P",{1.3148331,0.3055389,0.0993707},{0.1559163,0.6076837,0.3919574}}
}]

6-31+G*

Workhorse[4,0,{{0.,0.,0.}},{4},
"Spherical","Be_6-31+g_star",
{
{1,"S",{0.1264585690 (10^4),0.1899368060 (10^3),0.4315908900 (10^2),0.1209866270 (10^2),0.3806323220 (10^1),0.1272890300 (10^1)},{0.1944757590 (10^(-2)),0.1483505200 (10^(-1)),0.7209054629 (10^(-1)),0.2371541500,0.4691986519,0.3565202279}},
{1,"S",{0.3196463098 (10^1),0.7478133038,0.2199663302},{-0.1126487285,-0.2295064079,0.1186916764 (10^1)}},
{1,"P",{0.3196463098 (10^1),0.7478133038,0.2199663302},{0.5598019980 (10^(-1)),0.2615506110,0.7939723389}},
{1,"S",{0.8230990070 (10^(-1))},{1.}},
{1,"P",{0.8230990070 (10^(-1))},{1.}},
{1,"S",{0.2070000000 (10^(-1))},{1.}},
{1,"P",{0.2070000000 (10^(-1))},{1.}},
{1,"D",{0.4000000000},{1.}}
}]

GOK-DFT made simple

With an exchange-only functional (i.e., without correlation), the ensemble energy is

(E^w)=(Subscript[T, s]^w)+(Subscript[E, H]^w)+(Subscript[E, x]^w)

where, in particular, the Hartree and LDA exchange energy reads

(Subscript[E, H]^w)=(1/2)∫∫(((n^w)[r](n^w)[r'])/(\[LeftBracketingBar]r-r'\[RightBracketingBar]))\[DifferentialD]r\[DifferentialD]r'	(Subscript[E, x]^w)=(Subscript[C, x])∫((((n^w)[r]))^(4/3))\[DifferentialD]r

and (n^w)[r] is the ensemble density. This energy is obviously not linear wrt w due to the approximate nature of the functionals.

Our aim is to define a new (weight-dependent) exchange functional

(Subscript[(OverTilde[E]), x]^w)=(Subscript[(OverTilde[C]), x])∫((((n^w)[r]))^(4/3))\[DifferentialD]r

such that the ensemble energy

((OverTilde[E])^w)=(Subscript[(OverTilde[T]), s]^w)+(Subscript[(OverTilde[E]), H]^w)+(Subscript[(OverTilde[E]), x]^w)

associated with the ensemble density ((OverTilde[n])^w)[r] is linear between w=0 and w=1, i.e.,

((OverTilde[E])^w)= ((E^(w=1))-(E^(w=0)))w+(E^(w=0))

Assuming that the ensemble densities obtained by the minimization of (E^w) and ((OverTilde[E])^w)are equal (which might be a strong approximation in certain cases), i.e., (n^w)[r]=((OverTilde[n])^w)[r], the non-linearity of the energy can be approximated as

Δ=(E^w)-((OverTilde[E])^w)≈(Subscript[E, x]^w)-(Subscript[(OverTilde[E]), x]^w)	\[Implies]	(Subscript[(OverTilde[E]), x]^w)=(Subscript[E, x]^w)+((OverTilde[E])^w)-(E^w)

because the kinetic and Hartree energies depend only on the ensemble density (i.e., they are not explicitly weight dependent).
Assuming a linear-in-w behavior of (Subscript[(OverTilde[C]), x]) in the vicinity of w=0 (which is actually not an assumption I believe), i.e.,

(Subscript[(OverTilde[C]), x])=(Subscript[C, x])(1+α w)

we have

(Subscript[(OverTilde[E]), x]^w)=(1+α w)(Subscript[E, x]^w)

yielding

(Subscript[E, x]^w)(1+α w)=(Subscript[E, x]^w)+((OverTilde[E])^w)-(E^w)	\[Implies]	α=((((OverTilde[E])^w)-(E^w))/(w (Subscript[E, x]^w)))

Now, let’s define an infinitesimal weight ϵ>0. Evaluating the latest expression at w=ϵ, we get

α=((((OverTilde[E])^(w=ϵ))-(E^(w=ϵ)))/( ϵ (Subscript[E, x]^(w=ϵ))))=((((E^(w=1))-(E^(w=0)))ϵ+(E^(w=0))-(E^(w=ϵ)))/(ϵ (Subscript[E, x]^(w=ϵ))))=(( ϵ (E^(w=1))-(1-ϵ)(E^(w=0))-(E^(w=ϵ)))/(ϵ (Subscript[E, x]^(w=ϵ))))

This is our final expression.

Because

(Subscript[((∂(Subscript[(OverTilde[E]), x]^w)[n])/(∂w))\[RightBracketingBar], n=(n^w)])=(∂/(∂w))((Subscript[C, x])∫(1+α w)((n[r])^(4/3))\[DifferentialD]r)=α (Subscript[E, x]^w)[n]

((∂(Subscript[(OverTilde[E]), x]^w))/(∂w))=α (Subscript[E, x]^(w=ϵ))=((ϵ (E^(w=1))-(1-ϵ)(E^(w=0))-(E^(w=ϵ)))/ϵ)

the previous expression is an estimate of the exchange part of the ensemble derivative at the vicinity of w=0.

Summarizing, in practice, one needs to perform:

a zero-weight (i.e., ground state) calculation to get (E^(w=0)). Easy.

a MOM-style calculation at w=1 to get (E^(w=1)). Less easy but doable with the new algorithms developed in Gill’s group and Head-Gordon’s group (MOM, SGM, etc).

a GOK-DFT calculation at w=ϵ to get (E^(w=ϵ)) and (Subscript[E, x]^(w=ϵ)). Because ϵ is small, in practice, the calculation is going to converge easily, which is a very important point, especially if one starts with the converged ground-state density matrix.

In such a way, one can approximate efficiently the ensemble derivative at w=0, and obtains excitation energies.

What do you guys think?

Results

I=0:	2H	1-(Subscript[w, 1])-(Subscript[w, 2])
I=1:	H+L      (Subscript[w, 1])
I=2:	2L	(Subscript[w, 2])

(1-(Subscript[w, 1])-(Subscript[w, 2]))2H+(Subscript[w, 1])(H+L)+(Subscript[w, 2])(2L)==(1-((Subscript[w, 1])/2)-(Subscript[w, 2]))2H+(((Subscript[w, 1])/2)+(Subscript[w, 2]))2L

Hydrogen molecule at equilibrium

SizeLegend=16;
SizeLabel=16;
OptionsPlot={
PlotTheme->"Detailed",
BaseStyle->16,
ImageSize->Medium,
PlotStyle->Thickness[0.01],
InterpolationOrder->2,
Joined->True
(*PlotMarkers->"OpenMarkers"*)
};

SVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RS51_RVWN5.dat"];
GIC=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RGIC_H.dat"];
GICVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RGIC_RVWN5.dat"];
GICMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RGIC_RMFL20.dat"];
S51=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RS51_H.dat"];
MFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RMFL20_H.dat"];
MFLVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RMFL20_RVWN5.dat"];
MFLMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RMFL20_RMFL20.dat"];
ExS51=Table[{S51[[k,1]],S51[[k,4]]},{k,Length[S51]}];
EwS51=Table[{S51[[k,1]],S51[[k,2]]},{k,Length[S51]}];
ΩS51=Table[{S51[[k,1]],S51[[k,3]]},{k,Length[S51]}];
ΔΩS51=Table[{S51[[k,1]],S51[[k,3]]HaToeV-28.75},{k,Length[S51]}];
EwSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,2]]},{k,Length[SVWN5]}];
ΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]},{k,Length[SVWN5]}];
ΔΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]HaToeV-28.75},{k,Length[SVWN5]}];
EwGIC=Table[{GIC[[k,1]],GIC[[k,2]]},{k,Length[GIC]}];
ΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]},{k,Length[GIC]}];
ΔΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]HaToeV-28.75},{k,Length[GIC]}];
EwGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,2]]},{k,Length[GICVWN5]}];
ΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]},{k,Length[GICVWN5]}];
ΔΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]HaToeV-28.75},{k,Length[GICVWN5]}];
EwGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,2]]},{k,Length[GICMFL]}];
ΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]},{k,Length[GICMFL]}];
ΔΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]HaToeV-28.75},{k,Length[GICMFL]}];
FCI={{0,28.75},{1/2,28.75}};

ListPlot[{
ΔEw[EwS51],ΔEw[EwGIC],ΔEw[EwSVWN5],ΔEw[EwGICVWN5],ΔEw[EwGICMFL]
},
PlotRange->{{0,0.5},Automatic},
OptionsPlot,
FrameLabel->{
MaTeX["\\text{$\\xi = w_1/2+w_2 \\land w_1 = w_2$}",FontSize->SizeLabel],
MaTeX["\\text{Dev. from linearity (a.u.)}",FontSize->SizeLabel]
},
PlotLegends->{
MaTeX["E^\\xi_\\text{S}",FontSize->SizeLegend],
MaTeX["E^\\xi_\\text{CC-S}",FontSize->SizeLegend],
MaTeX["E^\\xi_\\text{SVWN5}",FontSize->SizeLegend],
MaTeX["E^\\xi_\\text{CC-SVWN5}",FontSize->SizeLegend],
MaTeX["E^\\xi_\\text{CC-SeVWN5}",FontSize->SizeLegend]
}
]
Export["~/Dropbox/Manuscripts/FarDFT/Manuscript/Ew_H2.pdf",%];

ListPlot[{
ΔΩS51,ΔΩGIC,ΔΩSVWN5,ΔΩGICVWN5,ΔΩGICMFL
},
OptionsPlot,PlotRange->{All,All},
FrameLabel->{
MaTeX["\\text{$W$}",FontSize->SizeLabel],
MaTeX["\\text{Error in excitation energy (eV)}",FontSize->SizeLabel]
},
PlotLegends->{
MaTeX["\\Omega_\\text{S}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{CC-S}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{SVWN5}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{CC-SVWN5}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{CC-SeVWN5}",FontSize->SizeLegend]
}
]
Export["~/Dropbox/Manuscripts/FarDFT/Manuscript/Om_H2.pdf",%];

Plot[{
1+w(w-1)(a+b (w-1/2)+c (((w-1/2))^2))/.{a->1.9015719148496788`,b->2.5236598782764412`,c->1.6652282199359842`},
1+w(w-1)(a+b (w-1/2)+c (((w-1/2))^2))/.{a->0.5739189000851961`,b->-0.0003469882157336496`,c->-0.2676338054343272`},
1+w(w-1)(a+b (w-1/2)+c (((w-1/2))^2))/.{a->0.01918229168254928`,b->-0.01545313842512261`,c->-0.012720073519142448`}
},{w,0,1/2},
PlotStyle->{{Thickness[0.01],Blue},{Thickness[0.01],Red},{Thickness[0.01],Darker[Green]}},PlotTheme->"Detailed",
BaseStyle->16,AspectRatio->3/2,
ImageSize->Medium,
PlotLegends->{
MaTeX["\\ce{He}",FontSize->SizeLegend],
MaTeX["\\text{\\ce{H2} with $R_{\\ce{H-H}} = 1.4$ bohr}",FontSize->SizeLegend],
MaTeX["\\text{\\ce{H2} with $R_{\\ce{H-H}} = 3.7$ bohr}",FontSize->SizeLegend]
},
FrameLabel->{
MaTeX["\\text{$w_1/2 + w_2$}",FontSize->SizeLabel],
MaTeX["C_\\text{x}^\\mathbold{w}/C_\\text{x}",FontSize->SizeLabel]
}]
Export["~/Dropbox/Manuscripts/FarDFT/Manuscript/Cxw.pdf",%];

EwGICS51=Table[{w,(S51[[41,2]]-S51[[1,2]])w+S51[[1,2]]},{w,0,1/2,0.025}];

FindFit[Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,21}],w(w-1)(a+b (w-1/2)+c (((w-1/2))^2)),{a,b,c},w,MaxIterations->1000]
(*EwGICSVWN5=Table[{w,(SVWN5[[41,2]]-SVWN5[[1,2]])w+SVWN5[[1,2]]},{w,0,1,0.025}];
FindFit[Table[{SVWN5[[k,1]],((SVWN5[[k,2]]-EwGICSVWN5[[k,2]])/(SVWN5[[k,2]]))},{k,Length[EwS51]}],w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),{a,b,c},w,MaxIterations->1000]*)

Show[{
ListPlot[
{
Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,Length[EwS51]}]
}
,
OptionsPlot,
PlotRange->{{0,1/2},Automatic},
FrameLabel->{
MaTeX["\\text{Weight $w$}",FontSize->SizeLabel],
MaTeX["\\text{Excitation energy (eV)}",FontSize->SizeLabel]
}
,
PlotLegends->Automatic
]
,
Plot[{
w(w-1)(a+b (w-1/2)+c (((w-1/2))^2))/.{a->0.5739189000851961`,b->-0.0003469882157336496`,c->-0.2676338054343272`}
},{w,0,1},PlotStyle->Red]
}]

Series[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.5751782560799208`,b->-0.021108186591137282`,c->-0.36718902716347124`},{w,0,1}]

D[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),w]==b/2+(2a+c/2) (w-1/2)-(1-w) w(3 b+4c ( w-1/2) )//FullSimplify

D[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),w]==(b/2+(2a+c/2)(w-1/2) - (1-w)w(3b+4c(w-1/2)))//Simplify

TableForm[
{
{"HF","","aug-cc-pVDZ",38.52,35.04,30.86,34.55,28.65},
{"HF","","aug-cc-pVTZ",38.58,35.77,35.82,35.68,28.65},
{"HF","","aug-cc-pVQZ",39.12,37.70,35.94,35.64,28.65},
{},
{"HF","VWN5","aug-cc-pVDZ",40.12,37.97,31.19,35.66,29.17},
{"HF","VWN5","aug-cc-pVTZ",40.84,38.88,36.98,37.23,29.17},
{"HF","VWN5","aug-cc-pVQZ",41.39,,37.07,37.21,29.17},
{},
{"HF","MFL20","aug-cc-pVDZ",40.37,,31.19,36.98,30.51},
{"HF","MFL20","aug-cc-pVTZ",41.09,,36.56,38.35,30.51},
{"HF","MFL20","aug-cc-pVQZ",41.64,,36.66,38.33,30.51},
{},
{"S","","aug-cc-pVDZ",38.40,,27.35,23.54,26.60},
{"S","","aug-cc-pVTZ",39.21,,27.42,23.62,26.67},
{"S","","aug-cc-pVQZ",39.78,,27.42,23.62,26.67},
{},
{"S","WN5","aug-cc-pVDZ",38.10,32.52,27.76,24.40,27.10},
{"S","WN5","aug-cc-pVTZ",38.54,32.68,27.81,24.46,27.17},
{"S","WN5","aug-cc-pVQZ",38.81,,27.81,24.46,27.17},
{},
{"S","MFL20","aug-cc-pVDZ",38.34,,27.74,25.80,28.42},
{"S","MFL20","aug-cc-pVTZ",38.79,32.85,27.80,25.87,28.49},
{"S","MFL20","aug-cc-pVQZ",39.06,,27.80,25.87,28.49},
{},
{"MFL20","","aug-cc-pVDZ",33.65,,24.69,20.12,23.45},
{"MFL20","","aug-cc-pVTZ",34.46,,24.73,20.17,23.52},
{"MFL20","","aug-cc-pVQZ",35.02,,24.72,20.17,23.52},
{},
{"MFL20","VWN5","aug-cc-pVDZ",33.28,,24.97,20.87,23.85},
{"MFL20","VWN5","aug-cc-pVTZ",33.72,28.93,25.00,20.91,23.91},
{"MFL20","VWN5","aug-cc-pVQZ",33.99,,24.99,20.90,23.91},
{},
{"MFL20","MFL20","aug-cc-pVDZ",16.47,21.34,25.14,20.96,24.04},
{"MFL20","MFL20","aug-cc-pVTZ",,,,,},
{"MFL20","MFL20","aug-cc-pVQZ",,,,,},
{},
{"B88","LYP","aug-mcc-pV8Z","","","","",28.42},
{"B3","LYP","aug-mcc-pV8Z","","","","",27.77},
{"HF","LYP","aug-mcc-pV8Z","","","","",29.18},
{"HF","","aug-mcc-pV8Z","","","","",28.65},
{},
{"HF","FCI","aug-mcc-pV8Z","","","","",28.75}
}
,TableHeadings->{None,{
"x","c","Basis",
"GOK(w=0)","GOK(w=1/4)","GOK(w=1/2)","LIM","MOM"}}]

Hydrogen molecule stretched

SizeLegend=16;
SizeLabel=16;
OptionsPlot={
PlotTheme->"Detailed",
BaseStyle->16,
ImageSize->Medium,
PlotStyle->Thickness[0.01],
InterpolationOrder->2,
Joined->True
(*PlotMarkers->"OpenMarkers"*)
};

SVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RS51_RVWN5.dat"];
GIC=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2st_aug-cc-pvtz_RGIC_H.dat"];
GICVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RGIC_RVWN5.dat"];
GICMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RGIC_RMFL20.dat"];
S51=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2st_aug-cc-pvtz_RS51_H.dat"];
MFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2st_aug-cc-pvtz_RMFL20_H.dat"];
MFLVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RMFL20_RVWN5.dat"];
MFLMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/H2_aug-cc-pvtz_RMFL20_RMFL20.dat"];
ExS51=Table[{S51[[k,1]],S51[[k,4]]},{k,Length[S51]}];
EwS51=Table[{S51[[k,1]],S51[[k,2]]},{k,Length[S51]}];
ΩS51=Table[{S51[[k,1]],S51[[k,3]]},{k,Length[S51]}];
ΔΩS51=Table[{S51[[k,1]],S51[[k,3]]HaToeV-8.69},{k,Length[S51]}];
EwSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,2]]},{k,Length[SVWN5]}];
ΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]},{k,Length[SVWN5]}];
ΔΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]HaToeV-8.69},{k,Length[SVWN5]}];
EwGIC=Table[{GIC[[k,1]],GIC[[k,2]]},{k,Length[GIC]}];
ΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]},{k,Length[GIC]}];
ΔΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]HaToeV-8.69},{k,Length[GIC]}];
EwGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,2]]},{k,Length[GICVWN5]}];
ΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]},{k,Length[GICVWN5]}];
ΔΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]HaToeV-8.69},{k,Length[GICVWN5]}];
EwGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,2]]},{k,Length[GICMFL]}];
ΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]},{k,Length[GICMFL]}];
ΔΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]HaToeV-8.69},{k,Length[GICMFL]}];
FCI={{0,8.69},{1/2,8.69}};

ΔEw[tab_]:=Table[{tab[[k,1]],tab[[k,2]]-(tab[[Length[tab],2]]-tab[[1,2]])tab[[k,1]]-tab[[1,2]]},{k,Length[tab]}]

ListPlot[{
ΔEw[EwS51],ΔEw[EwGIC](*,ΔEw[EwSVWN5],ΔEw[EwGICVWN5],ΔEw[EwGICMFL]*)
},
OptionsPlot,
FrameLabel->{
MaTeX["\\text{$w$}",FontSize->SizeLabel],
MaTeX["\\text{Dev. from linearity (a.u.)}",FontSize->SizeLabel]
},
PlotLegends->{
MaTeX["E^w_\\text{S}",FontSize->SizeLegend],
MaTeX["E^w_\\text{GIC-S}",FontSize->SizeLegend],
MaTeX["E^w_\\text{SVWN5}",FontSize->SizeLegend],
MaTeX["E^w_\\text{GIC-SVWN5}",FontSize->SizeLegend],
MaTeX["E^w_\\text{GIC-SeVWN5}",FontSize->SizeLegend]
}
]

ListPlot[{
ΔΩS51,ΔΩGIC(*,ΔΩSVWN5,ΔΩGICVWN5,ΔΩGICMFL*)
},
OptionsPlot,PlotRange->{All,All},
FrameLabel->{
MaTeX["\\text{$w$}",FontSize->SizeLabel],
MaTeX["\\text{Error in excitation energy (eV)}",FontSize->SizeLabel]
},
PlotLegends->{
MaTeX["\\Omega_\\text{S}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{GIC-S}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{SVWN5}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{GIC-SVWN5}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{GIC-SeVWN5}",FontSize->SizeLegend]
}
]

Plot[{
1+w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.5751782560799208`,b->-0.021108186591137282`,c->-0.36718902716347124`}
},{w,0,1},
PlotStyle->{Thickness[0.01],Red},PlotTheme->"Detailed",
BaseStyle->16,
ImageSize->Medium,PlotLegends->None,
FrameLabel->{
MaTeX["\\text{$w$}",FontSize->SizeLabel],
MaTeX["C_\\text{x}^w/C_\\text{x}^{w=0}",FontSize->SizeLabel]
}]

EwGICS51=Table[{w,(S51[[41,2]]-S51[[1,2]])w+S51[[1,2]]},{w,0,1/2,0.025}];

FindFit[Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,21}],w(w-1)(a+b (w-1/2)+c (((w-1/2))^2)),{a,b,c},w,MaxIterations->1000]
(*EwGICSVWN5=Table[{w,(SVWN5[[41,2]]-SVWN5[[1,2]])w+SVWN5[[1,2]]},{w,0,1,0.025}];
FindFit[Table[{SVWN5[[k,1]],((SVWN5[[k,2]]-EwGICSVWN5[[k,2]])/(SVWN5[[k,2]]))},{k,Length[EwS51]}],w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),{a,b,c},w,MaxIterations->1000]*)

(-a+(b/2)-(c/4)) w+((O[w])^2)/.{a->0.5094314403492467`,b->0.9182410165442554`,c->0.8477210543586113`}

Show[{
ListPlot[
{
Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,Length[EwS51]}]
}
,
OptionsPlot,
PlotRange->{{0,1},Automatic},
FrameLabel->{
MaTeX["\\text{Weight $w$}",FontSize->SizeLabel],
MaTeX["\\text{Excitation energy (eV)}",FontSize->SizeLabel]
}
,
PlotLegends->Automatic
]
,
Plot[{
w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.5094314403492467`,b->0.9182410165442554`,c->0.8477210543586113`}|
},{w,0,1},PlotStyle->Red]
}]

Series[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.5751782560799208`,b->-0.021108186591137282`,c->-0.36718902716347124`},{w,0,1}]

D[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),w]==b/2+(2a+c/2) (w-1/2)-(1-w) w(3 b+4c ( w-1/2) )//FullSimplify

D[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),w]==(b/2+(2a+c/2)(w-1/2) - (1-w)w(3b+4c(w-1/2)))//Simplify

Helium atom

SizeLegend=20;
SizeLabel=20;
OptionsPlot={
PlotTheme->"Detailed",
BaseStyle->20,
ImageSize->Large,
Joined->True,
PlotMarkers->"OpenMarkers"
};

SVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RS51_RVWN5.dat"];
GIC=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RGIC_H.dat"];
GICVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RGIC_RVWN5.dat"];
GICMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RGIC_RMFL20.dat"];
S51=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RS51_H.dat"];
MFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RMFL20_H.dat"];
MFLVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RMFL20_RVWN5.dat"];
MFLMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/He_aug-cc-pvtz_RMFL20_RMFL20.dat"];
ExS51=Table[{S51[[k,1]],S51[[k,4]]},{k,Length[S51]}];
EwS51=Table[{S51[[k,1]],S51[[k,2]]},{k,Length[S51]}];
ΩS51=Table[{S51[[k,1]],S51[[k,3]]},{k,Length[S51]}];
ΔΩS51=Table[{S51[[k,1]],S51[[k,3]]-2.126 },{k,Length[S51]}];
EwSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,2]]},{k,Length[SVWN5]}];
ΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]},{k,Length[SVWN5]}];
ΔΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]-2.126 },{k,Length[SVWN5]}];
EwGIC=Table[{GIC[[k,1]],GIC[[k,2]]},{k,Length[GIC]}];
ΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]},{k,Length[GIC]}];
ΔΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]-2.126 },{k,Length[GIC]}];
EwGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,2]]},{k,Length[GICVWN5]}];
ΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]},{k,Length[GICVWN5]}];
ΔΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]-2.126 },{k,Length[GICVWN5]}];
EwGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,2]]},{k,Length[GICMFL]}];
ΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]},{k,Length[GICMFL]}];
ΔΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]-2.126},{k,Length[GICMFL]}];
FCI={{0,2.126},{1/2,2.126}};

ListPlot[{
EwS51,EwGIC,EwGICVWN5,EwGICMFL
},
OptionsPlot,
FrameLabel->{
MaTeX["\\text{Weight $w$}",FontSize->SizeLabel],
MaTeX["\\text{Ensemble energy (au)}",FontSize->SizeLabel]
},
PlotLegends->Automatic
]

ListPlot[{
ΔΩS51,ΔΩGIC,ΔΩGICVWN5,ΔΩGICMFL
},
OptionsPlot,
FrameLabel->{
MaTeX["\\text{Weight $w$}",FontSize->SizeLabel],
MaTeX["\\text{Excitation energy (au)}",FontSize->SizeLabel]
},
PlotLegends->Automatic
]

EwGICS51=Table[{w,(S51[[41,2]]-S51[[1,2]])w+S51[[1,2]]},{w,0,1,0.025}];

Table[((EwS51[[k+1,2]]-EwS51[[k-1,2]])/(2 0.025))HaToeV,{k,2,Length[EwS51]-1}]
Table[((EwGICS51[[k+1,2]]-EwGICS51[[k-1,2]])/(2 0.025))HaToeV,{k,2,Length[EwGICS51]-1}]

EwGICS51=Table[{w,(S51[[41,2]]-S51[[1,2]])w+S51[[1,2]]},{w,0,1/2,0.025}];

FindFit[Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,21}],w(w-1)(a+b (w-1/2)+c (((w-1/2))^2)),{a,b,c},w,MaxIterations->1000]

Show[{
ListPlot[
{
Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,Length[EwS51]}]
}
,
OptionsPlot,
PlotRange->{{0,1},Automatic},
FrameLabel->{
MaTeX["\\text{Weight $w$}",FontSize->SizeLabel],
MaTeX["\\text{Excitation energy (eV)}",FontSize->SizeLabel]
}
,
PlotLegends->Automatic
]
,
Plot[{
w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->1.9125735895875828`,b->2.715266992840757`,c->2.1634223380633086`}
},{w,0,1},PlotStyle->Red]
}]

Series[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->1.9125735895875828`,b->2.715266992840757`,c->2.1634223380633086`},{w,0,2}]
Series[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.5751782560799208`,b->-0.021108186591137282`,c->-0.36718902716347124`},{w,0,2}]

HaToeV(∫_0^(1/2))w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))\[DifferentialD]w/.{a->1.9125735895875828`,b->2.715266992840757`,c->2.1634223380633086`}
HaToeV(∫_0^(1/2))w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))\[DifferentialD]w/.{a->0.5751782560799208`,b->-0.021108186591137282`,c->-0.36718902716347124`}

27.35-1.2716256877926888`

28.65

38.52-2.525301906413077`

Beryllium atom

TableForm[
{
{"GOK","HF","","aug-cc-pVDZ",17.25,8.23,14.00,8.29},
{"GOK","HF","","aug-cc-pVTZ",17.04,8.17,13.91,8.22},
{"GOK","HF","","aug-cc-pVQZ",26.17,8.15,13.90,8.20},
{},
{"GOK","HF","VWN5","aug-cc-pVDZ",,,,},
{"GOK","HF","VWN5","aug-cc-pVTZ",,,,},
{"GOK","HF","VWN5","aug-cc-pVQZ",,,,},
{},
{"GOK","HF","MFL20","aug-cc-pVDZ",16.20,8.69,18.79,12.89},
{"GOK","HF","MFL20","aug-cc-pVTZ",,,,},
{"GOK","HF","MFL20","aug-cc-pVQZ",,,,},
{},
{"GOK","S","","aug-cc-pVDZ",,,,},
{"GOK","S","","aug-cc-pVTZ",,,,},
{"GOK","S","","aug-cc-pVQZ",,,,},
{},
{"GOK","S","VWN5","aug-cc-pVDZ",8.57,7.15,7.11,7.87,7.07},
{"GOK","S","VWN5","aug-cc-pVTZ",8.49,7.06,7.04,7.77,6.97},
{"GOK","S","VWN5","aug-cc-pVQZ",8.48,7.04,7.02,7.75},
{},
{"GOK","S","MFL20","aug-cc-pVDZ",9.12,,,},
{"GOK","S","MFL20","aug-cc-pVTZ",9.04,,,},
{"GOK","S","MFL20","aug-cc-pVQZ",,,,},
{},
{"GOK","MFL20","","aug-cc-pVDZ",,,,},
{"GOK","MFL20","","aug-cc-pVTZ",,,,},
{"GOK","MFL20","","aug-cc-pVQZ",,,,},
{},
{"GOK","MFL20","VWN5","aug-cc-pVDZ",,,,},
{"GOK","MFL20","VWN5","aug-cc-pVTZ",,,,},
{"GOK","MFL20","VWN5","aug-cc-pVQZ",,,,},
{},
{"GOK","MFL20","MFL20","aug-cc-pVDZ",,,,,8.12},
{"GOK","MFL20","MFL20","aug-cc-pVTZ",,,,,},
{"GOK","MFL20","MFL20","aug-cc-pVQZ",,,,,},
{},
{"MOM","S","PW92","aug-cc-pVTZ",6.97},
{"MOM","PBE","PBE","aug-cc-pVTZ",6.98},
{"MOM","SCAN","SCAN","aug-cc-pVTZ",7.11},
{"MOM","PBE0","PBE0","aug-cc-pVTZ",7.23},
{"FCI","","","aug-cc-pVTZ",7.15},
{},
{"FCI","","aug-cc-pVQZ",7.11},
{"TBE","","",7.06}
}
,TableHeadings->{None,{
"Method","Exchange","Correlation","Basis",
"w=0","w=1/2","LIM","GIC-LIM","MOM"}}]

HNO

SizeLegend=16;
SizeLabel=16;
OptionsPlot={
PlotTheme->"Detailed",
BaseStyle->16,
ImageSize->Medium,
PlotStyle->Thickness[0.01],
InterpolationOrder->2,
Joined->True
(*PlotMarkers->"OpenMarkers"*)
};

SVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO2_aug-cc-pvtz_RS51_RVWN5.dat"];
GIC=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug-cc-pvtz_RGIC_H.dat"];
GICVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug-cc-pvtz_RGIC_RVWN5.dat"];
GICMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug-cc-pvtz_RGIC_RMFL20.dat"];
S51=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug_cc-pvdz_RS51_H.dat"];
MFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug-cc-pvtz_RMFL20_H.dat"];
MFLVWN5=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug-cc-pvtz_RMFL20_RVWN5.dat"];
MFLMFL=Import["~/Dropbox/Manuscripts/FarDFT/Data/HNO_aug-cc-pvdz_RMFL20_RMFL20.dat"];
ExS51=Table[{S51[[k,1]],S51[[k,4]]},{k,Length[S51]}];
EwS51=Table[{S51[[k,1]],S51[[k,2]]},{k,Length[S51]}];
ΩS51=Table[{S51[[k,1]],S51[[k,3]]},{k,Length[S51]}];
ΔΩS51=Table[{S51[[k,1]],S51[[k,3]]HaToeV-4.32},{k,Length[S51]}];
EwSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,2]]},{k,Length[SVWN5]}];
ΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]},{k,Length[SVWN5]}];
ΔΩSVWN5=Table[{SVWN5[[k,1]],SVWN5[[k,3]]HaToeV-4.32},{k,Length[SVWN5]}];
EwGIC=Table[{GIC[[k,1]],GIC[[k,2]]},{k,Length[GIC]}];
ΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]},{k,Length[GIC]}];
ΔΩGIC=Table[{GIC[[k,1]],GIC[[k,3]]HaToeV-4.32},{k,Length[GIC]}];
EwGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,2]]},{k,Length[GICVWN5]}];
ΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]},{k,Length[GICVWN5]}];
ΔΩGICVWN5=Table[{GICVWN5[[k,1]],GICVWN5[[k,3]]HaToeV-4.32},{k,Length[GICVWN5]}];
EwGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,2]]},{k,Length[GICMFL]}];
ΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]},{k,Length[GICMFL]}];
ΔΩGICMFL=Table[{GICMFL[[k,1]],GICMFL[[k,3]]HaToeV-4.32},{k,Length[GICMFL]}];
FCI={{0,4.32},{1/2,4.32}};

ΔEw[tab_]:=Table[{tab[[k,1]],tab[[k,2]]-(tab[[Length[tab],2]]-tab[[1,2]])tab[[k,1]]-tab[[1,2]]},{k,Length[tab]}]

ListPlot[{
ΔEw[EwS51](*,ΔEw[EwGIC],ΔEw[EwSVWN5],ΔEw[EwGICVWN5],ΔEw[EwGICMFL]*)
},
OptionsPlot,
FrameLabel->{
MaTeX["\\text{$w$}",FontSize->SizeLabel],
MaTeX["\\text{Dev. from linearity (a.u.)}",FontSize->SizeLabel]
},
PlotLegends->{
MaTeX["E^w_\\text{S}",FontSize->SizeLegend],
MaTeX["E^w_\\text{GIC-S}",FontSize->SizeLegend],
MaTeX["E^w_\\text{SVWN5}",FontSize->SizeLegend],
MaTeX["E^w_\\text{GIC-SVWN5}",FontSize->SizeLegend],
MaTeX["E^w_\\text{GIC-SeVWN5}",FontSize->SizeLegend]
}
]

ListPlot[{
ΔΩS51(*,ΔΩGIC,ΔΩSVWN5,ΔΩGICVWN5,ΔΩGICMFL*)
},
OptionsPlot,PlotRange->{All,All},
FrameLabel->{
MaTeX["\\text{$w$}",FontSize->SizeLabel],
MaTeX["\\text{Error in excitation energy (eV)}",FontSize->SizeLabel]
},
PlotLegends->{
MaTeX["\\Omega_\\text{S}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{GIC-S}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{SVWN5}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{GIC-SVWN5}",FontSize->SizeLegend],
MaTeX["\\Omega_\\text{GIC-SeVWN5}",FontSize->SizeLegend]
}
]

(∫_0^1)w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))\[DifferentialD]w/.{a->0.0061158387543040335`,b->-0.00005968703047293955`,c->-0.00001692245714408755`}

Plot[{
1+w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.0061158387543040335`,b->-0.00005968703047293955`,c->-0.00001692245714408755`}
},{w,0,1},
PlotStyle->{Thickness[0.01],Red},PlotTheme->"Detailed",
BaseStyle->16,
ImageSize->Medium,PlotLegends->None,
FrameLabel->{
MaTeX["\\text{$w$}",FontSize->SizeLabel],
MaTeX["C_\\text{x}^w/C_\\text{x}^{w=0}",FontSize->SizeLabel]
}]

EwGICS51=Table[{w,(S51[[41,2]]-S51[[1,2]])w+S51[[1,2]]},{w,0,1,0.025}];

FindFit[Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,Length[EwS51]}],w(w-1)(a+b (w-1/2)+c(((w-1/2))^2)),{a,b,c},w,MaxIterations->1000]

Show[{
ListPlot[
{
Table[{EwS51[[k,1]],((EwGICS51[[k,2]]-EwS51[[k,2]])/(ExS51[[k,2]]))},{k,Length[EwS51]}]
}
,
OptionsPlot,PlotMarkers->"OpenMarkers",
PlotRange->{{0,1},Automatic},
FrameLabel->{
MaTeX["\\text{Weight $w$}",FontSize->SizeLabel],
MaTeX["\\text{Excitation energy (eV)}",FontSize->SizeLabel]
}
,
PlotLegends->Automatic
]
,
Plot[{
w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.0061158387543040335`,b->-0.00005968703047293955`,c->-0.00001692245714408755`}
},{w,0,1},PlotStyle->Red]
}]

Series[w(w-1)(a+b (w-1/2)+c(((w-1/2))^2))/.{a->0.0061158387543040335`,b->-0.00005968703047293955`,c->-0.00001692245714408755`},{w,0,1}]