Initialization

HaToeV=27.21138602;

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

HaTokJ=2625.5002;

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

StatIP[set_,ref_]:=SetAccuracy[{Mean[Abs[set-ref]],RootMeanSquare[set-ref],Max[Abs[set-ref]]},3];

SOPPA

Definition

In the spin-orbital basis, the SOPPA equations read

(({{(Subscript[A, 11]^((0+1+2))), (Subscript[A, 12]^((1))), (Subscript[B, 11]^((0+1+2))), 0}, {(Subscript[A, 21]^((1))), (Subscript[A, 22]^((0))), 0, 0}, {(Subscript[B, 11]^((0+1+2))), 0, (Subscript[A, 11]^((0+1+2))), (Subscript[A, 12]^((1)))}, {0, 0, (Subscript[A, 21]^((1))), (Subscript[A, 22]^((0)))}})).(({{(Subscript[X, 1])}, {(Subscript[X, 2])}, {(Subscript[Y, 1])}, {(Subscript[Y, 2])}}))=ω(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}})).(({{(Subscript[X, 1])}, {(Subscript[X, 2])}, {(Subscript[Y, 1])}, {(Subscript[Y, 2])}}))

with the following definitions for the various blocks

(Subscript[((Subscript[A, 11]^((0+1)))), ai,jb])=((Subscript[ϵ, a])-(Subscript[ϵ, i]))(Subscript[δ, ij])(Subscript[δ, ab])+(ai|jb)-(ab|ij)

(Subscript[((Subscript[B, 11]^((0+1)))), ai,jb])=(ai|bj)-(aj|bi)

(Subscript[((Subscript[A, 11]^((0+1+2)))), ia,jb])=(Subscript[((Subscript[A, 11]^((0+1)))), ia,jb])+((Subscript[δ, ab])/2)(∑_kcd)(((kd||jc)(ck||di))/((Subscript[ϵ, c])+(Subscript[ϵ, d])-(Subscript[ϵ, i])-(Subscript[ϵ, k])))-((Subscript[δ, ij])/2)(∑_klc)(((kc||lb)(ck||la))/((Subscript[ϵ, a])+(Subscript[ϵ, c])-(Subscript[ϵ, k])-(Subscript[ϵ, l])))

(Subscript[((Subscript[B, 11]^((0+1+2)))), ia,jb])=(Subscript[((Subscript[B, 11]^((0+1)))), ia,jb])-(∑_md)(((cd||ad)(di||em))/((Subscript[ϵ, d])+(Subscript[ϵ, e])-(Subscript[ϵ, i])-(Subscript[ϵ, m])))-(∑_me)(((ce||mi)(ak||me))/((Subscript[ϵ, d])+(Subscript[ϵ, e])-(Subscript[ϵ, i])-(Subscript[ϵ, m])))-(1/2)(∑_de)(((cd||ad)(dk||ei))/((Subscript[ϵ, d])+(Subscript[ϵ, e])-(Subscript[ϵ, i])-(Subscript[ϵ, k])))-(1/2)(∑_ml)(((ik||ml)(ac||ml))/((Subscript[ϵ, a])+(Subscript[ϵ, c])-(Subscript[ϵ, l])-(Subscript[ϵ, m]))) (* typos *)

(Subscript[((Subscript[A, 12]^((1)))), kc,iajb])=-(Subscript[δ, ik])(bc||aj)+(Subscript[δ, jk])(bc||ai)-(Subscript[δ, bc])(ai||kj)+(Subscript[δ, ac])(bi||kj)

(Subscript[((Subscript[A, 22]^((0)))), iajb,kcld])=((Subscript[ϵ, a])+(Subscript[ϵ, b])-(Subscript[ϵ, i])-(Subscript[ϵ, j]))(Subscript[δ, ik])(Subscript[δ, ac])(Subscript[δ, bd])(Subscript[δ, jl])

Folding the SOPPA kernel

Let’s first enforce the Tamm-Dancoff approximation by setting the (Subscript[B, 11]^((0+1+2))) block to zero

(({{(Subscript[A, 11]^((0+1+2))), (Subscript[A, 12]^((1))), 0, 0}, {(Subscript[A, 21]^((1))), (Subscript[A, 22]^((0))), 0, 0}, {0, 0, (Subscript[A, 11]^((0+1+2))), (Subscript[A, 12]^((1)))}, {0, 0, (Subscript[A, 21]^((1))), (Subscript[A, 22]^((0)))}})).(({{(Subscript[X, 1])}, {(Subscript[X, 2])}, {(Subscript[Y, 1])}, {(Subscript[Y, 2])}}))=ω(({{1, 0, 0, 0}, {0, 1, 0, 0}, {0, 0, -1, 0}, {0, 0, 0, -1}})).(({{(Subscript[X, 1])}, {(Subscript[X, 2])}, {(Subscript[Y, 1])}, {(Subscript[Y, 2])}}))

This matrix is block diagonal now, so we can consider the upper block only without loss of generality

(({{(Subscript[A, 11]^((0+1+2))), (Subscript[A, 12]^((1)))}, {(Subscript[A, 21]^((1))), (Subscript[A, 22]^((0)))}}))(({{(Subscript[X, 1])}, {(Subscript[X, 2])}}))=ω(({{(Subscript[X, 1])}, {(Subscript[X, 2])}}))

Now let’s fold the (Subscript[A, 22]^((0))) block (which has a diagonal structure) onto the (Subscript[A, 11]^((0+1+2))) block

(Subscript[(OverTilde[A]), 11]^((0+1+2)))(Subscript[X, 1])=ω (Subscript[X, 1])

with

(Subscript[(OverTilde[A]), 11]^((0+1+2)))=(Subscript[A, 11]^((0+1+2)))+(Subscript[A, 12]^((1))).(((ω 1-(Subscript[A, 22]^((0)))))^(-1)).(Subscript[A, 21]^((1)))

The elements of this effective Hamiltonian are then

(Subscript[((Subscript[(OverTilde[A]), 11]^((0+1+2)))), ia,jb])=(Subscript[((Subscript[A, 11]^((0+1+2)))), ia,jb])+(∑_kcld)(((Subscript[((Subscript[A, 12]^((1)))), ia,kcld])(Subscript[((Subscript[A, 21]^((1)))), kcld,jb]))/(ω-((Subscript[ϵ, c])+(Subscript[ϵ, d])-(Subscript[ϵ, k])-(Subscript[ϵ, l]))))

(Subscript[((Subscript[(OverTilde[A]), 11]^((0+1+2)))), ia,jb])=(Subscript[((Subscript[A, 11]^((0+1+2)))), ia,jb])+(∑_kcld)(1/(ω-((Subscript[ϵ, c])+(Subscript[ϵ, d])-(Subscript[ϵ, k])-(Subscript[ϵ, l]))))((dk||il) (Subscript[δ, c,a])-(ck||il) (Subscript[δ, d,a])-(da||cl) (Subscript[δ, k,i])+(da||ck) (Subscript[δ, l,i])) ((dk||jl) (Subscript[δ, c,b])-(ck||jl) (Subscript[δ, d,b])-(db||cl) (Subscript[δ, k,j])+(db||ck) (Subscript[δ, l,j]))

(∑_kld)(((dk||il) (dk||jl) )/(ω-(2(Subscript[ϵ, a])-(Subscript[ϵ, k])-(Subscript[ϵ, l]))))(Subscript[δ, c,a])(Subscript[δ, c,b])-(∑_kl)(((bk||il) (ak||jl) )/(ω-((Subscript[ϵ, a])+(Subscript[ϵ, b])-(Subscript[ϵ, k])-(Subscript[ϵ, l]))))(Subscript[δ, c,a])(Subscript[δ, d,b])

(∑_kl)((<bi||kl><kl||aj> )/(ω-((Subscript[ϵ, a])+(Subscript[ϵ, b])-(Subscript[ϵ, k])-(Subscript[ϵ, l]))))(Subscript[δ, c,a])(Subscript[δ, d,b])

A12[{k_,c_},{i_,a_,j_,b_}]:=-(Subscript[δ, i,k])ERI[b,c,a,j]+(Subscript[δ, j,k])ERI[b,c,a,i]-(Subscript[δ, b,c])ERI[a,i,k,j]+(Subscript[δ, a,c])ERI[b,i,k,j]

A12[{i,a},{k,c,l,d}]A12[{j,b},{k,c,l,d}]

Dynamic correction for BSE

Derivation

We have the following dynamic BSE eigenproblem

(({{A[ω], B[ω]}, {-B[ω], -A[ω]}}))(({{X[ω]}, {Y[ω]}}))=ω(({{X[ω]}, {(Y^x)[ω]}}))

