Initialization

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

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;

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

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_,w_]:=ϵxLDA[n]+w[[1]](ϵx1[n]-ϵx0[n])+w[[2]](ϵx2[n]-ϵx0[n])

dϵxLDAdn[n_]:=D[ϵxLDA[n],n];
dϵx0dn[n_]:=D[ϵx0[n],n];
dϵx1dn[n_]:=D[ϵx1[n],n]
dϵx2dn[n_]:=D[ϵx2[n],n]
dϵxeLDAdn[n_,w_]:=dϵxLDAdn[n]+w[[1]](dϵx1dn[n]-dϵx0dn[n])+w[[2]](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.005409938954959674`,c=0.08307660433135138`},(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.05060204953178828`,c=0.033141734861014335`},(a/(1+b (n^(-1/6))+c (n^(-1/3))))]
ϵceLDA[n_,w_]:=ϵcLDA[n]+w[[1]](ϵc1[n]-ϵc0[n])+w[[2]](ϵc2[n]-ϵc0[n])

dϵcLDAdn[n_]:=D[ϵcLDA[n],n];
dϵc0dn[n_]:=D[ϵc0[n],n];
dϵc1dn[n_]:=D[ϵc1[n],n]
dϵc2dn[n_]:=D[ϵc2[n],n]
dϵceLDAdn[n_,w_]:=dϵcLDAdn[n]+w[[1]](dϵc1dn[n]-dϵc0dn[n])+w[[2]](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_,w_]:=ϵxeLDA[n,w]+ϵceLDA[n,w]

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_,w_]:=dϵxeLDAdn[n,w]+dϵceLDAdn[n,w]

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]

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