Multiplets for non-equivalent electrons

Quantum Mechanics
                     K. Konishi, G. Paffuti

Compute the possible multiplets for a shell of equivalent electrons and one electron in another shell. Study the diagonalization of the Hamiltonian.

Warning: This notebook has several global variables. It is better to restart from this point if some parameters are changed. It is suggested  to quit and restart the Kernel in the first use of this program.

ClearAll["Global`*"];
Off[General::spell, General::spell1];
<<Combinatorica`

Description of the method

This notebook is quite similar to notebook [NB3].  It is assumed that the reader has already acquainted with that notebook. The comments and some instructions specific to variables of the second shell are emphasized with a yellow background.

Let us consider a number numE of electrons with angular momentum Lorb. Starting with these orbitals a certain number of multiplets, with given total angular momentum L and total spin S, can be constructed. Each (antisymmetric) state with given quantum numbers for each electron is represented in x-space by a Slater determinant.  The number of different states which can be constructed with numE equivalent electrons is given by the number of combinations of numE elements constructed with the 2 (2 Lorb +1) single particle states, C(2 (2 Lorb+1); numE). These states are collected in linear combinations with given L and S. S is integer or semi-integer if numE is even or odd. Each multiplet has a degeneration (2L+1)(2S+1) and the sum of the degeneracies must reproduce the total number of states Σ (2L+1)(2S+1) = C(2 (2 Lorb +1); numE).

In each multiplet states can be classified by the quantum numbers L,S,(Subscript[L, z]),(Subscript[S, z]). Each multiplet has a state with (Subscript[L, z]) = 0 , (Subscript[S, z])=0 or (Subscript[L, z]) = 0 , (Subscript[S, z])=1/2, following the parity of numE. To classify multiplets is sufficient to compute the eigenvalues of (L^2) and (S^2) on this subset of states.

To these states, at variance with  [NB3] we add states for a non equivalent electron, all the variable for this particle have the suffix NE.

Notation

We will use the notation

v[{(Subscript[m, 1]),(Subscript[S, 1]),(Subscript[m, 2]),(Subscript[S, 2]),... ,(Subscript[m, k]),(Subscript[S, k])}]

a state with fixed projections along z for L and S of each electron, i.e. a Slater determinant. (Subscript[m, i]),(Subscript[S, i])are the quantum numbers of the i-th electron. The list in v[] is ordered and this fixes the sign of the determinant. The quantum numbers are understood completely antisymmetric, in agreement with the Pauli principle.

The last two entries refer to the non equivalent electron.

Input data

The program needs the parameters numE and Lorb. mLtot and mStot stand for the selected values of total orbital momentum and spin along z, following the above criteria.

(* numE is the number of equivalent electrons, with angular momentum Lorb *)
numE=2; Lorb =2;
 LorbNE =1; (* the angular momentum of the non-equivalent electron *)
 Lmin = Abs[Lorb-LorbNE]; Lmax = Lorb+LorbNE;
 mLtot=0;mStot=Mod[numE +1,2]/2;

The definition of some conventional names for states

lStati = {"S","P","D","F","G","H","I","K","L","M","N"};
nomorb={"s","p","d","f"};
nomeNE=If[LorbNE== Lorb, StringJoin[nomorb[[LorbNE+1]],"'"],nomorb[[LorbNE+1]]];
str1 = nomorb[[Lorb+1]];
str2=StringJoin[nomorb[[Lorb+1]],nomeNE];

The generation of the basis

Let us call numq a list with the quantum numbers of a single particle. We have -Lorb ≤ lz ≤ Lorb  and  sz = ± 1/2:

numq = Flatten[Table[{i,j},{i,-Lorb,Lorb},{j,-1/2,1/2}],1]

another list numqNE is generated for the non-equivalent electron.

numqNE = Flatten[Table[{i,j},{i,-LorbNE,LorbNE},{j,-1/2,1/2}],1]

To construct a basis we use the Mathematica function KSubsets[ list, k]  which generate all subsets of length k from a list, see examples in  [NB3].

For an instruction as part = KSubsets[Range[Length[numq]], numE] the element part[[i,j]] gives the component j of the i-th partition. numq[[ part[i,j], 1]] will be its angular momentum and numq[[ part[i,j], 2]] its spin

To extract the basis we compute the total Lz and Sz, if these satisfy the constraint the element will be added to the base.

At the end of a loop the function makeOrthogonal construct a vector indice composed of 0,1 which shows if a given state is filled in a given Slater determinant. This vector will be useful in the computation of  matrix elements of the Hamiltonian.

faibase[numq_,numqNE_,numE_,mLtot_,mStot_]:=
Module[{base,indice,iE,iv,index,m1,s1,i1},
(base={};indice={};
nL=Length[numq];
nLNE = Length[numqNE];
part=KSubsets[Range[Length[numq]],numE];
m1 = Table[0,{k,1,numE}];s1 = Table[0,{k,1,numE}];
i1 = Table[0,{k,1,numE}];
Do[
(Table[m1[[k]] = numq[[part[[iv,k]],1]],{k,1,numE}];
Table[s1[[k]] = numq[[part[[iv,k]],2]],{k,1,numE}];
mNE = numqNE[[iNE,1]];sNE = numqNE[[iNE,2]];
If[ (Total[m1]+mNE == mLtot) && (Total[s1]+sNE==mStot),
(index=Table[0,{kk,1,Length[numq]}];
index[[part[[iv]]]]=1;
indice=Append[indice,index];
 base = Append[base,v[Flatten[{numq[[part[[iv]]]],numqNE[[iNE]]}]]];),];),
{iv,1,Length[part]},{iNE,1,nLNE}]);
Return[{base,indice}]]

{base,indice} = faibase[numq,numqNE,numE,mLtot,mStot];
Length[base]

Example

{base[[1]],indice[[1]]}

Our states are assumed to be ordered and completely antisymmetric. In some cases below it is possible that some non ordered vectors are generated, these can be set in the canonical order by a permutation, and we must take into account the sign of this permutation. This is done with the function XXX. The call to the function Partition is made to construct sublists of length 2 (a single particle) which will be ordered in the canonical form

orderVector[v[kin_]]:=Module[{k,klast,k1,signk},(
k = Drop[kin,-2];klast = Take[kin,-2];
k1 = Sort[Partition[k,2]];lk1=Length[k1];
testk1 = MemberQ[k1[[Range[1,lk1-1]]]-k1[[Range[2,lk1]]],{0,0}];
If[testk1,0,(k1=Partition[k,2]; signk = Signature[k1];
signk*v[Flatten[Join[Sort[k1],klast]]])])];
orderVector[0]=0;

Angular momentum and Spin

We define rising and lowering operators for orbital angular momentum and spin. As an instance  lp is the rising operator for the base vector. We extend the definition to linear combinations of base vectors by the Mathematica function Distribute.

In the definition of lp and other operators we use the instruction x_. v[k_] The point means that without the argument (x in this case) the variable acquire its default value, which is 1 for the product . With this notation we will have lp[ 2 v] = 2 lp[v] while  lp[v] will be computed as lp[1 * v].
Due to our notation for the vectors v[] the index k[[i]] in the sums below will select angular momentum or spin in odd/even places of v[].
The function ReplacePart[ k, new, i] substitutes the i-th value of the list k with the value new.

lz[x_.  v[k_]]:= x* Sum[k[[i]]*v[k],{i,1,Length[k]-1,2}]
lz[0]=0; lp[0,Lorb]=0;lm[0,Lorb]=0;
lm[x_.  v[k_],Lorb_,LorbNE_]:= (
x*Sum[(Sqrt[Lorb*(Lorb+1) - k[[i]](k[[i]]-1)])*
orderVector[v[ReplacePart[k,k[[i]]-1,i]]],{i,1,Length[k]-3,2}]+
 x*(Sqrt[LorbNE*(LorbNE+1) - k[[-2]](k[[-2]]-1)])v[ReplacePart[k,k[[-2]]-1,-2]])
lm[0]=0;
lp[x_.  v[k_],Lorb_,LorbNE_]:= (
x*Sum[(Sqrt[Lorb*(Lorb+1) - k[[i]](k[[i]]+1)])*
orderVector[v[ReplacePart[k,k[[i]]+1,i]]],{i,1,Length[k]-3,2}]+
x*(Sqrt[LorbNE*(LorbNE+1) - k[[-2]](k[[-2]]+1)])v[ReplacePart[k,k[[-2]]+1,-2]])
lp[0]=0;
Lp[a_,Lorb_,LorbNE_]:= Expand[Distribute[lp[a,Lorb,LorbNE]]];
Lm[a_,Lorb_,LorbNE_]:= Expand[Distribute[lm[a,Lorb,LorbNE]]];
Lz[a_]:= Expand[Distribute[lz[a]]];
L2[a_,Lorb_,LorbNE_]:= Lp[Lm[a,Lorb,LorbNE],Lorb,LorbNE] + Lz[Lz[a]] - Lz[a];
Lp[0,Lorb_,LorbNE_]:= 0;

sz[x_.  v[k_]]:= x* Sum[k[[i]]*v[k],{i,2,Length[k],2}]
sz[0]=0;
sm[x_.  v[k_]]:=( x*Sum[(Sqrt[(1/2)*((1/2)+1) - k[[i]](k[[i]]-1)])*
orderVector[v[ReplacePart[k,k[[i]]-1,i]]],{i,2,Length[k]-2,2}]+
x*(Sqrt[(1/2)*((1/2)+1) - k[[-1]](k[[-1]]-1)])v[ReplacePart[k,k[[-1]]-1,-1]])
sm[0]:=0;
sp[x_.  v[k_]]:= (x*Sum[(Sqrt[(1/2)*((1/2)+1) - k[[i]](k[[i]]+1)])*
orderVector[v[ReplacePart[k,k[[i]]+1,i]]],{i,2,Length[k]-2,2}]+
x*(Sqrt[(1/2)*((1/2)+1) - k[[-1]](k[[-1]]+1)])v[ReplacePart[k,k[[-1]]+1,-1]])
sp[0]=0;
Sp[a_]:= Expand[Distribute[sp[a]]]
Sm[a_]:= Expand[Distribute[sm[a]]]
Sz[a_]:= Expand[Distribute[sz[a]]]
S2[a_]:= Sp[Sm[a]] + Sz[Sz[a]] - Sz[a]

In the following we will produce a basis of eigenvectors for (L^2) ed (S^2) . For degenerate eigenvalues the basis is not orthogonal in general. The function  makeOrthogonal provide the orthogonalization procedure. For speed is written to act on single subspace with fixed eigenvalue. This function works only if the difference between eigenvalues is an integer.

makeOrthogonal[eigvL_,autovL_]:=Module[{eigvLnew,newBdeg,degL},
(eigvLnew = Table[0,{k,1,Length[eigvL]}];
Do[
(degL={};
Table[If[autovL[[k]]==i,degL=Append[degL,k],],{k,1,Length[autovL]}];
If[Length[degL]≥ 1,newBdeg = Simplify[Orthogonalize[Table[eigvL[[k]],{k,First[degL],Last[degL]}]]];
eigvLnew[[degL]] = newBdeg;]),{i,Max[autovL],Min[autovL],-1}];
eigvLnew)]

Diagonalization of L and S

Operator L2

Let us compute the matrix cL defined by L2[ (Subscript[b, i])] = (Subscript[C, ki ])(Subscript[b, k]) ;  (Subscript[b, i]) are basis vectors.  Using Eigensystem to diagonalize the matrix one obtains a matrix whose column are the component of eigenvectors. Let call Eig this matrix.  The new eigenvectors are (Subscript[F, i])= (Subscript[Eig, ij])(Subscript[b, j]) (below we give a check on the order of indexes).
The matrix cL can be computed extracting the coefficient of base[[k]] from L2[ base[[i]] ]

cL=Table[Coefficient[Table[L2[base[[i]],Lorb,LorbNE],{i,1,Length[base]}],base[[k]]],{k,1,Length[base]}];

Here we list the values for L, related to the eigenvalues λ of (L^2)by L = (Sqrt[λ +1/4])-1/2:

{autovL2primabase,eigvL}=Eigensystem[cL];
autovLprimabase = Sqrt[autovL2primabase + 1/4] - 1/2

We construct an orthonormal basis using the function makeOrthogonal

eigvL=makeOrthogonal[eigvL,autovLprimabase];

Operator S2

Let call S2 the matrix representing (S^2) in the original basis S2[ (Subscript[b, i])] = (Subscript[S, ji])(Subscript[b, j]). We want to diagonalize the operator S2 in the basis in which L2 is also diagonal, in this way we can have a simultaneous diagonalization for both operators. We can switch to the basis F = eigvL with

S2[(Subscript[Eig, ji])(Subscript[b, i])]=(Subscript[Eig, ji]) (Subscript[S, si])(Subscript[((Eig^(-1))), sk])(Subscript[F, k])

Then with respect to the basis eigvL of the operator (L^2):

S2[(Subscript[F, j])]=(Subscript[((Eig^(T-1))), ks])(Subscript[S, si])(Subscript[((Eig^T)), ij]) (Subscript[F, k]) =(Subscript[(Eig), ks])(Subscript[S, si])(Subscript[((Eig^T)), ij])(Subscript[F, k])       (Eig is an orthogonal basis)

The expression of the new matrix for S2 is  Eig S2 (Eig^T).
Let eigvSconL the new eigenvectors, in the old representation they are given by:

eigvLS=eigvSconL.eigvL

First we compute the matrix S2 in the original basis again by extracting the coefficient of base[[k]] in  S2[base[[i]]]=S2_ki b_k

cSbase=Table[Coefficient[Table[Expand[Simplify[S2[base[[i]]]]],{i,1,Length[base]}],base[[k]]],{k,1,Length[base]}];

Here we find the eigenvalues and the eigenvectors with respect to the basis eigvL (eigenvectors for (L^2))

{autovS2,eigvSconL}=Eigensystem[Simplify[eigvL.cSbase.Transpose[eigvL]]];
autovS = Sqrt[autovS2 + 1/4] - 1/2

We diagonalize the basis

eigvSconL = makeOrthogonal[eigvSconL,autovS2];

and finally write the eigenvectors of (L^2)and (S^2) in the original basis

eigvLS = Simplify[eigvSconL.eigvL];

Let us check that the new vectors are eigenvectors of  L2 by computing the mean value of L2 on them:

autovL2=Expand[Table[eigvLS[[i]].Simplify[cL.eigvLS[[i]]],{i,1,Length[base]}]];
autovL = Sqrt[autovL2 + 1/4] - 1/2

The multiplets obtained in this computation are then

nomistati = Table[""^ToString[2 autovS[[k]]+1]lStati[[autovL[[k]]+1]],{k,1,Length[autovS]}]

Matrix elements of H for equivalent electrons

The general form of matrix elements for the interaction between two Slater determinants can be reduced to matrix elements of the two particle operators (written here g ). The matrix elements have the form (direct term and exchange term, see text)

⟨(Li,mi)(Lj,mj)|g|(Li,mr)(Lj,mt)⟩ = (Subscript[δ, si,sr])(Subscript[δ, sj,st]) (∑_k)c[k,Li,mi,Li,mr]c[k,Lj,mj,Lj,mt](F^k)[Li,Lj];
⟨(Li,mi)(Lj,mj)|g|(Lj,mt)(Li,mr)⟩ =  (Subscript[δ, si,st])(Subscript[δ, sj,sr]) (∑_k)c[k,Li,mi,Lj,mt]c[k,Lj,mj,Lr,mr](G^k)[Li,Lj];

For equivalent electrons all angular momenta L are equal and the G functions coincide with the F functions (see text).

For diagonal matrix elements, i.e. between the same Slater determinant, mi=mr e mj = mt  and previous formulas simplify to

⟨(Li,mi)(Lj,mj)|g|(Li,mr)(Lj,mt)⟩ =  (∑_k)c[k,Li,mi,Li,mi]c[k,Lj,mj,Lj,mi](F^k)[Li,Lj];
⟨(Li,mi)(Lj,mj)|g|(Lj,mt)(Li,mr)⟩ =   (∑_k)(((c[k,Li,mi,Lj,mt]))^2)(F^k)[Li,Lj];

For diagonal matrix elements one has to sum over all possible couples of electrons.

Non diagonal matrix elements are different from zero only if all other states are unchanged, i.e. one has to select between the couples of Slater determinants those which differ by only two states.

The functions

The function c[k, l1, m1, l2, m2] compute the necessary coefficients, see notebook [NB2]. The function is wrapped with Check and Quiet to give 0 without producing error messages if k, for example, is not compatible with angular momenta addition rules.

The function calD[va,base] compute the diagonal matrix element  ⟨base[[va]] | H | base[[va]]⟩ by making a sum on couples of orbitals.

The function intEne[L] gives the mean energy (see text) for orbitals with angular momentum L, (the mean energy of a couple).

c[k_,l1_,m1_,l2_,m2_]:= 
Quiet[Check[(((-1))^m2)(Sqrt[(((2l1+1)(2l2+1))/(((2k+1))^2))])*
ClebschGordan[{l1,-m1},{l2,m2},{k,m2-m1}]*
ClebschGordan[{l1,0},{l2,0},{k,0}],0]]

Diagonal matrix elements

We define calcD

calcD[va_,base_]:=
(k1 = Flatten[Partition[base[[va]][[1]],2]];
mL = k1[[Range[1,Length[k1]-1,2]]];
mS = k1[[Range[2,Length[k1],2]]];
fK=0; fJ=0;
Do[
(fK = fK+Sum[c[k,Lorb,mL[[i]],Lorb,mL[[i]]]*c[k,Lorb,mL[[j]],Lorb,mL[[j]]]*
(F^ToString[k])[str1],{k,0,2Lorb,2}];
If[(mS[[i]]== mS[[j]]),fJ = fJ + Sum[c[k,Lorb,mL[[i]],Lorb,mL[[j]]]*c[k,Lorb,mL[[i]],Lorb,mL[[j]]]*
(F^ToString[k])[str1],{k,0,2*Lorb,2}];,fJ=fJ+0;];),{i,1,numE-1},{j,i+1,numE}];
Do[
(fK = fK+Sum[c[k,Lorb,mL[[i]],Lorb,mL[[i]]]*c[k,LorbNE,Last[mL],LorbNE,Last[mL]]*
(F^ToString[k])[str2],{k,0,Max[Lmax,2Lorb],2}];
If[mS[[i]]== Last[mS],fJ = fJ + Sum[c[k,Lorb,mL[[i]],LorbNE,Last[mL]]*c[k,Lorb,mL[[i]],LorbNE,Last[mL]]*
(G^ToString[k])[str2],{k,Lmin,Lmax,2}];,fJ=fJ+0;];),{i,1,numE}];
fK-fJ)

Off diagonal matrix elements

To these matrix elements contribute only states which differ by 2 orbitals.  Let inda  and   indb the two vectors with occupation numbers for orbitals a e b. The product inda*indb works as a logic AND so its trace gives the number of states common to the two orbitals, the variable diff gives the orbitals different and must be equal to 2 to have a non vanishing matrix element.

In the loop below orba,orbb are the filled orbitals not in common between the two states.

The vectors psia e psib are constructed with the two changing orbitals at the first position and computing the sign in the permutation.

The notation used for matrix elements is hijrt and hijtr corresponding to  ⟨ i j | H | r t ⟩  and  ⟨ i j | H | t r ⟩ in the matrix elements on the two states (i.e. the direct element and the exchange integral).

ham = Table[0,{i,1,Length[base]},{j,1,Length[base]}];
Do[
inda = indice[[va]];indb=indice[[vb]];
lasta = Take[base[[va]][[1]],-2];
lastb = Take[base[[vb]][[1]],-2];
diff = numE - Total[inda*indb];
Which[(diff==2) &&  lasta== lastb,
(orba = inda-inda*indb;orbb = indb-inda*indb;
numa =Select[Range[1,Length[numq]]orba,#1>0&];
numb =Select[Range[1,Length[numq]]orbb,#1>0&];
numCom = Select[Range[1,Length[numq]]inda*indb,#1>0&];
psia = Join[numq[[numa]],numq[[numCom]]];
psib = Join[numq[[numb]],numq[[numCom]]];
segno = Signature[psia]*Signature[psib];
si =psia[[1,2]];mi =psia[[1,1]];
sj = psia[[2,2]];mj =psia[[2,1]];
sr = psib[[1,2]];mr =psib[[1,1]];
st = psib[[2,2]];mt =psib[[2,1]];
hijrt=0;hijtr=0;
If[(mi+mj == mr + mt)&& (si== sr) && (sj== st),
hijrt = Sum[c[k,Lorb,mi,Lorb,mr]*c[k,Lorb,mt,Lorb,mj]*(F^ToString[k])[str1],{k,0,2Lorb,2}];,hijrt=0;];
If[(mi+mj == mr + mt) && (si== st) && (sj== sr),
hijtr = Sum[c[k,Lorb,mi,Lorb,mt]*c[k,Lorb,mr,Lorb,mj]*(F^ToString[k])[str1],{k,0,2Lorb,2}];,hijtr=0;];ham[[va,vb]]=segno*(hijrt -hijtr)),
(diff==1) &&  lasta≠  lastb,
(orba = inda-inda*indb;orbb = indb-inda*indb;
numa =Select[Range[1,Length[numq]]orba,#1>0&];
numb =Select[Range[1,Length[numq]]orbb,#1>0&];
numCom = Select[Range[1,Length[numq]]inda*indb,#1>0&];
psia = Join[numq[[numa]],numq[[numCom]]];
psib = Join[numq[[numb]],numq[[numCom]]];
segno = Signature[psia]*Signature[psib];
si =psia[[1,2]];mi =psia[[1,1]];
sj =lasta[[2]];mj =lasta[[1]];
sr = psib[[1,2]];mr =psib[[1,1]];
st = lastb[[2]];mt =lastb[[1]];
If[(mi+mj == mr + mt) && (si== sr) && (sj== st),
hijrt = Sum[c[k,Lorb,mi,Lorb,mr]*c[k,LorbNE,mt,LorbNE,mj]*(F^ToString[k])[str2],{k,0,Max[Lmax,2Lorb],2}];,hijrt=0;];
If[(mi+mj == mr + mt) && (si== st) && (sj== sr),
hijtr = Sum[c[k,Lorb,mi,LorbNE,mt]*c[k,Lorb,mr,LorbNE,mj]*(G^ToString[k])[str2],{k,Lmin,Lmax,2}];,hijtr=0;];
ham[[va,vb]]=segno*(hijrt - hijtr)),
1== 1,ham[[va,vb]]= 0;],
{va,1,Length[base]-1},{vb,va+1,Length[base]}];
ham = ham + Transpose[ham];

Full hamiltonian

In the previous loop we have computed the matrix ham off diagonal,  let us now add the diagonal terms:

Table[ham[[va,va]]=calcD[va,base],{va,1,Length[base]}];

Mean energy

The mean energy  (for each couple) is given by (see notebook [NB-15.2-Mean energy.nb] )

intEne[mom_]:= Expand[(F^ToString[0])[str1] - (1/(4mom+1))Sum[c[k,mom,0,mom,0]*(F^ToString[k])[str1],{k,2,2*mom,2}]]

intEneNE[mom1_,mom2_]:= Expand[(F^ToString[0])[str2] - (1/(2(Sqrt[(2mom1+1)(2mom2+1)])))Sum[c[k,mom1,0,mom2,0]*(G^ToString[k])[str2],{k,Lmin,Lmax,2}]]

eneav = Expand[Binomial[numE,2]*intEne[Lorb]+numE*intEneNE[Lorb,LorbNE]]

Diagonalization of H

In absence of additional degeneracies H must be diagonal in the base in which L2 e S2 are diagonal, i.e. eigvLS. For use in other program we print H - eneav. If everything is correct the matrix ⟨eigvLS[[i]] |ham | eigvLS[[j] ⟩ must be diagonal.

ham0 = ham;
ham = Simplify[ham - eneav*IdentityMatrix[Length[base]]];

hamdiag = Expand[Table[eigvLS[[i]].Simplify[ham.eigvLS[[j]]],{i,1,Length[base]},{j,1,Length[base]}]];

nomistati = Table[""^ToString[2 autovS[[k]]+1]lStati[[autovL[[k]]+1]],{k,1,Length[autovS]}]

We print the first maxPrint rows and columns of the matrix (in general can be quite large)

maxPrint = 3; 
numPrint = Min[maxPrint,Length[hamdiag[[1]]]];
indPrint = Range[numPrint];
nomistatiC=Join[{""},nomistati[[indPrint]]];
pip1= Join[{nomistati[[indPrint]]},hamdiag[[indPrint,indPrint]],1];
pip2 = Join[Transpose[{nomistatiC}],pip1,2];
MatrixForm[pip2]

This is the mean energy

Print["Emed = ",eneav]

We can have a graphical representation of the matrix, which can help for a large number of multiplets. In the graphical representation one can see that for (d^3)configurations, for instance, some 2×2 blocks are non diagonal, if need this block must be treated separately. This situation happens when appear several multiplets with the same (L,S) quantum numbers.

MatrixPlot[hamdiag]

Textual representation

The function RunLengthEncode[ L ] lists the equal contiguous components and the multiplicity of the list L, as an example

RunLengthEncode[{1,1,2,2,2,7,7}]  ->  {{1,2},{2,3},{7,2}}

This function is used to print the various non zero Hamiltonian terms, even in case of a non completely diagonal H.

RunLengthEncode[x_List]:={First[#],Length[#]}&/@Split[x]

runnomistati=RunLengthEncode[nomistati];
i1=1;j1=1;
Do[
j1 = i1+runnomistati[[k]][[2]] -1;
Print[runnomistati[[k]][[1]],":   ",MatrixForm[Simplify[hamdiag[[Range[i1,j1],Range[i1,j1]]]]]];
Print[];
i1 = i1+runnomistati[[k]][[2]];,
{k,1,Length[runnomistati]}]

The mean energy is given by

Print["Emed = ",eneav]

Some checks on H

Check 1

If our computation is correct the weighted mean of the subtracted hamiltonian (the matrix hamdiag) must be zero. The instruction Tr[A,List] extract the diagonal from a  matrix A:

Expand[Tr[(2autovL+1)*(2autovS+1)*Tr[hamdiag,List]]/Tr[(2autovL+1)*(2autovS+1)]]