where the frequency-dependent  A[ω] and B[ω] matrices are given by the usual expression (for singlets)

(Subscript[A, ia,jb])[ω]=(Subscript[δ, ij])(Subscript[δ, ab])+2(ia|jb)-(Subscript[W, ij,ab])[ω]

(Subscript[B, ia,jb])[ω]=2(ia|bj)-(Subscript[W, ib,aj])[ω]

Thanks to Xavier’s magic tricks, we know that he screened Coulomb potential reads

(Subscript[W, ij,ab])[ω]=(ij|ab)+2(∑_m)[ij|m]((1/(ω-(Subscript[Ω, m]^RPA)-((Subscript[ϵ, b])-(Subscript[ϵ, i]))+I η))+(1/(ω-(Subscript[Ω, m]^RPA)-((Subscript[ϵ, a])-(Subscript[ϵ, j]))+I η)))[ab|m]

with

[pq|m]=(∑_ia)(pq|ia)(Subscript[((X^RPA)+(Y^RPA)), m]^ia)

The static version of the screened Coulomb potential reads

(Subscript[W, ij,ab]^stat)=(ij|ab)+2(∑_m)[ij|m]((1/(-(Subscript[Ω, m]^RPA)+I η))+(1/(-(Subscript[Ω, m]^RPA)+I η)))[ab|m]

Let us decompose this eigenproblem as a zeroth-order static part and a first-order dynamic part

(({{A[ω], B[ω]}, {-B[ω], -A[ω]}}))=(({{(A^((0))), (B^((0)))}, {-(B^((0))), -(A^((0)))}}))+(({{(A^((1)))[ω], (B^((1)))[ω]}, {-(B^((1)))[ω], -(A^((1)))[ω]}}))

where

(Subscript[A, ia,jb]^((0)))=(Subscript[δ, ij])(Subscript[δ, ab])((Subscript[ϵ, a])-(Subscript[ϵ, i]))+2(ia|jb)-(Subscript[W, ij,ab]^state)

(Subscript[B, ia,jb]^((0)))=2(ia|bj)-(Subscript[W, ib,aj]^stat)

and

(Subscript[A, ia,jb]^((1)))=-(Subscript[W, ij,ab])[ω]+(Subscript[W, ij,ab]^stat)

(Subscript[B, ia,jb]^((1)))=-(Subscript[W, ib,aj])[ω]+(Subscript[W, ib,aj]^stat)

The mth BSE excitation energy and its corresponding eigenvector can then decomposed as

(Subscript[Ω, m])=(Subscript[Ω, m]^((0)))+(Subscript[Ω, m]^((1)))+…

(({{(Subscript[X, m])}, {(Subscript[Y, m])}}))=(({{(Subscript[X, m]^((0)))}, {(Subscript[Y, m]^((0)))}}))+(({{(Subscript[X, m]^((1)))}, {(Subscript[Y, m]^((1)))}}))+…

Solving the zeroth-order static problem yields

(({{(A^((0))), (B^((0)))}, {-(B^((0))), -(A^((0)))}}))(({{(Subscript[X, m]^((0)))}, {(Subscript[Y, m]^((0)))}}))=(Subscript[Ω, m]^((0)))(({{(Subscript[X, m]^((0)))}, {(Subscript[Y, m]^((0)))}}))

Thanks to first-order perturbation theory, the first-order correction to the mth excitation energy is

(Subscript[Ω, m]^((1)))=(({{(Subscript[X, m]^((0)))}, {(Subscript[Y, m]^((0)))}}))\[Transpose](({{(A^((1)))[(Subscript[Ω, m]^((0)))], (B^((1)))[(Subscript[Ω, m]^((0)))]}, {-(B^((1)))[-(Subscript[Ω, m]^((0)))], -(A^((1)))[-(Subscript[Ω, m]^((0)))]}}))(({{(Subscript[X, m]^((0)))}, {(Subscript[Y, m]^((0)))}}))
        =(({{(Subscript[X, m]^((0)))}, {(Subscript[Y, m]^((0)))}}))\[Transpose](({{(A^((1)))[(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))+(B^((1)))[(Subscript[Ω, m]^((0)))](Subscript[Y, m]^((0)))}, {-(B^((1)))[-(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))-(A^((1)))[-(Subscript[Ω, m]^((0)))](Subscript[Y, m]^((0)))}}))
        =(Subscript[X, m]^((0)))\[Transpose]((A^((1)))[(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))+(B^((1)))[(Subscript[Ω, m]^((0)))](Subscript[Y, m]^((0))))+(Subscript[Y, m]^((0)))\[Transpose](-(B^((1)))[-(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))-(A^((1)))[-(Subscript[Ω, m]^((0)))](Subscript[Y, m]^((0))))
        =(Subscript[X, m]^((0)))\[Transpose](A^((1)))[(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))-(Subscript[Y, m]^((0)))\[Transpose](A^((1)))[-(Subscript[Ω, m]^((0)))](Subscript[Y, m]^((0)))+(Subscript[X, m]^((0)))\[Transpose](B^((1)))[(Subscript[Ω, m]^((0)))](Subscript[Y, m]^((0)))-(Subscript[Y, m]^((0)))\[Transpose](B^((1)))[-(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))

If one enforces the TDA, we recover Rohlfing’s expression (I think), i.e.,

(Subscript[Ω, m]^((1)))=(Subscript[X, m]^((0)))\[Transpose](A^((1)))[(Subscript[Ω, m]^((0)))](Subscript[X, m]^((0)))

This correction can be renormalized by computing

(((Subscript[(Subscript[Z, m])=(1-(Subscript[X, m]^((0)))\[Transpose]((∂(A^((1)))[ω])/(∂ω))\[RightBracketingBar], ω=(Subscript[Ω, m]^((0)))])(Subscript[X, m]^((0)))))^(-1))

with

((∂(A^((1)))[ω])/(∂ω))=-((∂W[ω])/(∂ω)) \[Implies] ((∂(Subscript[W, ij,ab])[ω])/(∂ω))=-2(∑_m)[ij|m]((1/(((ω-(Subscript[Ω, m]^RPA)-((Subscript[ϵ, b])-(Subscript[ϵ, i]))+I η))^2))+(1/(((ω-(Subscript[Ω, m]^RPA)-((Subscript[ϵ, a])-(Subscript[ϵ, j]))+I η))^2)))[ab|m]

Results

Singlets

Ω={
{"^1\\Pi\\text{(CT)}"},
{"^1B_1\\text{(R)}","^1A_2\\text{(R)}","^1A_1\\text{(R)}"},
{"^1\\Pi_g(n \\to \\pi^*)","^1\\Sigma_u^-(\\pi \\to \\pi^*)","^1\\Delta_u(\\pi \\to \\pi^*)","^1\\Sigma_g^+\\text{(R)}","^1\\Pi_u\\text{(R)}","^1\\Sigma_u^+\\text{(R)}","^1\\Pi_u\\text{(R)}"},
{"^1\\Pi(n \\to \\pi^*)","^1\\Sigma^-(\\pi \\to \\pi^*)","^1\\Delta(\\pi \\to \\pi^*)","^1\\Sigma^+\\text{(R)}","^1\\Sigma^+\\text{(R)}","^1\\Pi\\text{(R)}"},
{"^1\\Sigma_{u}^-(\\pi \\to \\pi^*)","^1\\Delta_{u}(\\pi \\to \\pi^*)"},
{"^1B_{3u}\\text{(R)}","^1B_{1u}(\\pi \\to \\pi^*)","^1B_{1g}\\text{(R)}"},
{"^1A_2(n \\to \\pi^*)","^1B_2\\text{(R)}","^1B_2\\text{(R)}","^1A_1\\text{(R)}","^1A_2\\text{(R)}","^1B_1(\\sigma \\to \\pi^*)","^1A_1(\\pi \\to \\pi^*)"}
};

{"^1A''(n \\to \\pi^*)","^1A'\\text{(R)}"},

