He 1(s^2) Perturbation Theory

(SuperPlus[He]) wavefunctions, atomic units

(Subscript[P, n_,l_])[r_]:=(Sqrt[(((n-l-1)! 2)/((n^2)((n+l)!)))])(((((4 r)/n)))^(l+1))(E^(-((2r)/n)))LaguerreL[n-l-1,2l+1,((4r)/n)]

Matrix element for L=0,  will need to be modified for the 1s2p states with L=1.

V[{n1_,n2_,l1_},{l_,lp_},{n3_,n4_,l4_}]:=
Sum[If[(a=(((-1))^ll)(Sqrt[((2lp+1)/(2l+1))])(((ClebschGordan[{ll,0},{lp,0},{l,0}]))^2))==0,0,
a NIntegrate[(Subscript[P, n1,l1])[r1](Subscript[P, n3,l4])[r1](Subscript[P, n2,l1])[r2](Subscript[P, n4,l4])[r2]Min[((r1^ll)/(r2^(ll+1))),((r2^ll)/(r1^(ll+1)))],{r2,0,∞},{r1,0,∞},AccuracyGoal->5]],{ll,Abs[l-lp],l+lp}]

Calculation of unperturbed energies.

E1[{n1_,n2_,l_}]:=((-2)/(n1^2))+((-2)/(n2^2))+V[{n1,n2,l},{l,l},{n1,n2,l}]

Quantum numbers for state |0⟩

s0={1,1,0};

E1[{1,1,0}]

List of states to include in the calculation.  Start small and increase until your answer stops changing.

st=Partition[Flatten[Table[(Piecewise[{{{n1,n2,l}, n1==n2}, {{{n1,n2,l},{n2,n1,l}}, True}}]),{l,0,2},{n1,l+1,11},{n2,n1,11}]],3];st=Delete[st,1];ls=st\[Transpose][[3]];

Calculation of energy to second order.  Don’t worry if the numerical integration gives “slwcon” errors.

E1[s0]+Sum[(((V[s0,{0,ls[[i]]},st[[i]]])^2)/(E1[s0]-E1[st[[i]]])),{i,Length[st]}]

%/-2.903