Step Down Operator Definition

Clear[ToLstFrm,SingleStepDown,CleanUpDeterm,l1,l2,l3,l4];
StepDown:=Module[{i,j,k,l=3,s=1/2,l1={},l2={},l3={},l4={},np=0,ToLstFrm,SingleStepDown,CleanUpDeterm,PermuteDeterm},
ToLstFrm:=
(If[MatchQ[Head[#],Plus],
List@@#,
{#}]//
(For[i=1,i≤Length[#],i++,
If[MatchQ[Head[#[[i]]],SD],{AppendTo[l1,1],AppendTo[l2,#[[i]]]},
{AppendTo[l1,1],AppendTo[l2,1],For[j=1,j≤Length[#[[i]]],j++,
If[NumericQ[#[[i,j]]],{l1[[i]]=l1[[i]]*#[[i,j]]},{l2[[i]]=#[[i,j]]}]]}]])&)&;

SingleStepDown:=
Sum[Sqrt[(l+#[[i]])(l-#[[i]]+1)]Apply[Head[#],{#[[1]]-KroneckerDelta[i,1],#[[2]]-KroneckerDelta[i,2],#[[3]]-KroneckerDelta[i,3],#[[4]],#[[5]],#[[6]]}],{i,1,3}]+Sum[Sqrt[(s+#[[i]])(s-#[[i]]+1)]Apply[Head[#],{#[[1]],#[[2]],#[[3]],#[[4]]-KroneckerDelta[i,4],#[[5]]-KroneckerDelta[i,5],#[[6]]-KroneckerDelta[i,6]}],{i,4,6}]&;

CleanUpDeterm:=
For[k=1,k≤Length[l2],k++,
For[i=1,i≤3,i++,
For[j=i+1,j≤3,j++,
If[l2[[k,i]]==l2[[k,j]]&&l2[[k,i+3]]==l2[[k,j+3]],l1[[k]]=0]]]];

(*to be run on a list of determinants, do ToLstFrm on it first*)
PermuteDeterm:=
For[k=1,k≤Length[l2],k++,
For[i=1,i≤3,i++,(*assume that three runs are enough (and it is :) *)
{
If[l2[[k,1]]<l2[[k,2]],{{l2[[k,1]],l2[[k,2]],l2[[k,4]],l2[[k,5]]}={l2[[k,2]],l2[[k,1]],l2[[k,5]],l2[[k,4]]};l1[[k]]=-l1[[k]]},
If[l2[[k,1]]==l2[[k,2]],
If[l2[[k,4]]<l2[[k,5]],{{l2[[k,4]],l2[[k,5]]}={l2[[k,5]],l2[[k,4]]};l1[[k]]=-l1[[k]]}]]];
If[l2[[k,2]]<l2[[k,3]],{{l2[[k,2]],l2[[k,3]],l2[[k,5]],l2[[k,6]]}={l2[[k,3]],l2[[k,2]],l2[[k,6]],l2[[k,5]]};l1[[k]]=-l1[[k]]},
If[l2[[k,2]]==l2[[k,3]],
If[l2[[k,5]]<l2[[k,6]],{{l2[[k,5]],l2[[k,6]]}={l2[[k,6]],l2[[k,5]]};l1[[k]]=-l1[[k]]}]]];
}]];

ToLstFrm[#];
For[j=1,j≤Length[l2],j++,l2[[j]]=SingleStepDown[l2[[j]]]];
l3=l1;l4=l2;l1={};l2={};

ToLstFrm[Expand[l3.l4]];
CleanUpDeterm;
PermuteDeterm;

l1.l2
]&;

State Determination

St={SD[3,2,1,1/2,1/2,1/2]};
JMFact[J_,M_]=1/Sqrt[(J+M)(J-M+1)];
For[i=0,i≤14,i++,
{AppendTo[St,Expand[JMFact[15/2,15/2-i]StepDown[ St[[i+1]] ]]]} ]

Matrix Element Calculation

MatElem[S1_,S2_,LV_,MV_,KV_]:=
Module[{i,j,k,st1,st2,SDs1,SD1,Fcs1,SDs2,SD2,Fcs2,SDpr,df1,df2,perm,skd,ints,res,ToLstFrm,NumDiff,SnglMatElem},

(*=== Local Modules ===*)
MatElem::argtype="Determinants with unequal number of arguments not allowed.";

(*transform algebraic sum to list form*)
ToLstFrm:=If[MatchQ[Head[#],Plus],List@@#,{#}]&;

(* determines number of different spin orbitals in two Slater determinants and the positions of them within the determinant *)
NumDiff[sd1_,sd2_]:=
Module[{i,n=0,ls={}},
If[Length[sd1]≠Length[sd2],Message[MatElem::argtype]];
For[i=1,i≤3,i++,If[sd1[[i]]≠sd2[[i]]||sd1[[i+3]]≠sd2[[i+3]],{n++,AppendTo[ls,i]}]];
Return[Prepend[ls,n]];
];

(* matrix element between two spin orbitals *)
SnglMatElem[l1_,m1_,l2_,m2_,lv_,mv_,kv_]:=(-1)^m1 Sqrt[(2l1+1)(2lv+1)(2l2+1)/(4Pi)]ThreeJSymbol[{l1,0},{lv,0},{l2,0}]ThreeJSymbol[{l1,-m1},{lv,mv},{l2,m2}]F[kv];


(*=== Main Module Body ===*)
st1=ToLstFrm[Expand[S1]];
st2=ToLstFrm[Expand[S2]];

SDs1=Extract[st1,Position[st1,SD[__]]];
Fcs1=ReplacePart[st1,Position[st1,SD[__]]->1];
SDs2=Extract[st2,Position[st2,SD[__]]];
Fcs2=ReplacePart[st2,Position[st2,SD[__]]->1];

res=0;
For[i=1,i≤Length[SDs1],i++,
For[j=1,j≤Length[SDs2],j++,

SD1={{SDs1[[i,1]],SDs1[[i,4]]},{SDs1[[i,2]],SDs1[[i,5]]},{SDs1[[i,3]],SDs1[[i,6]]}};
SD2={{SDs2[[j,1]],SDs2[[j,4]]},{SDs2[[j,2]],SDs2[[j,5]]},{SDs2[[j,3]],SDs2[[j,6]]}};
ints=Intersection[SD1,SD2];

skd=Sum[KroneckerDelta[SD1[[k,1]]-SD2[[k,1]]],{k,1,3}];
If[Length[ints]==0,Continue[]];
If[Length[ints]==1,Continue[]];
If[Length[ints]==2,
df1=Delete[SD1,{Position[SD1,ints[[1]]][[1]],Position[SD1,ints[[2]]][[1]]}];
df2=Delete[SD2,{Position[SD2,ints[[1]]][[1]],Position[SD2,ints[[2]]][[1]]}];

If[skd==2,perm=0];
If[skd==1,perm=1];
If[skd==0,
If[Position[SD1,df1[[1]]]==Position[SD2,df2[[1]]],perm=1,perm=2];
];

res+=(-1)^perm KroneckerDelta[ df1[[1,2]],df2[[1,2]] ]Fcs1[[i]]Fcs2[[j]]SnglMatElem[3,df1[[1,1 ]],3,df2[[1,1]],LV,MV,KV];
];
If[Length[ints]==3,
If[skd==3,perm=0];
If[skd==2,perm=0];
If[skd==1,perm=1];
If[skd==0,perm=2];

For[k=1,k≤3,k++,
res+=(-1)^perm KroneckerDelta[ SD1[[k,2]],SD2[[k,2]] ]Fcs1[[i]]Fcs2[[j]]SnglMatElem[3,SD1[[k,1]],3,SD2[[k,1]],LV,MV,KV];
];
];
];
];
Return[res];
];

Off[ClebschGordan::"phy"]
mat0=Table[Simplify[(*12Sqrt[(Pi/11)]*)MatElem[Expand[St[[i]]],St[[j]],0,0,0]],{i,1,Length[St]},{j,1,Length[St]}];
mat2=Table[Simplify[-I(*(2/121)Sqrt[((30Pi)/11)]*)MatElem[Expand[St[[i]]],St[[j]],2,-2,2]],{i,1,Length[St]},{j,1,Length[St]}];
mat2-=Table[Simplify[-I(*(2/121)Sqrt[((30Pi)/11)]*)MatElem[Expand[St[[i]]],St[[j]],2,2,2]],{i,1,Length[St]},{j,1,Length[St]}];

Simplify[mat2//MatrixForm]

Determination of Eigenenergies

Dat1=Import["G:/Documents/AAU/Mathematica/HF_ExcStates/EigV_Er_3+_4f_U.txt","List"];
Dat2 =Import["G:/Documents/AAU/Mathematica/HF_ExcStates/Bas_Er_f.txt","List"];
EigV=Table[Dat1[[i]],{i,4,Dat1[[3]]+3}];
Bas=Table[Dat2[[i]],{i,3,Dat2[[2]]+2}];

l=3;k=0;
RadInt0=Sum[Sum[EigV[[i]]EigV[[j]]Integrate[r^(2l+k+2)Exp[-(Bas[[i]]+Bas[[j]])r^2],{r,0,∞}],{i,1,Length[Bas]}],{j,1,Length[Bas]}];
l=3;k=2;
RadInt2=Sum[Sum[EigV[[i]]EigV[[j]]Integrate[r^(2l+k+2)Exp[-(Bas[[i]]+Bas[[j]])r^2],{r,0,∞}],{i,1,Length[Bas]}],{j,1,Length[Bas]}];

Sort[Eigenvalues[(mat0/.F[0]->RadInt0)*0+(mat2/.F[2]->1)]//N]

Testing

StepDown[2 (Sqrt[3]) SD[3,2,0,(1/2),(1/2),(1/2)]+SD[3,2,1,-(1/2),(1/2),(1/2)]+SD[3,2,1,(1/2),-(1/2),(1/2)]+SD[3,2,1,(1/2),(1/2),-(1/2)]]

St[[16]]

St[[16]]

SnglMatElem[l1_,m1_,l2_,m2_,lv_,mv_,kv_]:=(-1)^m1 Sqrt[(2l1+1)(2lv+1)(2l2+1)/(4Pi)]ThreeJSymbol[{l1,0},{lv,0},{l2,0}]ThreeJSymbol[{l1,-m1},{lv,mv},{l2,m2}]F[kv];

(1/I)(Table[SnglMatElem[3,i,3,j,2,-2,2],{i,-3,3},{j,-3,3}]-Table[SnglMatElem[3,i,3,j,2,2,2],{i,-3,3},{j,-3,3}])//MatrixForm

(({{0, 0, ((I F[2])/(3 (Sqrt[2 π]))), 0, 0, 0, 0}, {0, 0, 0, -((I F[2])/(3 (Sqrt[π]))), 0, 0, 0}, {-((I F[2])/(3 (Sqrt[2 π]))), 0, 0, 0, I (Sqrt[(2/(15 π))]) F[2], 0, 0}, {0, ((I F[2])/(3 (Sqrt[π]))), 0, 0, 0, -((I F[2])/(3 (Sqrt[π]))), 0}, {0, 0, -I (Sqrt[(2/(15 π))]) F[2], 0, 0, 0, ((I F[2])/(3 (Sqrt[2 π])))}, {0, 0, 0, ((I F[2])/(3 (Sqrt[π]))), 0, 0, 0}, {0, 0, 0, 0, -((I F[2])/(3 (Sqrt[2 π]))), 0, 0}}))

St[[3]]

2Sqrt[(3/35)](SnglMatElem[3,-1,3,1,2,-2,2]-SnglMatElem[3,-1,3,1,2,2,2])-Sqrt[(2/7)](SnglMatElem[3,0,3,2,2,-2,2]-SnglMatElem[3,0,3,2,2,2,2])

SnglMatElem[3,-1,3,1,2,2,2]-SnglMatElem[3,-1,3,1,2,-2,2]
SnglMatElem[3,0,3,2,2,2,2]-SnglMatElem[3,0,3,2,2,-2,2]

Simplify[MatElem[St[[4]],St[[2]],2,-2,2]]

St[[4]]
St[[2]]

MatElem[((2 SD[3,2,0,(1/2),(1/2),(1/2)])/(Sqrt[5])),((2 SD[2,1,0,(1/2),(1/2),(1/2)])/(Sqrt[91]))+4 (Sqrt[(2/91)]) SD[3,1,-1,(1/2),(1/2),(1/2)]*0+2 (Sqrt[(2/91)]) SD[3,2,-2,(1/2),(1/2),(1/2)]*0,2,-2,2]-MatElem[((2 SD[3,2,0,(1/2),(1/2),(1/2)])/(Sqrt[5])),((2 SD[2,1,0,(1/2),(1/2),(1/2)])/(Sqrt[91]))+4 (Sqrt[(2/91)]) SD[3,1,-1,(1/2),(1/2),(1/2)]*0+2 (Sqrt[(2/91)]) SD[3,2,-2,(1/2),(1/2),(1/2)]*0,2,2,2]

St[[1]]

Simplify[MatElem[St[[4]],St[[2]],2,-2,2]-MatElem[St[[4]],St[[2]],2,2,2]]

{{0,(1/2)},{2,(1/2)},{1,(1/2)}}[[1,1]]

Testing2

St[[1]]
St[[4]]

Simplify[MatElem[St[[1]],St[[1]],0,0,0]]

Integrate[Integrate[SphericalHarmonic

(mat0/.F[0]->1)//MatrixForm

SnglMatElem[l1_,m1_,l2_,m2_,lv_,mv_,kv_]:=(-1)^m1 Sqrt[(2l1+1)(2lv+1)(2l2+1)/(4Pi)]ThreeJSymbol[{l1,0},{lv,0},{l2,0}]ThreeJSymbol[{l1,-m1},{lv,mv},{l2,m2}]F[kv];
Integrate[Integrate[SphericalHarmonicY[3,3,t,p]\[Conjugate]SphericalHarmonicY[3,3,t,p]Sin[t],{t,0,Pi}],{p,0,2π}]
SnglMatElem[3,3,3,3,0,0,0]

Three