ΩTBE={
{7.84},
{7.17,8.92,9.52},
{9.34,9.88,10.29,12.98,13.03,13.09,13.46},
{8.49,9.92,10.06,10.95,11.52,11.72},
(*{1.74,6.27},*)
{7.10,7.44},
{7.39,7.93,8.08},
{3.98,7.23,8.13,8.23,8.67,9.22,9.43}
};
ΩdynBSE={
{8.19},
{8.00,9.72,10.35},
{9.99,9.66,10.33,13.57,13.94,13.91,14.21},
{9.19,9.90,10.39,11.85,12.37,12.32},
(*{1.98,7.01},*)
{7.05,7.46},
{7.62,8.03,8.26},
{4.68,7.85,8.72,8.84,8.85,9.77,9.81}
};
ΩstatBSE={
{8.30},
{8.09,9.79,10.42},
{10.42,10.11,10.75,13.60,13.98,13.98,14.24},
{9.54,10.25,10.71,11.88,12.39,12.37},
(*{2.46,7.05},*)
{7.37,7.74},
{7.64,8.18,8.29},
{5.03,7.87,8.76,8.85,8.87,10.18,10.05}
};
ΩCC2={
{7.96},
{7.23,8.89,9.58},
{9.44,10.32,10.86,12.83,13.15,12.89,13.96},
{8.64,10.30,10.60,11.11,11.63,11.83},
(*{1.74,5.72},*)
{7.26,7.59},
{7.29,7.92,7.95},
{4.07,6.56,7.57,7.52,8.04,9.32,9.54}
};
ΩCCSD={
{7.91},
{7.60,9.36,9.96},
{9.41,10.00,10.44,13.15,13.43,13.26,13.67},
{8.59,9.99,10.12,11.22,11.75,11.96},
(*{1.76,6.30},*)
{7.15,7.48},
{7.42,8.02,8.08},
{4.01,7.23,8.12,8.21,8.65,9.28,9.67}};
ΩADC2={
{7.97},
{7.18,8.84,9.52},
{9.48,10.26,10.79,12.99,13.32,13.07,14.00},
{8.69,10.03,10.30,11.32,11.83,12.03},
(*{1.68,5.73},*)
{7.24,7.56},
{7.34,7.91,7.99},
{3.92,6.50,7.53,7.47,7.99,9.17,9.46}
};
ΩCISD={
{6.07},
{7.62,9.41,9.99},
{9.66,10.31,10.85,13.67,13.64,13.75,14.52},
{8.78,10.13,10.41,11.48,11.71,12.06},
(*{1.80,5.81},*)
{7.28,7.62},
{7.35,7.95,8.01},
{4.04,6.64,7.56,8.16,8.04,9.38,9.08}
};

SetAccuracy[{Mean[Abs[Flatten[ΩstatBSE-ΩTBE]]],Mean[Abs[Flatten[ΩdynBSE-ΩTBE]]],Mean[Abs[Flatten[ΩCISD-ΩTBE]]],Mean[Abs[Flatten[ΩADC2-ΩTBE]]],Mean[Abs[Flatten[ΩCCSD-ΩTBE]]],Mean[Abs[Flatten[ΩCC2-ΩTBE]]]},3]

SetAccuracy[{Mean[Flatten[ΩstatBSE-ΩTBE]],Mean[Flatten[ΩdynBSE-ΩTBE]],Mean[Flatten[ΩCISD-ΩTBE]],Mean[Flatten[ΩADC2-ΩTBE]],Mean[Flatten[ΩCCSD-ΩTBE]],Mean[Flatten[ΩCC2-ΩTBE]]},3]

SetAccuracy[{RootMeanSquare[Flatten[ΩstatBSE-ΩTBE]],RootMeanSquare[Flatten[ΩdynBSE-ΩTBE]],RootMeanSquare[Flatten[ΩCISD-ΩTBE]],RootMeanSquare[Flatten[ΩADC2-ΩTBE]],RootMeanSquare[Flatten[ΩCCSD-ΩTBE]],RootMeanSquare[Flatten[ΩCC2-ΩTBE]]},3]

SetAccuracy[{Max[Flatten[ΩstatBSE-ΩTBE]],Max[Flatten[ΩdynBSE-ΩTBE]],Max[Flatten[ΩCISD-ΩTBE]],Max[Flatten[ΩADC2-ΩTBE]],Max[Flatten[ΩCCSD-ΩTBE]],Max[Flatten[ΩCC2-ΩTBE]]},3]

SetAccuracy[{Min[Flatten[ΩstatBSE-ΩTBE]],Min[Flatten[ΩdynBSE-ΩTBE]],Min[Flatten[ΩCISD-ΩTBE]],Min[Flatten[ΩADC2-ΩTBE]],Min[Flatten[ΩCCSD-ΩTBE]],Min[Flatten[ΩCC2-ΩTBE]]},3]

SizeFrameLabel=18;
SizeChartLabel=18;
Options$Chart={
PlotTheme->"Classic",
BarSpacing->{0.1,3/4},
Joined->False,
Frame->True,
ImageSize->1100,
GridLines->{None,Automatic},
BaseStyle->18,
ImagePadding->{{All,All},{All,All}},
BarSpacing->{0,0},
FrameTicks->{Automatic,Automatic},
(*PlotLabel->MaTeX["\\text{Ionization potentials from $G_0W_0$@PBE}",FontSize->SizeFrameLabel],*)
FrameLabel->{None,MaTeX["\\text{Error wrt TBE (eV)}",FontSize->SizeFrameLabel]},
PlotRange->{Automatic,Automatic},
(*Epilog->{Inset[ChemicalData["Adenine","MoleculePlot"],Scaled[{0.17,0.4}],Automatic,4],
Inset[ChemicalData["Cytosine","MoleculePlot"],Scaled[{0.36,0.4}],Automatic,4],
Inset[ChemicalData["Guanine","MoleculePlot"],Scaled[{0.55,0.4}],Automatic,4],
Inset[ChemicalData["Thymine","MoleculePlot"],Scaled[{0.75,0.4}],Automatic,4]
,Inset[ChemicalData["Uracil","MoleculePlot"],Scaled[{0.95,0.4}],Automatic,4]},*)
AspectRatio->1/6
};
Show[{
BarChart[Table[(ΩstatBSE-ΩTBE)[[k]],{k,Length[ΩTBE]}],(*ChartLegends->Placed[{MaTeX["\\text{Static BSE}",FontSize->SizeFrameLabel]},Top],*)PlotRange->{{1,33},{-0.2,1.3}},ChartElementFunction->"Rectangle",ChartStyle->Directive[Red,Opacity[0.5]],ChartLabels->{Placed[{"HCl",Style["(Subscript[H, 2])O",SingleLetterItalics->False],Style["(Subscript[N, 2])",SingleLetterItalics->False],"CO",(*"HNO",*)Style["(Subscript[C, 2])(Subscript[H, 2])",SingleLetterItalics->False],Style["(Subscript[C, 2])(Subscript[H, 4])",SingleLetterItalics->False],Style["(Subscript[CH, 2])O",SingleLetterItalics->False]},Above],
Placed[Table[
Rotate[MaTeX[Ω[[k,l]],FontSize->SizeFrameLabel],90 Degree],{k,Length[ΩTBE]},{l,Length[ΩTBE[[k]]]}],Axis]},Options$Chart]
,
BarChart[Table[(ΩdynBSE-ΩTBE)[[k]],{k,Length[ΩTBE]}],(*ChartLegends->Placed[{MaTeX["\\text{Dynamic BSE}",FontSize->SizeFrameLabel]},Top],*)PlotRange->{-0.2,1.1},ChartElementFunction->"Rectangle",ChartStyle->Directive[Blue,Opacity[0.5]],Options$Chart]
}]
Export["fig1a.pdf",%];

Triplets

{"^3A''(n \\to \\pi^*)","^3A'(\\pi \\to \\pi^*)"},

Ω={
{"^3B_1\\text{(R)}","^3A_2\\text{(R)}","^3A_1\\text{(R)}"},
{"^3\\Sigma_u^+(\\pi \\to \\pi^*)","^3\\Pi_g(n \\to \\pi^*)","^3\\Delta_u(\\pi \\to \\pi^*)","^3\\Sigma_u^-(\\pi \\to \\pi^*)"},
{"^3\\Pi(n \\to \\pi^*)","^3\\Sigma^+(\\pi \\to \\pi^*)","^3\\Delta(\\pi \\to \\pi^*)","^3\\Sigma_u^-(\\pi \\to \\pi^*)","^3\\Sigma_u^+\\text{(R)}"},
{"^3\\Sigma_{u}^+(\\pi \\to \\pi^*)","^3\\Delta_{u}(\\pi \\to \\pi^*)","^3\\Sigma_{u}^-(\\pi \\to \\pi^*)"},
{"^3B_{1u}(\\pi \\to \\pi^*)","^3B_{3u}\\text{(R)}","^3B_{1g}\\text{(R)}"},
{"^3A_2(n \\to \\pi^*)","^3A_1(\\pi \\to \\pi^*)","^3B_2\\text{(R)}"}
};
ΩTBE={
{6.92,8.91,9.30},
{7.70,8.01,8.87,9.66},
{6.28,8.45,9.27,9.80,10.47},
(*{0.88,5.61},*)
{5.53,6.40,7.08},
{4.54,7.23,7.98},
{3.58,6.06,7.06}
};
ΩdynBSE={
{7.48,9.50,9.66},
{7.38,8.10,8.48,9.66},
{6.25,8.06,8.96,9.90,11.07},
(*{0.67,4.87},*)
{5.32,6.23,7.05},
{4.49,7.42,8.19},
{3.87,5.75,7.56}
};
ΩstatBSE={
{7.62,9.61,9.80},
{8.02,8.66,9.04,10.11},
{6.80,8.56,9.39,10.25,11.17},
(*{1.27,5.55},*)
{5.83,6.64,7.37},
{4.95,7.46,8.23},
{4.28,6.31,7.60}
};
ΩCC2={
{6.91,8.77,9.20},
{8.19,8.19,9.30,10.29},
{6.42,8.72,9.56,10.27,10.60},
(*{0.84,5.44},*)
{5.76,6.60,7.29},
{4.59,7.19,7.91},
{3.59,6.30,6.44}
};
ΩCCSD={
{7.20,9.20,9.49},
{7.66,8.09,8.91,9.83},
{6.36,8.34,9.23,9.81,10.71},
(*{0.85,5.49},*)
{5.45,6.41,7.12},
{4.46,7.29,8.03},
{3.56,5.97,7.08}
};
ΩADC2={
{6.86,8.72,9.15},
{8.15,8.20,9.25,10.23},
{6.45,8.54,9.33,10.01,10.83},
(*{0.78,5.46},*)
{5.75,6.57,7.27},
{4.59,7.23,7.95},
{3.46,6.20,6.39}
};
ΩCISD={
{7.25,9.24,9.54},
{8.20,8.33,9.30,10.29},
{6.51,8.63,9.44,10.10,10.98},
(*{0.91,5.72},*)
{5.79,6.62,7.31},
{4.62,7.26,7.97},
{3.58,6.27,6.66}
};

SetAccuracy[{Mean[Abs[Flatten[ΩstatBSE-ΩTBE]]],Mean[Abs[Flatten[ΩdynBSE-ΩTBE]]],Mean[Abs[Flatten[ΩCISD-ΩTBE]]],Mean[Abs[Flatten[ΩADC2-ΩTBE]]],Mean[Abs[Flatten[ΩCCSD-ΩTBE]]],Mean[Abs[Flatten[ΩCC2-ΩTBE]]]},3]

SetAccuracy[{Mean[Flatten[ΩstatBSE-ΩTBE]],Mean[Flatten[ΩdynBSE-ΩTBE]],Mean[Flatten[ΩCISD-ΩTBE]],Mean[Flatten[ΩADC2-ΩTBE]],Mean[Flatten[ΩCCSD-ΩTBE]],Mean[Flatten[ΩCC2-ΩTBE]]},3]

SetAccuracy[{RootMeanSquare[Flatten[ΩstatBSE-ΩTBE]],RootMeanSquare[Flatten[ΩdynBSE-ΩTBE]],RootMeanSquare[Flatten[ΩCISD-ΩTBE]],RootMeanSquare[Flatten[ΩADC2-ΩTBE]],RootMeanSquare[Flatten[ΩCCSD-ΩTBE]],RootMeanSquare[Flatten[ΩCC2-ΩTBE]]},3]

SetAccuracy[{Max[Flatten[ΩstatBSE-ΩTBE]],Max[Flatten[ΩdynBSE-ΩTBE]],Max[Flatten[ΩCISD-ΩTBE]],Max[Flatten[ΩADC2-ΩTBE]],Max[Flatten[ΩCCSD-ΩTBE]],Max[Flatten[ΩCC2-ΩTBE]]},3]

SetAccuracy[{Min[Flatten[ΩstatBSE-ΩTBE]],Min[Flatten[ΩdynBSE-ΩTBE]],Min[Flatten[ΩCISD-ΩTBE]],Min[Flatten[ΩADC2-ΩTBE]],Min[Flatten[ΩCCSD-ΩTBE]],Min[Flatten[ΩCC2-ΩTBE]]},3]

SizeFrameLabel=18;
SizeChartLabel=18;
Options$Chart={
PlotTheme->"Classic",
BarSpacing->{0.1,3/4},
Joined->False,
Frame->True,
ImageSize->1100,
GridLines->{None,Automatic},
BaseStyle->18,
ImagePadding->{{All,All},{All,All}},
BarSpacing->{0,0},
FrameTicks->{Automatic,Automatic},
(*PlotLabel->MaTeX["\\text{Ionization potentials from $G_0W_0$@PBE}",FontSize->SizeFrameLabel],*)
FrameLabel->{None,MaTeX["\\text{Error wrt TBE (eV)}",FontSize->SizeFrameLabel]},
PlotRange->{Automatic,Automatic},
AspectRatio->1/6
};
Show[{
BarChart[Table[(ΩstatBSE-ΩTBE)[[k]],{k,Length[ΩTBE]}],ChartLegends->Placed[{MaTeX["\\text{Static BSE}",FontSize->SizeFrameLabel]},Top],PlotRange->{{1,24.5},{-0.4,0.9}},ChartElementFunction->"Rectangle",ChartStyle->Directive[Red,Opacity[0.5]],ChartLabels->{Placed[{Style["(Subscript[H, 2])O",SingleLetterItalics->False],Style["(Subscript[N, 2])",SingleLetterItalics->False],"CO",(*"HNO",*)Style["(Subscript[C, 2])(Subscript[H, 2])",SingleLetterItalics->False],Style["(Subscript[C, 2])(Subscript[H, 4])",SingleLetterItalics->False],Style["(Subscript[CH, 2])O",SingleLetterItalics->False]},Above],
Placed[Table[
Rotate[MaTeX[Ω[[k,l]],FontSize->SizeFrameLabel],90 Degree],{k,Length[ΩTBE]},{l,Length[ΩTBE[[k]]]}],Axis]},Options$Chart]
,
BarChart[Table[(ΩdynBSE-ΩTBE)[[k]],{k,Length[ΩTBE]}],ChartLegends->Placed[{MaTeX["\\text{Dynamic BSE}",FontSize->SizeFrameLabel]},Top],PlotRange->{-0.2,1.1},ChartElementFunction->"Rectangle",ChartStyle->Directive[Blue,Opacity[0.5]],Options$Chart]
}]
Export["fig1b.pdf",%];

Larger molecules

Ω={
{"^1A''(n \\to \\pi^*)","^1A'(n \\to \\pi^*)","^1A'\\text{(R)}","^3A''(n \\to \\pi^*)","^3A'(\\pi \\to \\pi^*)"},
{"^1B_u(\\pi \\to \\pi^*)","^1A_g(\\pi \\to \\pi^*)","^3B_u(\\pi \\to \\pi^*)","^3A_g(\\pi \\to \\pi^*)","^3B_g\\text{(R)}"},
{"^1\\Sigma_u^-(\\pi \\to \\pi^*)","^1\\Delta_u(\\pi \\to \\pi^*)","^3\\Sigma_u^+(\\pi \\to \\pi^*)","^3\\Delta_u(\\pi \\to \\pi^*)"},
{"^1A_u(n \\to \\pi^*)","^1B_g(n \\to \\pi^*)","^1B_u\\text{(R)}","^3A_u(n \\to \\pi^*)","^3B_g(n \\to \\pi^*)","^3B_u(\\pi \\to \\pi^*)"},
{"^1B_2(\\pi \\to \\pi^*)"}
};

ΩCC3={
{3.77,6.67,6.99,3.47,3.95},
{6.25,6.68,3.36,5.21,6.20},
{5.44,5.69,4.06,4.86},
{2.90,4.30,7.55,2.49,3.91,5.20},
{7.14}
};
ΩdynBSE={
{4.28,6.70,7.53,3.54,3.61},
{6.13,6.86,3.25,5.01,6.25},
{5.35,5.63,3.82,4.68},
{3.14,4.55,7.86,2.38,3.75,4.47},
{7.51}
};
ΩstatBSE={
{4.62,6.86,7.57,3.97,4.03},
{6.25,6.88,3.68,5.51,6.29},
{5.62,5.87,4.30,5.04},
{3.46,4.96,7.90,2.77,4.23,5.01},
{7.66}
};

SetAccuracy[{Mean[Abs[Flatten[ΩstatBSE-ΩCC3]]],Mean[Abs[Flatten[ΩdynBSE-ΩCC3]]]},3]

SetAccuracy[{Mean[Flatten[ΩstatBSE-ΩCC3]],Mean[Flatten[ΩdynBSE-ΩCC3]]},3]

SetAccuracy[{RootMeanSquare[Flatten[ΩstatBSE-ΩCC3]],RootMeanSquare[Flatten[ΩdynBSE-ΩCC3]]},3]

SetAccuracy[{Max[Flatten[ΩstatBSE-ΩCC3]],Max[Flatten[ΩdynBSE-ΩCC3]]},3]

SetAccuracy[{Min[Flatten[ΩstatBSE-ΩCC3]],Min[Flatten[ΩdynBSE-ΩCC3]]},3]

SizeFrameLabel=18;
SizeChartLabel=18;
Options$Chart={
PlotTheme->"Classic",
BarSpacing->{1/5,3/4},
Joined->False,
Frame->True,
ImageSize->1100,
GridLines->{None,Automatic},
BaseStyle->18,
ImagePadding->{{All,All},{All,All}},
BarSpacing->{0,0},
FrameTicks->{Automatic,Automatic},
(*PlotLabel->MaTeX["\\text{Ionization potentials from $G_0W_0$@PBE}",FontSize->SizeFrameLabel],*)
FrameLabel->{None,MaTeX["\\text{Error wrt CC3 (eV)}",FontSize->SizeFrameLabel]},
PlotRange->{Automatic,Automatic},
AspectRatio->1/6
};
Show[{
BarChart[Table[(ΩstatBSE-ΩCC3)[[k]],{k,Length[ΩCC3]}],ChartLegends->Placed[{MaTeX["\\text{Static BSE}",FontSize->SizeFrameLabel]},Top],PlotRange->{{1,24},{-0.2,1.0}},ChartElementFunction->"Rectangle",ChartStyle->Directive[Red,Opacity[0.5]],ChartLabels->{Placed[{"acrolein","butadiene","diacetylene","glyoxal","streptocyanine"},Above],
Placed[Table[
Rotate[MaTeX[Ω[[k,l]],FontSize->SizeFrameLabel],90 Degree],{k,Length[ΩCC3]},{l,Length[ΩCC3[[k]]]}],Axis]},Options$Chart]
,
BarChart[Table[(ΩdynBSE-ΩCC3)[[k]],{k,Length[ΩCC3]}],ChartLegends->Placed[{MaTeX["\\text{Dynamic BSE}",FontSize->SizeFrameLabel]},Top],PlotRange->{-0.2,1.1},ChartElementFunction->"Rectangle",ChartStyle->Directive[Blue,Opacity[0.5]],Options$Chart]
}]
Export["fig2.pdf",%];

Examples

Two-level system

Adiabatic connection

vvvv=1.0269071688766602;
cccc=0.76636289621137521;
vccv=0.22767049526677599;
vvcc=0.85813333443746587;
cvvv=0.31649047035074229;
vccc=0.25555354578313411;
ϵvHF=-0.91412663;	ϵcHF=1.39985934;	ΔϵHF=ϵcHF-ϵvHF;
(*ϵvGW=-0.863700318396;	ϵcGW=1.37364001917;	ΔϵGW=ϵcGW-ϵvGW;
ϵvGF2=-0.8827034015226544`;	ϵcGF2=1.3830817721794242`;	ΔϵGF2=ϵcGF2-ϵvGF2;
*)

Ω[λ_]:=ΔϵHF+2 λ vccv;

ϵvGW[λ_]:=ϵvHF+ZvGW[λ] ΣvGW[λ];
ϵcGW[λ_]:=ϵcHF+ZcGW[λ] ΣcGW[λ];
ΔϵGW[λ_]:=ϵcGW[λ]-ϵvGW[λ]

ZvGW[λ_]:=(((1+(λ^2)((2 (cvvv^2))/(((ϵvHF-ϵvHF+Ω[λ]))^2))+(λ^2)((2 (vccv^2))/(((ϵvHF-ϵcHF-Ω[λ]))^2))))^(-1))
ZcGW[λ_]:=(((1+(λ^2)((2 (vccv^2))/(((ϵcHF-ϵvHF+Ω[λ]))^2))+(λ^2)((2 (vccc^2))/(((ϵcHF-ϵcHF-Ω[λ]))^2))))^(-1))

ΣvGW[λ_]:=(λ^2)((2 (cvvv^2))/(ϵvHF-ϵvHF+Ω[λ]))+(λ^2)((2 (vccv^2))/(ϵvHF-ϵcHF-Ω[λ]))
ΣcGW[λ_]:=(λ^2)((2 (vccv^2))/(ϵcHF-ϵvHF+Ω[λ]))+(λ^2)((2 (vccc^2))/(ϵcHF-ϵcHF-Ω[λ]))

WA[ω_,λ_]:= λ vvcc+((2(λ^2)cvvv vccc)/(ω-Ω[λ]-ΔϵGW[λ]))+((2 (λ^2)cvvv vccc)/(ω-Ω[λ]-ΔϵGW[λ]));
WB[ω_,λ_]:= λ vccv+((2(λ^2)vccv vccv)/(ω-Ω[λ]))+((2(λ^2)vccv vccv)/(ω-Ω[λ]));

H1BSE[λ_]:=(({{ΔϵGW[λ]+2λ vccv-WA[ΔϵGW,λ], 2λ vccv-WB[0,λ]}, {-(2λ vccv-WB[0,λ]), -(ΔϵGW[λ]+2λ vccv-WA[ΔϵGW,λ])}}));
H3BSE[λ_]:=(({{ΔϵGW[λ]-WA[ΔϵGW,λ], -WB[0,λ]}, {-(-WB[0,λ]), -(ΔϵGW[λ]-WA[ΔϵGW,λ])}}));

H1dBSE[ω_,λ_ ]:=(({{ΔϵGW[λ]+2λ vccv-WA[ω,λ ], 2 λ vccv-WB[ω,λ ]}, {-(2 λ vccv- WB[-ω,λ ]), -(ΔϵGW[λ]+2 λ vccv-WA[-ω,λ ])}}))
H3dBSE[ω_,λ_]:=(({{ΔϵGW[λ]-WA[ω,λ ], -WB[ω,λ ]}, {-(-WB[-ω,λ ]), -(ΔϵGW[λ]-WA[-ω,λ ])}}))

H1dRPAx[ω_,λ_]:=(({{ΔϵGF2+2λ vccv-λ vvcc, 2λ vccv-λ vccv}, {-(2λ vccv-λ vccv), -(ΔϵGF2[λ]+2λ vccv-λ vvcc)}}))+(λ^2)(({{-((4cvvv vccc-(vccc^2)-(cvvv^2))/(ω-2ΔϵGF2)), -((4(vccv^2)- cccc vccv- vvvv vccv)/(2ΔϵGF2))}, {+((4(vccv^2)- cccc vccv- vvvv vccv)/(2ΔϵGF2)), +((4cvvv vccc-(vccc^2)-(cvvv^2))/(-ω-2ΔϵGF2))}}));
H3dRPAx[ω_,λ_]:=(({{ΔϵGF2[λ]-λ vvcc, -λ vccv}, {-(-λ vccv), -(ΔϵGF2[λ]-λ vvcc)}}))+(λ^2)(({{-(((vccc^2)+(cvvv^2))/(ω-2ϵcGF2+2ϵvGF2)), -(( cccc vccv+vvvv vccv)/(2ϵcGF2-2ϵvGF2))}, {+(( cccc vccv+vvvv vccv)/(2ϵcGF2-2ϵvGF2)), +(((vccc^2)+(cvvv^2))/(-ω-2ϵcGF2+2ϵvGF2))}}));

dat=Table[
Flatten[{λ,ω/.NSolve[{ω>0,Det[H1dBSE[ω,λ]-ω IdentityMatrix[2]]==0},ω]//Sort}]
,{λ,0.01,1,0.01}];

ListPlot[{dat[[;;,{1,2}]],dat[[;;,{1,3}]],dat[[;;,{1,4}]]},Joined->True]

dat=Table[
Flatten[{λ,ω/.NSolve[{ω>0,Det[H3dBSE[ω,λ]-ω IdentityMatrix[2]]==0},ω]//Sort}]
,{λ,0.01,1,0.01}];

ListPlot[{dat[[;;,{1,2}]],dat[[;;,{1,3}]],dat[[;;,{1,4}]]},Joined->True]

dat=Table[
Flatten[{λ,ω/.NSolve[{ω>0,Det[H1dRPAx[ω,λ]-ω IdentityMatrix[2]]==0},ω]//Sort}]
,{λ,0.01,1,0.01}];

ListPlot[{dat[[;;,{1,2}]],dat[[;;,{1,3}]]},Joined->True]

dat=Table[
Flatten[{λ,ω/.NSolve[{ω>0,Det[H3dRPAx[ω,λ]-ω IdentityMatrix[2]]==0},ω]//Sort}]
,{λ,0.01,1,0.01}];

ListPlot[{dat[[;;,{1,2}]],dat[[;;,{1,3}]]},Joined->True]

Fully dynamical approach

vvvv=1.0269071688766602;
cccc=0.76636289621137521;
vccv=0.22767049526677599;
vvcc=0.85813333443746587;
cvvv=0.31649047035074229;
vccc=0.25555354578313411;
ϵvHF=-0.91412663;	ϵcHF=1.39985934;	ΔϵHF=ϵcHF-ϵvHF;
(*ϵvGW=-0.863700318396;	ϵcGW=1.37364001917;	ΔϵGW=ϵcGW-ϵvGW;
ϵvGF2=-0.8827034015226544`;	ϵcGF2=1.3830817721794242`;	ΔϵGF2=ϵcGF2-ϵvGF2;
*)

Ω[λ_]:=ΔϵHF+2 λ vccv;
ϵvGW[ω_,λ_]:=ϵvHF+ΣvGW[ω,λ];
ϵcGW[ω_,λ_]:=ϵcHF+ΣcGW[ω,λ];
ΔϵGW[ω_,λ_]:=ϵcGW[ω,λ]-ϵvGW[ω,λ]

ΣvGW[ω_,λ_]:=(λ^2)((2 (cvvv^2))/(ω-ϵvHF+Ω[λ]))+(λ^2)((2 (vccv^2))/(ω-ϵcHF-Ω[λ]))
ΣcGW[ω_,λ_]:=(λ^2)((2 (vccv^2))/(ω-ϵvHF+Ω[λ]))+(λ^2)((2 (vccc^2))/(ω-ϵcHF-Ω[λ]))

WA[ω_,λ_]:= λ vvcc+((2(λ^2)cvvv vccc)/(ω-Ω[λ]-ΔϵGW[ω,λ]))+((2 (λ^2)cvvv vccc)/(ω-Ω[λ]-ΔϵGW[ω,λ]));
WB[ω_,λ_]:= λ vccv+((2(λ^2)vccv vccv)/(ω-Ω[λ]))+((2(λ^2)vccv vccv)/(ω-Ω[λ]));

H1dBSE[ω_,λ_ ]:=(({{ΔϵGW[ω,λ]+2λ vccv-WA[ω,λ ], 2 λ vccv-WB[ω,λ ]}, {-(2 λ vccv- WB[-ω,λ ]), -(ΔϵGW[-ω,λ]+2 λ vccv-WA[-ω,λ ])}}))
H3dBSE[ω_,λ_]:=(({{ΔϵGW[ω,λ]-WA[ω,λ ], -WB[ω,λ ]}, {-(-WB[-ω,λ ]), -(ΔϵGW[-ω,λ]-WA[-ω,λ ])}}))

With[{λ=1},NSolve[{ω>0,Det[ΔϵGW[ω,λ]+2λ vccv-WA[ω,λ ]-ω IdentityMatrix[1]]==0},ω]]

With[{λ=1},NSolve[{ω>0,Det[H1dBSE[ω,λ]-ω IdentityMatrix[2]]==0},ω]]

With[{λ=1},NSolve[{ω>0,Det[H3dBSE[ω,λ]-ω IdentityMatrix[2]]==0},ω]]

dat=Table[
Flatten[{λ,ω/.NSolve[{ω>0,Det[ΔϵGW[ω,λ]+2λ vccv-WA[ω,λ ]-ω IdentityMatrix[1]]==0},ω]//Sort}]
,{λ,0.01,1,0.01}];

ListPlot[{dat[[;;,{1,2}]],dat[[;;,{1,3}]],dat[[;;,{1,4}]],dat[[;;,{1,5}]]},Joined->True]

Parameters

(* parameters for (Subscript[H, 2]) *)
vvvv=0.67459407408316485;
cccc=0.69749533680191667;
vvcc=0.66356398189533250;
vccv=0.18125791779638803;
cvvv=0;
vccc=0;
ϵvHF=-0.57820296;	ϵcHF=0.67026771;	ΔϵHF=ϵcHF-ϵvHF;
ϵvGW=-0.6009992283369916`;	ϵcGW=0.6930639985092535`;	ΔϵGW=ϵcGW-ϵvGW;
Ω=ΔϵHF+2 vccv;
ϵvGF2=-0.5912918580543514`;	ϵcGF2=0.6833566282266132`;	ΔϵGF2=ϵcGF2-ϵvGF2;

(*(* parameters for stretched (Subscript[H, 2]) *)
vvvv=0.45759218516850819;
cccc=0.46005248040625285;
vvcc=0.45881335974065185;
vccv=0.31579709471060158;
cvvv=0;
vccc=0;
ϵvHF=-0.15415523;	ϵcHF=0.00561232;	ΔϵHF=ϵcHF-ϵvHF;
ϵvGW=-0.3285878930763851`;	ϵcGW=0.1688203605881594`;	ΔϵGW=ϵcGW-ϵvGW;
Ω=ΔϵHF+2 vccv;
ϵvGF2=-0.31175971682459713`;	ϵcGF2=0.15199214758704896`;	ΔϵGF2=ϵcGF2-ϵvGF2;*)

(* parameters for HeH+ *)
vvvv=0.94309859025505716;
cccc=0.75252617021781365;
vvcc=0.66025406610737403;
vccv=0.14539723245724573;
cvvv=-0.17296841369167024;
vccc=0.037282147343908126;
ϵvHF=-1.63280248;	ϵcHF=-0.17248357;	ΔϵHF=ϵcHF-ϵvHF;
ϵvGW=-1.6122824823312696`;	ϵcGW=-0.16096309819649535`;	ΔϵGW=ϵcGW-ϵvGW;
Ω=ΔϵHF+2 vccv;
ϵvGF2=-1.6197686500645219`;	ϵcGF2=-0.16621674458903582`;	ΔϵGF2=ϵcGF2-ϵvGF2;

(* parameters for He *)
vvvv=1.0269071688766602;
cccc=0.76636289621137521;
vvcc=0.85813333443746587;
vccv=0.22767049526677599;
cvvv=0.31649047035074229;
vccc=0.25555354578313411;
ϵvHF=-0.91412663;	ϵcHF=1.39985934;	ΔϵHF=ϵcHF-ϵvHF;
ϵvGW=-0.863700318396;	ϵcGW=1.37364001917;	ΔϵGW=ϵcGW-ϵvGW;
Ω=ΔϵHF+2 vccv;
ϵvGF2=-0.8827034015226544`;	ϵcGF2=1.3830817721794242`;	ΔϵGF2=ϵcGF2-ϵvGF2;

Exact Hamiltonian

SortEigensystem[2ϵvHF-vvvv+Eigensystem[(({{0, 0, 0, vccv}, {0, ΔϵHF-vvcc+vccv, vccv, vccc-cvvv}, {0, vccv, ΔϵHF-vvcc+vccv, vccc-cvvv}, {vccv, vccc-cvvv, vccc-cvvv, 2ΔϵHF+vvvv+cccc-4vvcc+2vccv}}))]]

Eigensystem[(({{0, 0, 0, vccv}, {0, ΔϵHF+vccv-vvcc, vccv, vccc-cvvv}, {0, vccv, ΔϵHF+vccv-vvcc, vccc-cvvv}, {vccv, vccc-cvvv, vccc-cvvv, 2ΔϵHF+vvvv+cccc-4vvcc+2vccv}}))]

{Eigenvalues[(({{0, 0, vccv}, {0, ΔϵHF+2vccv-vvcc, (Sqrt[2])(vccc-cvvv)}, {vccv, (Sqrt[2])(vccc-cvvv), 2ΔϵHF+vvvv+cccc-4vvcc+2vccv}}))],ΔϵHF-vvcc}

Solve[Det[(({{0, 0}, {0, ΔϵHF+2vccv-vvcc}}))+(1/(ω-(2ΔϵHF+vvvv+cccc-4vvcc+2vccv)))(({{(vccv^2), (Sqrt[2]) (vccc-cvvv) vccv}, {(Sqrt[2]) (vccc-cvvv) vccv, 2 (((vccc-cvvv))^2)}}))-ω IdentityMatrix[2]]==0,ω]

Ec=-0.009629680837542426`;

Solve[Det[(({{-Ec, 0}, {0, -Ec+ΔϵHF+2vccv-vvcc}}))+(1/(ω-(-Ec+2ΔϵHF+vvvv+cccc-4vvcc+2vccv)))(({{(vccv^2), (Sqrt[2]) (vccc-cvvv) vccv}, {(Sqrt[2]) (vccc-cvvv) vccv, 2 (((vccc-cvvv))^2)}}))-ω IdentityMatrix[2]]==0,ω]
HaToeV ω/.%

Maitra’s kernel

fxcA[ω_,λ_]=((2 (λ^2)(((vccc-cvvv))^2))/(ω-(2ΔϵHF+λ(vvvv+cccc-4vvcc+2vccv))));
fxcB[ω_,λ_]=((2(λ^2)(((vccc-cvvv))^2))/(ω-(2ΔϵHF+λ(vvvv+cccc-4vvcc+2vccv))));

H1[ω_]=(({{ΔϵHF+λ (2vccv-vvcc), λ(2vccv-vccv)}, {-λ(2vccv-vccv), -(ΔϵHF+λ(2vccv-vvcc))}}));
H3[ω_]=(({{ΔϵHF-λ vvcc, -λ  vccv}, {-λ (-vccv), -(ΔϵHF-λ vvcc)}}));
H1TDA[ω_]=(({{ΔϵHF+2vccv-vvcc, 0}, {0, -(ΔϵHF+2vccv-vvcc)}}));
H3TDA[ω_]=(({{ΔϵHF-vvcc, 0}, {0, -(ΔϵHF-vvcc)}}));
H1D[ω_]=(({{ΔϵHF+λ (2vccv-vvcc)+fxcA[ω,λ], λ(2vccv-vccv)+fxcB[ω,λ]}, {-(λ(2vccv-vccv)+fxcB[-ω,λ]), -(ΔϵHF+λ(2vccv-vvcc)+fxcA[-ω,λ])}}));
H1DTDA[ω_]=(({{ΔϵHF+2vccv-vvcc+fxcA[ω,λ], 0}, {0, -(ΔϵHF+2vccv-vvcc+fxcA[-ω,λ])}}));

NSolve[Det[H1[ω]-ω IdentityMatrix[2]]==0/.λ->1,ω]
HaToeV ω/.%
NSolve[Det[H1D[ω]-ω IdentityMatrix[2]]==0/.λ->1,ω]
HaToeV ω/.%
NSolve[Det[H3[ω]-ω IdentityMatrix[2]]==0/.λ->1,ω]
HaToeV ω/.%
NSolve[Det[H1TDA[ω]-ω IdentityMatrix[2]]==0/.λ->1,ω]
HaToeV ω/.%
NSolve[Det[H1DTDA[ω]-ω IdentityMatrix[2]]==0/.λ->1,ω]
HaToeV ω/.%
NSolve[Det[H3TDA[ω]-ω IdentityMatrix[2]]==0/.λ->1,ω]
HaToeV ω/.%

Plot[{
Det[H1[ω]-ω IdentityMatrix[2]]/.λ->1,Det[H1D[ω]-ω IdentityMatrix[2]]/.λ->1,
Det[H3[ω]-ω IdentityMatrix[2]]/.λ->1
},{ω,0,3},
PlotStyle->{{Thick,Gray,Dashed},{Thick,Black},{Thick,Orange,Dashed}},
PlotRange->{Automatic,{-4,10}},ImageSize->Medium,BaseStyle->12,PlotTheme->"Scientific",PlotStyle->Thick,FrameLabel->{"ω","det[H(ω)-ω1]"},
PlotLegends->Placed[{"Singlets (TDHF)","Singlets (D-TDHF)","Triplets (TDHF)"},{Left,Top}]
]

Plot[{
Det[H1BSE-ω IdentityMatrix[2]],Det[H1dBSE[ω]-ω IdentityMatrix[2]],
Det[H3BSE-ω IdentityMatrix[2]],Det[H3dBSE[ω]-ω IdentityMatrix[2]]
},{ω,0,7},
PlotStyle->{{Thick,Gray,Dashed},{Thick,Black},{Thick,Orange,Dashed},{Thick,Red}},
PlotRange->{-30,60},ImageSize->Large,BaseStyle->18,PlotTheme->"Scientific",PlotStyle->Thick,FrameLabel->{"ω","det[H(ω)-ω1]"},
PlotLegends->Placed[{"Singlets (stat.)","Singlets (dyn.)","Triplets (stat.)","Triplets (dyn.)"},{Left,Top}]
]

BSE

Let’s suppose a two-level system with a valence and conduction orbital (HOMO and LUMO basically).
We’re looking for the neutral excitations of the system at the dynamical BSE level enforcing the Tamm-Dancoff approximation.
There’s then one singlet and one triplet singly-excited state, and the BSE Hamiltonian is

WA[ω_]=vvcc+((2 cvvv vccc)/(ω-Ω-ΔϵGW))+((2 cvvv vccc)/(ω-Ω-ΔϵGW));
WB[ω_]=vccv+((2vccv vccv)/(ω-Ω))+((2vccv vccv)/(ω-Ω));

The static BSE Hamiltonian for the singlet states are

H1BSE=(({{ΔϵGW+2vccv-WA[ΔϵGW], 2vccv-WB[0]}, {-(2vccv-WB[0]), -(ΔϵGW+2vccv-WA[ΔϵGW])}}));
H3BSE=(({{ΔϵGW-WA[ΔϵGW], -WB[0]}, {-(-WB[0]), -(ΔϵGW-WA[ΔϵGW])}}));
H1BSETDA=(({{ΔϵGW+2vccv-WA[ΔϵGW], 0}, {0, -(ΔϵGW+2vccv-WA[ΔϵGW])}}));
H3BSETDA=(({{ΔϵGW-WA[ΔϵGW], 0}, {0, -(ΔϵGW-WA[ΔϵGW])}}));

which yields the following excitation energies

Eigenvalues[H1BSE] HaToeV

Eigenvalues[H3BSE]HaToeV

Eigenvalues[H1BSETDA] HaToeV

Eigenvalues[H3BSETDA]HaToeV

The dynamical BSE Hamiltonian is

H1dBSE[ω_]:=(({{ΔϵGW+2vccv-WA[ω], 2vccv-WB[0]}, {-(2vccv-WB[-0]), -(ΔϵGW+2vccv-WA[-ω])}}))
H3dBSE[ω_]:=(({{ΔϵGW-WA[ω], -WB[0]}, {-(-WB[-0]), -(ΔϵGW-WA[-ω])}}))
H1dBSETDA[ω_]:=(({{ΔϵGW+2vccv-WA[ω], 0}, {0, -(ΔϵGW+2vccv-WA[-ω])}}))
H3dBSETDA[ω_]:=(({{ΔϵGW-WA[ω], 0}, {0, -(ΔϵGW-WA[-ω])}}))

HaToeV ω/.NSolve[Det[H1dBSE[ω]-ω IdentityMatrix[2]]==0,ω]//Sort

HaToeV ω/.NSolve[Det[H3dBSE[ω]-ω IdentityMatrix[2]]==0,ω]//Sort

HaToeV ω/.NSolve[Det[H1dBSETDA[ω]-ω IdentityMatrix[2]]==0,ω]//Sort

HaToeV ω/.NSolve[Det[H3dBSETDA[ω]-ω IdentityMatrix[2]]==0,ω]//Sort

Plot[{
Det[H1BSE-ω IdentityMatrix[2]],Det[H1dBSE[ω]-ω IdentityMatrix[2]],
Det[H3BSE-ω IdentityMatrix[2]],Det[H3dBSE[ω]-ω IdentityMatrix[2]]
},{ω,0,4.5},
PlotStyle->{{Thick,Gray,Dashed},{Thick,Black},{Thick,Orange,Dashed},{Thick,Red}},
PlotRange->{-10,30},ImageSize->Medium,BaseStyle->12,PlotTheme->"Scientific",PlotStyle->Thick,FrameLabel->{"ω","det[H(ω)-ω1]"},
PlotLegends->Placed[{"Singlets (BSE)","Singlets (dBSE)","Triplets (BSE)","Triplets (dBSE)"},{Left,Top}]
]

Plot[{
Det[H1BSETDA-ω IdentityMatrix[2]],Det[H1dBSETDA[ω]-ω IdentityMatrix[2]],
Det[H3BSETDA-ω IdentityMatrix[2]],Det[H3dBSETDA[ω]-ω IdentityMatrix[2]]
},{ω,0,5},
PlotStyle->{{Thick,Gray,Dashed},{Thick,Black},{Thick,Orange,Dashed},{Thick,Red}},
PlotRange->{-10,30},ImageSize->Medium,BaseStyle->12,PlotTheme->"Scientific",PlotStyle->Thick,FrameLabel->{"ω","det[H(ω)-ω1]"},
PlotLegends->Placed[{"Singlets [BSE(TDA)]","Singlets [dBSE(TDA)]","Triplets [BSE(TDA)]","Triplets [dBSE(TDA)]"},{Left,Top}]
]

Partitioning for perturbation theory

V1[ω_]:=H1dBSE[ω]-H1BSE;
V3[ω_]:=H3dBSE[ω]-H3BSE;
V1TDA[ω_]:=H1dBSETDA[ω]-H1BSETDA;
V3TDA[ω_]:=H3dBSETDA[ω]-H3BSETDA;

{ω1,Z1}=SortEigensystem[Eigensystem[H1BSE]]
{ω3,Z3}=SortEigensystem[Eigensystem[H3BSE]]

With[{n=2},
{ω,(((1-((Z1[[n]].D[V1[ω],ω].Z1[[n]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)),ω+((((1-((Z1[[n]].D[V1[ω],ω].Z1[[n]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)))((Z1[[n]].V1[ω].Z1[[n]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))}/.ω->ω1[[n]]
]

With[{n=2},
{ω,(((1-((Z3[[n]].D[V3[ω],ω].Z3[[n]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)),ω+((((1-((Z3[[n]].D[V3[ω],ω].Z3[[n]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)))((Z3[[n]].V3[ω].Z3[[n]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))}/.ω->ω3[[n]]
]

HaToeV With[{n=2},
{ω,(((1-((Z1[[n,1]]D[V1[ω][[1,1]],ω]Z1[[n,1]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)),ω+((((1-((Z1[[n,1]]D[V1[ω][[1,1]],ω]Z1[[n,1]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)))((Z1[[n,1]]V1[ω][[1,1]]Z1[[n,1]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))}/.ω->ω1[[n]]
]

HaToeV With[{n=2},
{ω,(((1-((Z3[[n,1]]D[V3[ω][[1,1]],ω]Z3[[n,1]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)),ω+((((1-((Z3[[n,1]]D[V3[ω][[1,1]],ω]Z3[[n,1]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)))((Z3[[n,1]]V3[ω][[1,1]]Z3[[n,1]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))}/.ω->ω3[[n]]
]

BSE2

H1BSE2=(({{ΔϵGF2+2vccv-vvcc, 2vccv-vccv}, {-(2vccv-vccv), -(ΔϵGF2+2vccv-vvcc)}}));
H3BSE2=(({{ΔϵGF2-vvcc, -vccv}, {-(-vccv), -(ΔϵGF2-vvcc)}}));
H1BSE2TDA=(({{ΔϵGF2+2vccv-vvcc, 0}, {0, -(ΔϵGF2+2vccv-vvcc)}}));
H3BSE2TDA=(({{ΔϵGF2-vvcc, 0}, {0, -(ΔϵGF2-vvcc)}}));

H1dBSE2[ω_]=(({{ΔϵGF2+2vccv-vvcc, 2vccv-vccv}, {-(2vccv-vccv), -(ΔϵGF2+2vccv-vvcc)}}))+(({{-((4cvvv vccc-(vccc^2)-(cvvv^2))/(ω-2ΔϵGF2)), +((4(vccv^2)- cccc vccv- vvvv vccv)/(2ΔϵGF2))}, {-((4(vccv^2)- cccc vccv- vvvv vccv)/(2ΔϵGF2)), +((4cvvv vccc-(vccc^2)-(cvvv^2))/(-ω-2ΔϵGF2))}}));
H3dBSE2[ω_]=(({{ΔϵGF2-vvcc, -vccv}, {-(-vccv), -(ΔϵGF2-vvcc)}}))+(({{-(((vccc^2)+(cvvv^2))/(ω-2ϵcGF2+2ϵvGF2)), +(( cccc vccv+vvvv vccv)/(2ϵcGF2-2ϵvGF2))}, {-(( cccc vccv+vvvv vccv)/(2ϵcGF2-2ϵvGF2)), +(((vccc^2)+(cvvv^2))/(-ω-2ϵcGF2+2ϵvGF2))}}));
H1dBSE2TDA[ω_]=(({{ΔϵGF2+2vccv-vvcc, 0}, {0, -(ΔϵGF2+2vccv-vvcc)}}))+(({{-((4cvvv vccc-(vccc^2)-(cvvv^2))/(ω-2ΔϵGF2)), 0}, {0, +((4cvvv vccc-(vccc^2)-(cvvv^2))/(-ω-2ΔϵGF2))}}));
H3dBSE2TDA[ω_]=(({{ΔϵGF2-vvcc, 0}, {0, -(ΔϵGF2-vvcc)}}))+(({{-(((vccc^2)+(cvvv^2))/(ω-2ϵcGF2+2ϵvGF2)), 0}, {0, +(((vccc^2)+(cvvv^2))/(-ω-2ϵcGF2+2ϵvGF2))}}));

Plot[{
Det[H1dBSE2-ω IdentityMatrix[2]],Det[H1dBSE2[ω]-ω IdentityMatrix[2]],
Det[H3dBSE2-ω IdentityMatrix[2]],Det[H3dBSE2[ω]-ω IdentityMatrix[2]]
},{ω,0,6},
PlotStyle->{{Thick,Gray,Dashed},{Thick,Black},{Thick,Orange,Dashed},{Thick,Red}},
PlotRange->{-10,40},ImageSize->Medium,BaseStyle->12,PlotTheme->"Scientific",PlotStyle->Thick,FrameLabel->{"ω","det[H(ω)-ω1]"},
PlotLegends->Placed[{"Singlets (BSE2)","Singlets (dBSE2)","Triplets (BSE2)","Triplets (dBSE2)"},{Left,Top}]
]

NSolve[Det[H1BSE2-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%
NSolve[Det[H3BSE2-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%

NSolve[Det[H1dBSE2[ω]-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%
NSolve[Det[H3dBSE2[ω]-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%

NSolve[Det[H1BSE2TDA-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%
NSolve[Det[H3BSE2TDA-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%

NSolve[Det[H1dBSE2TDA[ω]-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%
NSolve[Det[H3dBSE2TDA[ω]-ω IdentityMatrix[2]]==0]
HaToeV  ω/.%

Partitioning for perturbation theory

V1[ω_]:=H1dBSE2[ω]-H1BSE2;
V3[ω_]:=H3dBSE2[ω]-H3BSE2;
V1TDA[ω_]:=H1dBSE2TDA[ω]-H1BSE2TDA;
V3TDA[ω_]:=H3dBSE2TDA[ω]-H3BSE2TDA;

{ω1,Z1}=SortEigensystem[Eigensystem[H1BSE2]]
{ω3,Z3}=SortEigensystem[Eigensystem[H3BSE2]]

With[{n=2},
{ω,(((1-((Z1[[n]].D[V1[ω],ω].Z1[[n]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)),ω+((((1-((Z1[[n]].D[V1[ω],ω].Z1[[n]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)))((Z1[[n]].V1[ω].Z1[[n]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))}/.ω->ω1[[n]]
]

With[{n=2},
{ω,(((1-((Z3[[n]].D[V3[ω],ω].Z3[[n]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)),ω+((((1-((Z3[[n]].D[V3[ω],ω].Z3[[n]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)))((Z3[[n]].V3[ω].Z3[[n]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))}/.ω->ω3[[n]]
]

HaToeV With[{n=2},
{ω,(((1-((Z1[[n,1]]D[V1[ω][[1,1]],ω]Z1[[n,1]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)),ω+((((1-((Z1[[n,1]]D[V1[ω][[1,1]],ω]Z1[[n,1]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))))^(-1)))((Z1[[n,1]]V1[ω][[1,1]]Z1[[n,1]])/(((Z1[[n,1]])^2)-((Z1[[n,2]])^2)))}/.ω->ω1[[n]]
]

HaToeV With[{n=2},
{ω,(((1-((Z3[[n,1]]D[V3[ω][[1,1]],ω]Z3[[n,1]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)),ω+((((1-((Z3[[n,1]]D[V3[ω][[1,1]],ω]Z3[[n,1]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))))^(-1)))((Z3[[n,1]]V3[ω][[1,1]]Z3[[n,1]])/(((Z3[[n,1]])^2)-((Z3[[n,2]])^2)))}/.ω->ω3[[n]]
]

Sangalli’s kernel (I think)

H1[ω_]=(({{ΔϵHF+((((cvvv+vccc))^2)/(ω-2Ω1)), 0}, {0, -(ΔϵHF+((((cvvv+vccc))^2)/(-ω-2Ω1)))}}));

NSolve[Δϵ+2((((cvvv+vccc))^2)/(ω-2 Ω1))-ω==0,ω]

Plot[{Δϵ+2((((cvvv+vccc))^2)/(ω-2 Ω1))-ω},{ω,0,7},PlotRange->{-15,10},ImageSize->Large,BaseStyle->20,PlotTheme->"Scientific",PlotStyle->Thick,FrameLabel->{"ω","det[(H^(dBSE (TDA)))(ω)-ω I]"},PlotLegends->Placed[{"Singlets","Triplets"},{Left,Top}]]