Hartree-Fock 1

Quantum Mechanics
                      K. Konishi, G. Paffuti

Compute the ground-state energy of the Carbon in the Hartree-Fock approximation. Give an estimate of the energies of the (L,S) terms and compare with the results of HF for single levels.

NOTE: In this program there will be several global variables which must not be changed inadvertently by the user, for this reason their name is protected. The complete list of these variables is contained as argument in the functions Unprotect below.
Different names follow usual Mathematica rules.

Global variables and initialization

Off[General::spell, General::spell1];

Unprotect[cyclesHF,errHF,hgrid,x0grid,psiInput,alphaHF,inizioHF,firstEHF,printHF,hPerHartreeHF,gridOldHF,si$EneHF];
Unprotect[grid,electronsnames,FkHF,GkHF,epsLSlevel,virialeLSlevel,listadf,listadg,namesEl,namesMult,numQ$HF];

We change the default output of Mathematica function ElementData (see help) for Copper as it appears that the electron configuration we use us more stable in Hartree approximation.

ClearAll["Global`*"];
Off[General::spell, General::spell1,InterpolatingFunction::"dmval"];
Off[CompiledFunction::cfta];
Unprotect[ElementData];
ElementData/:ElementData["Cu","ElectronConfiguration"]={{2},{2,6},{2,6,9},{2}};
 ElementData/:ElementData[29,"ElectronConfiguration"]={{2},{2,6},{2,6,9},{2}}; 
ElementData/:ElementData["Cu","ElectronConfigurationString"]="[Ar]3(d^("9"))"<>"4(s^("2"))";
 ElementData/:ElementData[29,"ElectronConfigurationString"]= "[Ar]3(d^("9"))"<>"4(s^("2"))";
Protect[ElementData];

Protect[cyclesHF,errHF,hgrid,x0grid,psiInput,alphaHF,inizioHF,firstEHF,printHF,hPerHartreeHF,gridOldHF,si$EneHF];
Protect[grid,electronsnames,FkHF,GkHF,epsLSlevel,virialeLSlevel,listadf,listadg,namesEl,namesMult,numQ$HF];

Use of the notebook

What can be done and how

This notebook can be used  for the computation of ground state energy and energy of the few first excited levels of an atom in the Hartree Fock approximation. The user can compute energies and wave functions. The fine structure (i.e. spin interaction) is neglected.
To simplify as much as possible the use of the notebook, even for the novice in Mathematica, the basic usage needs a very limited number of functions (essentially two). Almost everything has a default value. The user can modify the behavior of the functions using options, in the usual Mathematica syntax var -> opz. The elementary approach used to realize this scheme require the use of some global variables, usually hidden to the user, but in any case accessible. These variables are used to allow the exchange of data between different functions. These variables are those listed at the beginning of the notebook. Their value is varied within some functions but not in the standard use of the program by the reader. The inexperienced reader can forget them.

Summary of the user interface

Basic input

The only input strictly necessary is the electronic configuration of the atom (and eventually its ionization degree).

The reader can choose between two options

degIoniz = 0;
electrons ={{"1s",2},{"2s",2},{"2p",4}};

electrons=ElementData[8,"ElectronConfiguration"];
electrons=ElementData["O","ElectronConfiguration"];

Both formats are accepted by the functions below. For excited states one must use the first option, as these are not contained in Mathematica database.

The next step is the determination of quantum numbers of the electrons. This is done with the function getquantumNum as shown below

quantumNum = getquantumNum[electrons,degIoniz];

The function choiceMultiplet[] gives a list of possible (L,S) levels, and the reader has to choose between them

listLS=choiceMultiplet[];

There are 3 possible choices (we call here LSlevel the variable which record the choice )

As an example for Oxygen listLS will have the form (({)^3)P((,)^1)S((,)^1)D}and the choice listLS[[2]] will correspond to the state 1S.

Limits on the computation

The default assignment is LSlevel = {ground}, than this instruction can be omitted in this case.

Functions for basic use

{ene,psi,enetot}= cycleHartreeFock[LSlevel];

compute the energy of LSlevel ( or of the default value, the ground state). There are some limitations on the computable terms (see above).

All energies are given in atomic units, except those printed during the cycle which are in Rydberg. We remember that 1 au = 2 Ry.

In the output psi[[k]] is a vector containing the values of the reduced radial wave function for the k-orbital on the grid of points used in the numerical computation.

A similar function is calcCore which compute the same quantities only for the closed shells.

The function

{eneMult,listaMult}=approxMultiplets[];

compute and print energies differences between computable energies. The omissions are due to the need of other quantum numbers to specify some levels, as in some excited states of the  (d^3) electron configuration.

plotPsi[psi,zoom];

This function gives a list of the graphs of different wave functions.

The reader can use  Show[] to see all functions at once, or  Show/@ to visualize each plot separately. The function has an optional list as argument, this can be used to select some specific orbitals. As an example Show[psi,zoom,{1,3}] visualize the orbitals 1 and 3.

The zoom parameter allow to enlarge or reduce the graph.

Options

The function cycleHartree has several options which can be used with the usual Mathematica mechanism  option → val.
The argument LSlevel is optional and both syntax are valid:

cycleHartreeFock[LSlevel, opz1-> val1, opz2-> val2];  cycleHartreeFock[ opz1-> val1, opz2-> val2];

The options and their default values are

Options[cycleHartreeFock]

Their meaning is the following

For calcCore there is a similar set of options

Options[calcCore]

si$EneHF: says if the total energy must be computed or not.

Functions

List for dfk and dgk

(* {{numE,La,Lb,L,S},{dfk, k=0,6}}  *)
Unprotect[listadf];

listadf={
{{2,1,1,1,1},{0,0,-(3/25),0,0,0,0}},
{{2,1,1,0,0},{0,0,(12/25),0,0,0,0}},
{{2,1,1,2,0},{0,0,(3/25),0,0,0,0}},
{{3,1,1,0,(3/2)},{0,0,-(9/25),0,0,0,0}},
{{3,1,1,1,(1/2)},{0,0,(6/25),0,0,0,0}},
{{3,1,1,2,(1/2)},{0,0,0,0,0,0,0}},
{{2,2,2,1,1},{0,0,(11/63),0,-(10/63),0,0}},
{{2,2,2,3,1},{0,0,-(58/441),0,(5/441),0,0}},
{{2,2,2,0,0},{0,0,(20/63),0,(20/63),0,0}},
{{2,2,2,2,0},{0,0,-(13/441),0,(50/441),0,0}},
{{2,2,2,4,0},{0,0,(50/441),0,(5/147),0,0}},
{{3,2,2,1,(3/2)},{0,0,(2/21),0,-(5/21),0,0}},
{{3,2,2,3,(3/2)},{0,0,-(31/147),0,-(10/147),0,0}},
{{3,2,2,1,(1/2)},{0,0,-(4/147),0,(10/147),0,0}},
{{3,2,2,2,(1/2)},{0,0,(5/21),0,(5/21),0,0}},
{{3,2,2,3,(1/2)},{0,0,(41/147),0,-(5/49),0,0}},
{{3,2,2,4,(1/2)},{0,0,-(19/147),0,(55/441),0,0}},
{{4,2,2,2,2},{0,0,-(5/21),0,-(5/21),0,0}},
{{4,2,2,1,1},{0,0,(1/21),0,(10/63),0,0}},
{{4,2,2,3,1},{0,0,(4/147),0,-(10/147),0,0}},
{{4,2,2,3,1},{0,0,(22/147),0,(5/49),0,0}},
{{4,2,2,0,0},{0,0,(380/1617),0,(940/1617),0,0}},
{{4,2,2,2,0},{0,0,(6221/20139),0,(3155/20139),0,0}},
{{5,2,2,0,(5/2)},{0,0,-(25/63),0,-(25/63),0,0}},
{{5,2,2,1,(3/2)},{0,0,-(16/63),0,(5/63),0,0}},
{{5,2,2,2,(3/2)},{0,0,-(22/441),0,-(85/441),0,0}},
{{5,2,2,3,(3/2)},{0,0,(23/441),0,-(40/441),0,0}},
{{5,2,2,4,(3/2)},{0,0,-(85/441),0,-(50/441),0,0}},
{{5,2,2,0,(1/2)},{0,0,(113/441),0,-(55/441),0,0}},
{{5,2,2,1,(1/2)},{0,0,(320/441),0,-(100/441),0,0}},
{{5,2,2,3,(1/2)},{0,0,-(85/441),0,(125/441),0,0}},
{{5,2,2,3,(1/2)},{0,0,(59/441),0,-(25/441),0,0}},
{{2,3,3,1,1},{0,0,(43/195),0,(19/429),0,-(875/5577)}},
{{2,3,3,3,1},{0,0,-(14/585),0,-(7/429),0,-(350/16731)}},
{{2,3,3,5,1},{0,0,-(53/585),0,-(155/4719),0,(2975/184041)}},
{{2,3,3,0,0},{0,0,(56/195),0,(28/143),0,(1400/5577)}},
{{2,3,3,2,0},{0,0,(307/2925),0,-(1/13),0,(175/1521)}},
{{2,3,3,4,0},{0,0,-(22/195),0,(1459/14157),0,(1750/61347)}},
{{2,3,3,6,0},{0,0,(77/585),0,(35/1573),0,(3325/184041)}},
{{3,3,3,0,(3/2)},{0,0,-(14/195),0,-(7/143),0,-(350/5577)}},
{{3,3,3,2,(3/2)},{0,0,(101/585),0,(5/429),0,-(3325/16731)}},{{3,3,3,3,(3/2)},{0,0,-(14/195),0,-(7/143),0,-(350/5577)}},
{{3,3,3,4,(3/2)},{0,0,(2/117),0,-(127/4719),0,-(20650/184041)}},
{{3,3,3,6,(3/2)},{0,0,-(133/585),0,-(413/4719),0,(4375/184041)}},
{{3,3,3,1,(1/2)},{0,0,-(29/585),0,(2/1287),0,(1225/16731)}},{{3,3,3,4,(1/2)},{0,0,(758/6825),0,(1301/9009),0,-(200/1859)}},
{{3,3,3,4,(1/2)},{0,0,(304/4095),0,-(512/33033),0,(8600/184041)}},{{3,3,3,5,(1/2)},{0,0,-(1189/8775),0,-(37/1287),0,(6475/50193)}}
};

Protect[listadf];

(* {{La,Lb,L,S},{dfk, k=0,6},{dgk, k=0,6}}  *)
Unprotect[listadg];

listadfdg = {
{{0,0,0,1},{0,0,0,0,0,0,0},{-(1/2),0,0,0,0,0,0}},{{0,0,0,0},{0,0,0,0,0,0,0},{(3/2),0,0,0,0,0,0}},
{{0,1,1,1},{0,0,0,0,0,0,0},{0,-(1/6),0,0,0,0,0}},{{0,1,1,0},{0,0,0,0,0,0,0},{0,(1/2),0,0,0,0,0}},
{{0,2,2,1},{0,0,0,0,0,0,0},{0,0,-(1/10),0,0,0,0}},{{0,2,2,0},{0,0,0,0,0,0,0},{0,0,(3/10),0,0,0,0}},
{{0,3,3,1},{0,0,0,0,0,0,0},{0,0,0,-(1/14),0,0,0}},{{0,3,3,0},{0,0,0,0,0,0,0},{0,0,0,(3/14),0,0,0}},
{{1,1,0,1},{0,0,(2/5),0,0,0,0},{-(5/6),0,-(1/3),0,0,0,0}},{{1,1,1,1},{0,0,-(1/5),0,0,0,0},{(7/6),0,-(2/15),0,0,0,0}},{{1,1,2,1},{0,0,(1/25),0,0,0,0},{-(5/6),0,(2/75),0,0,0,0}},{{1,1,0,0},{0,0,(2/5),0,0,0,0},{(7/6),0,(7/15),0,0,0,0}},{{1,1,1,0},{0,0,-(1/5),0,0,0,0},{-(5/6),0,(4/15),0,0,0,0}},{{1,1,2,0},{0,0,(1/25),0,0,0,0},{(7/6),0,(8/75),0,0,0,0}},
{{1,2,1,1},{0,0,(1/5),0,0,0,0},{0,0,0,-(3/14),0,0,0}},{{1,2,2,1},{0,0,-(1/5),0,0,0,0},{0,(4/15),0,-(3/70),0,0,0}},{{1,2,3,1},{0,0,(2/35),0,0,0,0},{0,-(1/3),0,(3/98),0,0,0}},{{1,2,1,0},{0,0,(1/5),0,0,0,0},{0,(2/15),0,(3/10),0,0,0}},{{1,2,2,0},{0,0,-(1/5),0,0,0,0},{0,-(2/15),0,(9/70),0,0,0}},{{1,2,3,0},{0,0,(2/35),0,0,0,0},{0,(7/15),0,(27/490),0,0,0}},
{{1,3,2,1},{0,0,(4/25),0,0,0,0},{0,0,(9/350),0,-(10/63),0,0}},{{1,3,3,1},{0,0,-(1/5),0,0,0,0},{0,0,(9/70),0,-(1/63),0,0}},{{1,3,4,1},{0,0,(1/15),0,0,0,0},{0,0,-(3/14),0,(5/189),0,0}},{{1,3,2,0},{0,0,(4/25),0,0,0,0},{0,0,(3/50),0,(2/9),0,0}},{{1,3,3,0},{0,0,-(1/5),0,0,0,0},{0,0,-(3/70),0,(5/63),0,0}},{{1,3,4,0},{0,0,(1/15),0,0,0,0},{0,0,(3/10),0,(1/27),0,0}},
{{2,2,0,1},{0,0,(2/7),0,(2/7),0,0},{-(9/10),0,-(9/35),0,-(9/35),0,0}},{{2,2,1,1},{0,0,(1/7),0,-(4/21),0,0},{(11/10),0,(6/35),0,-(17/105),0,0}},{{2,2,2,1},{0,0,-(3/49),0,(4/49),0,0},{-(9/10),0,(22/245),0,-(13/245),0,0}},{{2,2,3,1},{0,0,-(8/49),0,-(1/49),0,0},{(11/10),0,-(33/245),0,(2/245),0,0}},{{2,2,4,1},{0,0,(4/49),0,(1/441),0,0},{-(9/10),0,-(13/245),0,(58/2205),0,0}},{{2,2,0,0},{0,0,(2/7),0,(2/7),0,0},{(11/10),0,(11/35),0,(11/35),0,0}},{{2,2,1,0},{0,0,(1/7),0,-(4/21),0,0},{-(9/10),0,-(4/35),0,(23/105),0,0}},{{2,2,2,0},{0,0,-(3/49),0,(4/49),0,0},{(11/10),0,-(8/245),0,(27/245),0,0}},{{2,2,3,0},{0,0,-(8/49),0,-(1/49),0,0},{-(9/10),0,(47/245),0,(12/245),0,0}},{{2,2,4,0},{0,0,(4/49),0,(1/441),0,0},{(11/10),0,(27/245),0,(68/2205),0,0}},
{{2,3,1,1},{0,0,(8/35),0,(2/21),0,0},{0,(1/70),0,-(2/35),0,-(15/77),0}},{{2,3,2,1},{0,0,(2/35),0,-(1/7),0,0},{0,(9/70),0,(16/105),0,-(20/231),0}},{{2,3,3,1},{0,0,-(11/105),0,(2/21),0,0},{0,-(9/70),0,-(13/315),0,-(10/693),0}},{{2,3,4,1},{0,0,-(1/7),0,-(2/63),0,0},{0,(23/70),0,-(29/315),0,(10/693),0}},{{2,3,5,1},{0,0,(2/21),0,(1/231),0,0},{0,-(27/70),0,-(4/315),0,(160/7623),0}},{{2,3,1,0},{0,0,(8/35),0,(2/21),0,0},{0,(1/14),0,(2/21),0,(5/21),0}},{{2,3,2,0},{0,0,(2/35),0,-(1/7),0,0},{0,-(3/70),0,-(4/35),0,(10/77),0}},{{2,3,3,0},{0,0,-(11/105),0,(2/21),0,0},{0,(3/14),0,(5/63),0,(40/693),0}},{{2,3,4,0},{0,0,-(1/7),0,-(2/63),0,0},{0,-(17/70),0,(41/315),0,(20/693),0}},{{2,3,5,0},{0,0,(2/21),0,(1/231),0,0},{0,(33/70),0,(16/315),0,(170/7623),0}}
};

Protect[listadg];

Unprotect[namesEl];

namesEl={"H","He","Li","Be","B","C","N","O","F","Ne","Na","Mg","Al","Si","P","S","Cl","Ar","K","Ca","Sc","Ti","V","Cr","Mn","Fe","Co","Ni","Cu","Zn","Ga","Ge","As","Se","Br","Kr","Rb","Sr","Y","Zr","Nb","Mo","Te","Ru","Rh","Pd","Ag","Cd","In","Sn","Sb",
"Te","I","Xe","Cs","Ba","La","Ce","Pr","Nd","Pm","Sm","Eu","Gd","Tb","Dy","Ho","Er","Tm","Yb","Lu","Hf","Ta","W","Re","Os","Ir","Pt","Au","Hg","Tl","Pb","Bi","Po","At","Rn"};

Protect[namesEl];

Unprotect[namesMult];

namesMult = {"S","P","D","F","G","H","I","L","M","N"};

Protect[namesMult];

Electron configuration and multiplets

getquantumNum[electrons_,degIoniz_]:=
Module[{quantumNum,nomi,numPrinc,LperOrb},
nomi = {"s","p","d","f","g","h","i","l","m","n"};
If[NumberQ[electrons[[1,1]]],
quantumNum=getquantumNumB[electrons,degIoniz];
{LperOrb,numPrinc}= quantumNum[[{4,5}]];
(* global variable electronsnames *)
Unprotect[electronsnames];
electronsnames = Inner[StringJoin,ToString[#]&/@numPrinc,nomi[[LperOrb +1]],List];
Protect[electronsnames];
,
quantumNum=getquantumNumA[electrons,degIoniz];
Unprotect[electronsnames];
electronsnames = electrons[[All,1]];
Protect[electronsnames];
];
Unprotect[numQ$HF];
numQ$HF = quantumNum;
Protect[numQ$HF];
Return[quantumNum]];

getquantumNumA[electrons_,degIoniz_]:= Module[
{Z,quantumNum,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
Lstring},norbitals=Length[electrons];
numEperOrb = electrons[[Range[norbitals],2]];
Lstring =StringTake[# ,-1]&/@electrons[[Range[norbitals],1]];
LperOrb=Which[#== "s",0,#=="p",1,#== "d",2,#== "f",3]&/@Lstring;
numPrinc= ToExpression[StringTake[# ,1]&/@electrons[[Range[norbitals],1]]];
numrad = numPrinc - LperOrb;
Z = Total[numEperOrb] + degIoniz;
sigma=(Z-Table[ Total[numEperOrb[[Range[1,i]]]] -( numEperOrb[[i]]+1)/2,
{i,1,norbitals}])/Z;
gusciNC = Flatten[Position[Positive[(2*(2*LperOrb+1) - numEperOrb)],True]];
quantumNum = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC};
quantumNum]

getquantumNumB[electrons_,degIoniz_]:= Module[
{Z,quantumNum,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
ind},
norbitals=Length[Flatten[electrons]];
numEperOrb = Flatten[electrons];
LperOrb = Flatten[Range[0,Length[#]-1]&/@electrons];
numPrinc= Table[Table[i,{Length[electrons[[i]]]}],{i,1,Length[electrons]}];
numPrinc = Flatten[numPrinc];
numrad = numPrinc - LperOrb;
Z = Total[numEperOrb] + degIoniz;
sigma=(Z-Table[ Total[numEperOrb[[Range[1,i]]]] -( numEperOrb[[i]]+1)/2,
{i,1,norbitals}])/Z;
gusciNC = Flatten[Position[Positive[(2*(2*LperOrb+1) - numEperOrb)],True]];
quantumNum = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC};
quantumNum]

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

fmed[la_]:= Module[{fk},
fk= Table[0,{k,0,6}];
Table[fk[[k+1]]= -(2 la+1)/(4 la+1)*(ThreeJSymbol[{la,0},{k,0},{la,0}])^2;,
{k,2,2 la,2}];fk[[1]]=1;fk]

calcdFHund[Lorb_,numETrue_]:=
Module[{fK1,fJ1,mL1,df,a,numE,Lmult,Smult},
numE = Min[numETrue,2*(2*Lorb+1)-numETrue];
df= Table[0,{k,0,6}];
mL1= Table[Lorb - i,{i,0,numE-1}];
Do[
Do[
df[[k+1]]= df[[k+1]]+
c[k,Lorb,mL1[[i]],Lorb,mL1[[i]]]*c[k,Lorb,mL1[[j]],Lorb,mL1[[j]]]-
c[k,Lorb,mL1[[i]],Lorb,mL1[[j]]]^2;
,{i,1,numE-1},{j,i+1,numE}];
,{k,2,2Lorb,2}];
df[[1]]=1;
df = df - (numE*(numE-1))/2*fmed[Lorb];
fatt = Max[(numETrue*(numETrue-1)^2)/2,1];
fatt = (numETrue-1)/fatt;
df=df*fatt;
df[[1]]=0;
Smult = numE/2; 
Lmult = Sum[Lorb-k,{k,0,numE-1}];
{Lmult,Smult,df }]

calcdFdGHund[Lorb_,numETrue_]:=
Module[{df,numE,Lmult,Smult,dg},
{Lmult,Smult,df } = calcdFHund[Lorb,numETrue];
numE = Min[numETrue,2*(2*Lorb+1)-numETrue];
dg = -(numE/(2Lorb+1))(1/2)(1/numETrue);
Smult = Smult+1/2;
Return[{Lmult,Smult,df ,dg}];];

trovaMultipletto[numEQ_,la_,lb_,L_,S_,listadf_]:=Module[{lista,mult},
If[numEQ > 0,lista = {numEQ,la,lb,L,S};,lista = {la,lb,L,S};];
mult=Select[listadf,#[[1]] == lista&];
If[Length[mult]>0,Return[mult];,Print["LSlevel non found"];Abort[];];
];

fgMedi[nOrb_,LperOrb_]:=
Module[{la,lb,fk,gk},
gk=Table[0,{i1,1,nOrb},{i2,1,nOrb},{k,0,6}];
fk= Table[0,{a,1,nOrb},{k,0,6}];
fk[[All,1]] = 1;
Do[
la = LperOrb[[a]]; lb = LperOrb[[b]];
Table[gk[[a,b,k+1]]= -1/2*(ThreeJSymbol[{la,0},{k,0},{lb,0}])^2,{k,Abs[la-lb],la+lb,2}];,
{a,1,nOrb},{b,1,nOrb}];
Do[
la = LperOrb[[a]];
Table[fk[[a,k+1]]= -(((2 la+1))/((4 la+1)))*(ThreeJSymbol[{la,0},{k,0},{la,0}])^2;,
{k,2,2 la,2}];,
{a,1,nOrb}];
{fk,gk}
];

df1GuscioInc[numETrue_,orbEQ_,LeQ_,Ltot_,Stot_,listadf_]:=
Module[{nEeff,numCoppie,df,LMult,SMult,dfTot},
df= Table[0,{k,0,6}];
nEeff = Min[numETrue,2(2*LeQ+1)-numETrue];
numCoppie = numETrue*(numETrue-1)/2;
dfTot = trovaMultipletto[nEeff,LeQ,LeQ,Ltot,Stot,listadf];
If[Length[dfTot]>0,
LMult=Ltot;
SMult=Stot;
df = dfTot[[1,2]]/numCoppie;,
{LMult,SMult,df}=calcdFHund[LeQ,numETrue];
];
{LMult,SMult,df}];

dfdg2GusciInc[numEperOrb_,gusciNC_,LperOrb_,Linput_,Sinput_,listadfdg_]:=
Module[{orb1,orb2,Lorb1,Lorb2,dfTot,dfCoppia,dg,num1,num2},
orb1=gusciNC[[1]];orb2=gusciNC[[2]];
num1 = Min[numEperOrb[[orb1]],2*(2LperOrb[[orb1]]+1) - numEperOrb[[orb1]]];
num2 = Min[numEperOrb[[orb2]],2*(2LperOrb[[orb2]]+1) - numEperOrb[[orb2]]];
{Lorb1,Lorb2} = Sort[Flatten[{LperOrb[[orb1]],LperOrb[[orb2]]}]];
dfTot = trovaMultipletto[-2,Lorb1,Lorb2,Linput,Sinput,listadfdg];
dfCoppia = dfTot[[1,2]]/(num1*num2);
dg = dfTot[[1,3]]/(num1*num2);
{dfCoppia,dg}];

parametriMultipletto[quantumNum_,LSlevelIn_,listadf_,listadfdg_]:=
Module[{gk,fk,la,lb,df,dg,dfTot,gusciNC,orbEQ,momL,numEeff,fatt,orbL,orbS,
Linput,Sinput,LMult,SMult,dfCoppia,nOrb,Ltot,Stot,
orb1,orb2,Lorb1,Lorb2,LperOrb,numEperOrb,nomi,StatoLista,num1,num2,LSlevel},

{nOrb,numEperOrb,LperOrb,gusciNC} = quantumNum[[{2,3,4,8}]]; 
nomi = {"S","P","D","F","G","H","I","L","M","N"};
LMult = -1;SMult = -1;
LSlevel = LSlevelIn;

StatoLista = True;

{fk,gk}=fgMedi[nOrb,LperOrb];
dfCoppia = Table[0.,{k,0,6}];

If[Length[gusciNC]==2,
orb1 = gusciNC[[1]];orb2=gusciNC[[2]];
num1 = Min[numEperOrb[[orb1]],2*(2LperOrb[[orb1]]+1) - numEperOrb[[orb1]]];
num2 = Min[numEperOrb[[orb2]],2*(2LperOrb[[orb2]]+1) - numEperOrb[[orb2]]];
]


Which[
Length[LSlevel]== 0,
StatoLista = False;
,
Length[LSlevel] == 1,
Which[
Length[gusciNC]==0,
{LMult,SMult} = {0,0};
StatoLista = False;
,
Length[gusciNC]==1,
orbEQ = gusciNC[[1]];
{LMult,SMult,df}=calcdFHund[LperOrb[[orbEQ]],numEperOrb[[orbEQ]]];
fk[[orbEQ,All]] = fk[[orbEQ,All]]+df;
StatoLista = False;
,
Length[gusciNC]==2 && LperOrb[[gusciNC[[1]]]]*LperOrb[[gusciNC[[2]]]]== 0,
{orb1,orb2} = Sort[gusciNC,LperOrb[[#1]]≥ LperOrb[[#2]] &];
{LMult,SMult,df,dg}=calcdFdGHund[LperOrb[[orb1]],numEperOrb[[orb1]]];
fk[[orb1,All]] = fk[[orb1,All]]+df;
gk[[orb1,orb2,LperOrb[[orb1]]+1]] =  gk[[orb1,orb2,LperOrb[[orb1]]+1]]+dg;
gk[[orb2,orb1,LperOrb[[orb1]]+1]] =  gk[[orb2,orb1,LperOrb[[orb1]]+1]]+dg;
StatoLista = False;
,
Length[gusciNC]==2 && num1*num2== 1,
LSlevel = {LperOrb[[gusciNC[[1]]]]+LperOrb[[gusciNC[[2]]]],1};
,
True,
StatoLista = False;];
,
Length[LSlevel] == 2,
,
True,
StatoLista = False;];

If[StatoLista== False, Return[{LMult,SMult,fk,gk,dfCoppia}];];

(*  LSlevel has now length 2  *)
{Ltot,Stot} = LSlevel;
Which[
Length[gusciNC]== 1,
orbEQ = gusciNC[[1]];
{LMult,SMult,df}=df1GuscioInc[numEperOrb[[orbEQ]],orbEQ,LperOrb[[orbEQ]],Ltot,Stot,listadf];
fk[[orbEQ,All]] = fk[[orbEQ,All]]+df;
,
Length[gusciNC]== 2,
LMult = Ltot;SMult = Stot;
{dfCoppia,dg}=
dfdg2GusciInc[numEperOrb,gusciNC,LperOrb,Ltot,Stot,listadfdg];
orb1=gusciNC[[1]];orb2=gusciNC[[2]];
gk[[orb1,orb2,All]] = gk[[orb1,orb2,All]]+ dg;
gk[[orb2,orb1,All]] = gk[[orb1,orb2,All]] ;
];
Return[{LMult,SMult,fk,gk,dfCoppia}];]

choiceMult[]:=
Module[{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
Smult,Lmult,orb,lista,listanomi,listLS,nomi,mult,nomefond,quantumNum,numE},
quantumNum = numQ$HF;
{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC} =
quantumNum;
nomi = {"S","P","D","F","G","H","I","L","M","N"};
Print["Choose between:   mean E: {}","  ground state: {ground} "];
listLS = {};
listanomi = {};
nomefond = "";
Which[
Length[gusciNC] == 0,

nomefond = ""^ToString[1]nomi[[1]];
listanomi = ""^ToString[1]nomi[[1]];
listLS = {0,0};
,
Length[gusciNC] == 1,

orb = gusciNC[[1]];
numE = Min[numEperOrb[[orb]],2*(2*LperOrb[[orb]]+1)-numEperOrb[[orb]]];
Smult = numE/2; 
Lmult = Sum[LperOrb[[orb]]-k,{k,0,numE-1}];
nomefond = ""^ToString[2 Smult +1]nomi[[Lmult+1]];
lista = {numE,LperOrb[[orb]],LperOrb[[orb]]};
mult=Select[listadf,#[[1]][[Range[3]]] == lista&];
Do[
 Smult = mult[[j]][[1,5]];
Lmult = mult[[j]][[1,4]];
listanomi = Append[listanomi,""^ToString[2 Smult +1]nomi[[Lmult+1]]];
,{j,1,Length[mult]}];
listLS = Table[{mult[[j]][[1,4]],mult[[j]][[1,5]]},{j,1,Length[mult]}];
,
Length[gusciNC] == 2,

lista = Sort[LperOrb[[gusciNC]]];
mult=Select[listadfdg,#[[1]][[Range[2]]] == lista&];
Do[
 Smult = mult[[j]][[1,4]];
Lmult = mult[[j]][[1,3]];
listanomi = Append[listanomi,""^ToString[2 Smult +1]nomi[[Lmult+1]]];
,{j,1,Length[mult]}];
listLS = Table[{mult[[j]][[1,3]],mult[[j]][[1,4]]},{j,1,Length[mult]}];
];
Print["or {L,S} or listLS[[k]] from: ",listanomi," fond = ",nomefond];
Return[listLS]];

Grid and initializations

parametrigrid[quantumNum_,hin_,x0in_]:=Module[
{x,x0,h,Z,nP,xmax,rho,rho2,rho32,grid,numrad},
h = hin;x0=x0in;
{Z,numrad }= quantumNum[[{1,6}]];
(* h = 1./10; x0=0.10; *)
xmax = 40.*Max[numrad]^2 * Z;
nP = Ceiling[ Log[xmax/x0]/h];
nP = nP + Mod[nP,2];
x = x0*(Exp[h*Range[nP]] - 1.);
rho2 = x0^2Exp[2*h*Range[nP]];
rho = x+x0; 
rho32 = rho2/Sqrt[rho];
grid = {x,x0,h,rho,rho2};
grid];

getgridFixed[maxnumrad_,Z_,hin_,x0in_]:= Module[
{x,x0,h,nP,xmax,rho,rho2,rho32,grid},
h = hin;x0=x0in;
xmax = 40.*maxnumrad^2 * Z;
nP = Ceiling[ Log[xmax/x0]/h];
nP = nP + Mod[nP,2];
x = x0*(Exp[h*Range[nP]] - 1.);
rho2 = x0^2Exp[2*h*Range[nP]];
rho = x+x0; 
rho32 = rho2/Sqrt[rho];
grid = {x,x0,h,rho,rho2};
grid];

updategridPsi[x0_,h_,quantumNum_,gridIn_,psiIn_]:=
Module[{psi,x,rho,rho2,funz,xTot,xIn,x0In,hIn,rhoIn,rho2In,norbitals,xmax},
{xIn,x0In,hIn,rhoIn,rho2In} = gridIn;
norbitals = quantumNum[[2]];
Unprotect[grid];
grid = parametrigrid[quantumNum,h,x0];
Protect[grid];
x = grid[[1]];
xmax = Max[xIn,x]+0.01;
xTot = Flatten[{0.,xIn,xmax}];
psi = Table[0.,{i,1,norbitals}];
Do[
funz = Interpolation[Transpose[{xTot,Flatten[{0.,psiIn[[i]],0.}]}]];
psi[[i]] = funz[x];,
{i,1,norbitals}];
Return[psi];];

(* Integral on the exponential grid, ξ, for functions which vanishes at the origin *)
integ = Compile[{{f,_Real,1},{h,_Real}},
Module[{nL,f0,val},
nL = Length[f];
f0=0.;val=0.;
val = h*(119./90.*f[[1]] + 32./45.f[[2]] + 113./90.f[[3]] +
113./90. f[[nL-1]] + 17./45.f[[nL]]);
val = val + h/90.*68.*Total[f[[Range[4,nL-2,2]]]] +
 h/90.*112.*Total[f[[Range[5,nL-3,2]]] ];
val]];

InitMatrici[quantumNum_,grid_]:=
Module[{Kin,delta2,matS,vec8,x,h,
idiag,iPdiag,iMdiag,rho2},
(* grid = {x,x0,h,rho,rho2}; *)
{x,h,rho2} = grid[[{1,3,5}]];
idiag = Table[{i,i},{i,1,Length[x]}];
    iPdiag = Table[{i,i+1},{i,1,Length[x]-1}];
    iMdiag = Table[{i,i-1},{i,2,Length[x]}];delta2= N[SparseArray[{{i_,i_} ->  -2.,{i_,j_}/;j-i==1-> 1.,
{i_,j_}/; i-j== 1->  1.},{Length[x],Length[x]}]];
     Kin = -0.5/h^2*delta2;
vec8 = Table[1./8.,{i,1,Length[x]}];
matS = SparseArray[{idiag-> 10./12.*rho2,iPdiag->Drop[rho2/12.,1],iMdiag->Drop[rho2/12.,-1]},
{Length[x],Length[x]}];
{Kin,matS,vec8}];

χ[r_, n_, L_] := r* (2/(n^2))(Sqrt[(((n-L-1)!)/((n+L)!))])*
(((((2 r)/n)))^L) * (E^(-(r/n)))* LaguerreL[n-L-1, 2 L + 1,((2 r)/n)] ;

InitPsi[quantumNum_,grid_]:= Module[
{norma,psi,rho2,rho,x0,Z,norbitals,LperOrb,numPrinc,sigma,h,x},
(* nQ = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}; *)
{Z,norbitals,LperOrb,numPrinc,sigma}=quantumNum[[{1,2,4,5,7}]];
{x,x0,h,rho,rho2} = grid;
psi = Table[0.,{i,1,norbitals}];
Table[psi[[i]] =(1/(Sqrt[rho]))(1/(Sqrt[Z])) χ[sigma[[i]]*x, numPrinc[[i]], LperOrb[[i]]],{i,1,norbitals}];
Table[(norma = integ[rho2*Chop[psi[[i]]^2],h];
psi[[i]] =Chop[psi[[i]]/Sqrt[norma]];),{i,1,norbitals}];
psi]

Initalpha[quantumNum_]:= Module[{alpha,E0,Z,norbitals,LperOrb,numPrinc},
{Z,norbitals,LperOrb,numPrinc}=quantumNum[[{1,2,4,5}]];
alpha=Table[Which[
LperOrb[[j]] > 1 || numPrinc[[j]] > Max[numPrinc]-1,0.8,
 numPrinc[[j]] > 4,0.8,True,0.5],{j,1,norbitals}];
Return[alpha];];

InitEne[quantumNum_]:= Module[{alpha,E0,Z,norbitals,LperOrb,numPrinc},
(* nQ = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC} *)
{Z,norbitals,LperOrb,numPrinc}=quantumNum[[{1,2,4,5}]];
E0 = Table[1./(2*numPrinc[[i]]^2),{i,1,norbitals}];
alpha=Table[Which[
LperOrb[[j]] > 1 || numPrinc[[j]] > Max[numPrinc]-1,0.8,
 numPrinc[[j]] > 4,0.8,True,0.5],{j,1,norbitals}];
{E0,alpha}];

Numerical functions Zk, Yk and potentials

(* input functions are rescaled with ψ=P/(Sqrt[ρ]) *)
funZ=
Compile[{{valk,_Integer},{La,_Integer},{Lb,_Integer},
{fa,_Real,1},{fb,_Real,1},{x,_Real,1},{x0,_Real},{h,_Real}},
Module[{nL,fH,fun,cost1,cost2,fZ,j,k,kind,asint},
(nL = Length[fa];
If[valk== 0,fH = (x+x0)^2*fa*fb,fH = x^valk*(x+x0)^2*fa*fb];
If[La == Lb && fa[[3]] == fb[[3]] && fa[[5]] == fb[[5]] && valk==0,
asint =1.;,asint = 0.];
fZ=Table[0.,{j,1,nL}];
cost1= N[(1./(La+Lb+valk+3.))];
cost2=N[(1./(La+Lb+valk+4.))((1/(La+1.))+(1/(Lb+1.)))];
(* here r* Sqrt[rho] ψ_a Sqrt[rho] ψ_b = r fun /(r+r0) *)
fZ[[1]] = cost1*x[[1]]*fH[[1]]/(x[[1]]+x0)*(1+cost2*x[[1]]);
fZ[[2]] = cost1*x[[2]]*fH[[2]]/(x[[2]]+x0)*(1+cost2*x[[2]]);
fZ[[3]] = fZ[[1]]+ h/3.*(fH[[1]]+4.*fH[[2]]+fH[[3]]);
Table[fZ[[j]] = fZ[[j-2]]+ h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,5,nL-1,2}];
Table[fZ[[j]] = fZ[[j-2]]+ h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,4,nL-2,2}];
fZ[[nL]] = fZ[[nL-2]]+ h/3.*(fH[[nL-2]]+4.*fH[[nL-1]]+fH[[nL]]);
If[asint > 0.1,
fZ[[Range[1,nL-1,2]]] = fZ[[Range[1,nL-1,2]]] + asint - fZ[[nL-1]];
fZ[[Range[2,nL,2]]] = fZ[[Range[2,nL,2]]] + asint - fZ[[nL]]; ,];
If[valk== 0,fZ = fZ ;,fZ = fZ/x^valk];
fZ)]];

funY=
Compile[{{valk,_Integer},{La,_Integer},{Lb,_Integer},
{fa,_Real,1},{fb,_Real,1},{fZ,_Real,1},{x,_Real,1},{x0,_Real},{h,_Real}},
Module[{nL,fH,j,corr,k,fY},
nL = Length[fa];
fH = (2.*valk+1.)*(x+x0)*fZ/x^(valk+2);
fY=Table[0.,{j,1,nL}];
fY[[nL]] = fZ[[nL]]/x[[nL]]^(valk+1);
fY[[nL-1]]=  fZ[[nL-1]]/x[[nL-1]]^(valk+1);
Table[fY[[j-2]] = fY[[j]]+ 
h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,nL-1,5,-2}];
fY[[1]] = fY[[3]] + h/3.*(fH[[1]]+4.*fH[[2]]+fH[[3]]);
fY[[nL-2]] = fY[[nL]] +h/3.*(fH[[nL-2]]+4.*fH[[nL-1]]+fH[[nL]]);
Table[fY[[j-2]] = fY[[j]]+
 h/90.*(-fH[[j+1]] + 34.fH[[j]] +114.*fH[[j-1]] +
34.*fH[[j-2]] - fH[[j-3]]);,{j,nL-2,4,-2}];
fY = fY*x^(valk+1);
fY]];

Potentials Y[k,a,b] due to  psi[a] and psi[b].

calcpotHartree[psi_,LperOrb_,x_,x0_,h_]:=
Module[{Y0,minL,maxL,norbitals,fZ1,fY1,k,nP},
nP = Length[psi[[1]]];
norbitals = Length[LperOrb];
k=0;
Y0 = Table[0.,{i,1,norbitals}];
Do[
fZ1 = funZ[k,LperOrb[[i]],LperOrb[[i]],psi[[i]],psi[[i]],x,x0,h];
fY1=funY[k,LperOrb[[i]],LperOrb[[i]],psi[[i]],psi[[i]],fZ1,x,x0,h];
Y0[[i]]=fY1/x;
,
{i,1,norbitals}];
Return[Y0]]

calcpotYk[psi_,LperOrb_,x_,x0_,h_]:=
Module[{Yk,minL,maxL,nP,norbitals,fY1,fZ1},
nP = Length[psi[[1]]];
norbitals = Length[LperOrb];
Yk = Table[0.,{i,1,norbitals},{j,1,norbitals},{k,0,6}];
Do[
minL = Abs[LperOrb[[i]]- LperOrb[[j]]];
maxL = LperOrb[[i]]+ LperOrb[[j]];
Table[
fZ1 = funZ[k,LperOrb[[i]],LperOrb[[j]],psi[[i]],psi[[j]],x,x0,h];
fY1=funY[k,LperOrb[[i]],LperOrb[[j]],psi[[i]],psi[[j]],fZ1,x,x0,h];
Yk[[i,j,k+1]]=fY1/x;,{k,minL,maxL,2}];
,
{i,1,norbitals},{j,i,norbitals}];
Return[Yk]];

calcpotYorb[psi_,Yk_,fk_,gk_,dfCoppia_,quantumNum_]:= Module[
{Z,norbitals,numEperOrb,LperOrb,ib,ia,i2,
La,VYq,Yorb,indice2,gusciNC,nP},
{Z,norbitals,numEperOrb,LperOrb,gusciNC} = quantumNum[[{1,2,3,4,8}]];
indice2=Flatten[{Range[norbitals],Range[norbitals]}];
nP = Length[Yk[[1,1,1]]];
Yorb = Table[0.,{j,1,norbitals}];
Do[
La = LperOrb[[i]];
VYq = (numEperOrb[[i]]-1)*Sum[Yk[[i,i,k]]*fk[[i,k]],{k,1,2*La+1,2}];
Do[
ib = indice2[[j]];
VYq = VYq + numEperOrb[[ib]]*Yk[[ib,ib,1]]*fk[[ib,1]];,
{j,i+1,i+norbitals-1}];
Yorb[[i]]=VYq/Z;
,
{i,1,norbitals}];

If[Length[gusciNC]==2,
indice2={1,2,1,2};
Do[
ia = gusciNC[[i]];
La = LperOrb[[ia]];
VYq = 0.;
Do[
ib = gusciNC[[indice2[[j]]]];
VYq = VYq + numEperOrb[[ib]]*
Sum[Yk[[ib,ib,k]]*dfCoppia[[k]],{k,1,2*La+1,2}];
,
{j,i+1,i+Length[gusciNC]-1}];
Yorb[[ia]]=Yorb[[ia]]+VYq/Z;
,
{i,1,2}];
];

Return[Yorb]]

calcpotXorb[psi_,Yk_,fk_,gk_,quantumNum_]:= Module[
{Z,norbitals,numEperOrb,LperOrb,ia,ib,i2,minL,maxL,
La,potYJ,Xorb,indice2,nP},
{Z,norbitals,numEperOrb,LperOrb} = quantumNum[[{1,2,3,4}]];
indice2=Flatten[{Range[norbitals],Range[norbitals]}];
nP = Length[Yk[[1,1,1]]];
Xorb = Table[0.,{j,1,norbitals}];
Do[
La = LperOrb[[i]];
potYJ = Table[0.,{k,1,nP}];
Do[
i2 = indice2[[irun]];
ia = Min[i,i2];
ib = Max[i,i2];
minL = Abs[LperOrb[[i]] - LperOrb[[i2]]];
maxL = LperOrb[[i2]] + LperOrb[[i]];
potYJ = potYJ + numEperOrb[[i2]]*Sum[ gk[[i,i2,k+1]]*psi[[i2]]*
Yk[[ia,ib,k+1]],{k,minL,maxL,2}];,
{irun,i+1,i+norbitals-1}];
Xorb[[i]]=potYJ/Z;
,
{i,1,norbitals}];
Return[Xorb]]

Energy functions

calcEneHartree[psi_,grid_,quantumNum_]:=
Module[{Z,norbitals,numEperOrb,LperOrb,gusciNC,
x,x0,h,rho2,Kin,Yk,enekinsing,enekinsing1,indice2,Fk,Gk,ia,ib,
enepot,enetot,viriale,matS,vec8,delta2,enekin},

{Z,norbitals,numEperOrb,LperOrb,gusciNC} = quantumNum[[{1,2,3,4,8}]];
{x,x0,h,rho2} = grid[[{1,2,3,5}]];
{Kin,matS,vec8} = InitMatrici[quantumNum,grid];

Yk=calcpotYk[psi,LperOrb,x,x0,h];
delta2= N[SparseArray[{{i_,i_} ->  -2.,{i_,j_}/;j-i==1-> 1.,
{i_,j_}/; i-j== 1->  1.},{Length[x],Length[x]}]];

enekinsing = enekinsing1=Table[0.,{i,1,norbitals}];
enekinsing1 = Table[integ[rho2*psi[[i]]^2*(1./8/rho2+ LperOrb[[i]]*(LperOrb[[i]]+1)/(2 x^2)),h],{i,1,norbitals}];
enekinsing =  Table[integ[psi[[i]]*(Kin .(psi[[i]] - 1./12. delta2.psi[[i]])),h],{i,1,norbitals}];
enekinsing = enekinsing + enekinsing1;
enekin = Total[numEperOrb*enekinsing]*Z^2;

Fk = Table[0.,{ia,1,norbitals},{ib,1,norbitals}];
Do[
Fk[[ia,ib]] = 1./Z*integ[rho2*psi[[ia]]^2*Yk[[ib,ib,1]],h];,
{ia,1,norbitals},{ib,ia,norbitals}];

enepot = 0.;
Do[
enepot = enepot-numEperOrb[[ia]]*integ[rho2*psi[[ia]]^2/x,h] +
 numEperOrb[[ia]]*(numEperOrb[[ia]]-1)/2*
Fk[[ia,ia]],{ia,1,norbitals}];
Do[
enepot= enepot + numEperOrb[[ia]]*numEperOrb[[ib]]*Fk[[ia,ib]];
,
{ia,1,norbitals-1},{ib,ia+1,norbitals}];

enepot = enepot*Z^2;
enetot = enekin + enepot;
viriale = -enepot/enekin;
Unprotect[virialeLSlevel];
virialeLSlevel = viriale;
Protect[virialeLSlevel];
{enetot,enekin,viriale}];

calcEne[psi_,grid_,quantumNum_,fk_,gk_,dfCoppia_]:=
Module[{Z,norbitals,numEperOrb,LperOrb,gusciNC,
x,x0,h,rho2,Kin,Yk,enekinsing,enekinsing1,indice2,ia,ib,
enepot,enetot,viriale,Fk,Gk,minL,maxL,matS,vec8,delta2,enekin},

{Z,norbitals,numEperOrb,LperOrb,gusciNC} = quantumNum[[{1,2,3,4,8}]];
{x,x0,h,rho2} = grid[[{1,2,3,5}]];
{Kin,matS,vec8} = InitMatrici[quantumNum,grid];

Yk=calcpotYk[psi,LperOrb,x,x0,h];
delta2= N[SparseArray[{{i_,i_} ->  -2.,{i_,j_}/;j-i==1-> 1.,
{i_,j_}/; i-j== 1->  1.},{Length[x],Length[x]}]];

enekinsing = enekinsing1=Table[0.,{i,1,norbitals}];
enekinsing1 = Table[integ[rho2*psi[[i]]^2*(1./8/rho2+ LperOrb[[i]]*(LperOrb[[i]]+1)/(2 x^2)),h],{i,1,norbitals}];
enekinsing =  Table[integ[psi[[i]]*(Kin .(psi[[i]] - 1./12. delta2.psi[[i]])),h],{i,1,norbitals}];
enekinsing = enekinsing + enekinsing1;
enekin = Total[numEperOrb*enekinsing]*Z^2;

indice2=Flatten[{Range[norbitals],Range[norbitals]}];
Fk = Table[0.,{ia,1,norbitals},{ib,1,norbitals},{k,0,6}];
Gk = Table[0.,{ia,1,norbitals},{ib,1,norbitals},{k,0,6}];
(* Integrals Fk e Gk *)
Do[
minL  =Abs[ LperOrb[[ia]] - LperOrb[[ib]]];
maxL = LperOrb[[ia]] + LperOrb[[ib]];
Table[Fk[[ia,ib,k+1]] = 1./Z*integ[rho2*psi[[ia]]^2*Yk[[ib,ib,k+1]],h];,{k,0,6}];,
{ia,1,norbitals},{ib,ia,norbitals}];
Do[
minL  =Abs[ LperOrb[[ia]] - LperOrb[[ib]]];
maxL = LperOrb[[ia]] + LperOrb[[ib]];
Table[Gk[[ia,ib,k+1]] = 1./Z*integ[rho2*psi[[ia]]psi[[ib]]*Yk[[ia,ib,k+1]],h];,{k,minL,maxL,2}];,
{ia,1,norbitals},{ib,ia+1,norbitals}];
enepot = 0.;
Do[
enepot = enepot-numEperOrb[[ia]]*integ[rho2*psi[[ia]]^2/x,h] +
 numEperOrb[[ia]]*(numEperOrb[[ia]]-1)/2*
Sum[fk[[ia,k+1]]*Fk[[ia,ia,k+1]],{k,0,2*LperOrb[[ia]],2}];,
{ia,1,norbitals}];
Do[
enepot= enepot + numEperOrb[[ia]]*numEperOrb[[ib]]*Fk[[ia,ib,1]];
minL = Abs[ LperOrb[[ia]] - LperOrb[[ib]]];
maxL = LperOrb[[ia]] + LperOrb[[ib]];
enepot = enepot+ numEperOrb[[ia]]*numEperOrb[[ib]]*Sum[gk[[ia,ib,k+1]]*Gk[[ia,ib,k+1]],{k,minL,maxL,2}];,
{ia,1,norbitals-1},{ib,ia+1,norbitals}];
If[Length[gusciNC]== 2 && 
numEperOrb[[gusciNC[[1]]]] *numEperOrb[[gusciNC[[2]]]] ==1 ,
{ia,ib} = Sort[gusciNC];
enepot = enepot + dfCoppia.Fk[[ia,ib,All]];,];
enepot = enepot*Z^2;
enetot = enekin + enepot;
viriale = -enepot/enekin;
(* definition of gloabl variables *)
Unprotect[virialeLSlevel,FkHF,GkHF];
virialeLSlevel = viriale;
FkHF = (Z^2)*Fk; GkHF = (Z^2)*Gk;
Protect[virialeLSlevel,FkHF,GkHF];

{enetot,enekin,viriale}];

approxMultipletsOld[]:=
Module[{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
nEQ,orbEQ,LeQ,nEQeff,lista,ene,eneMult,mult,listaMult,Smult,Lmult,Fk,Gk,
quantumNum,au2Ry,Ry2cm},

Fk = FkHF; Gk = GkHF; 
quantumNum = numQ$HF;

{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC} = quantumNum;
If[Length[gusciNC]==1 ,
orbEQ = gusciNC[[1]];
nEQ = numEperOrb[[orbEQ]];
nEQeff = Min[nEQ,2(2LeQ+1)-nEQ];
 LeQ =LperOrb[[orbEQ]];
,nEQ=0;];
eneMult = {};
Which[
Length[gusciNC]==1 && Min[nEQ,2(2 LeQ+1)-nEQ]>1,
lista = {nEQeff,LeQ};
mult=Select[listadf,#[[1]][[{1,2}]] == lista&];
(* Warning  in listadf there is the total energy, not for a couple *)
ene=Table[mult[[j,2]].Fk[[orbEQ,orbEQ,All]],{j,1,Length[mult]}];
eneMult=Table[{(ene[[i]]-Min[ene]),mult[[i,1]][[4]],
mult[[i,1]][[5]]},{i,1,Length[mult]}];
eneMult= Sort[eneMult,#1[[1]] < #2[[1]] &];
];
listaMult = Table[
""^ToString[2 eneMult[[j,3]] +1]namesMult[[eneMult[[j,2]]+1]],{j,1,Length[eneMult]}];
au2Ry = 2;
Ry2cm = 109737.31568525;
If[Length[eneMult] == 0,
Print["No computable structure "];,
Print[Table[{listaMult[[i]],eneMult[[i,1]]*au2Ry*Ry2cm *"cm"^ToString[-1], "=" ,eneMult[[i,1]]"au" } ,{i,1,Length[eneMult]}]//TableForm];
];
Return[eneMult];];

approxMultiplets[]:=
(* Module[{Z,nOrbitali,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
nEQ,orbEQ,LeQ,nEQeff,lista,ene,eneMult,mult,listaMult,Smult,Lmult,Fk,Gk,
numQuantici,orb1,orb2,num2,num1,Lorb1,Lorb2}, *)
Module[{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
nEQ,orbEQ,LeQ,nEQeff,lista,ene,eneMult,mult,listaMult,Smult,Lmult,Fk,Gk,
quantumNum,orb1,orb2,num2,num1,Lorb1,Lorb2,au2Ry,Ry2cm},

Fk = FkHF; Gk = GkHF; 
quantumNum = numQ$HF;

{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC} = quantumNum;

If[Length[gusciNC]==1 ,
orbEQ = gusciNC[[1]];
nEQ = numEperOrb[[orbEQ]];
nEQeff = Min[nEQ,2(2LeQ+1)-nEQ];
 LeQ =LperOrb[[orbEQ]];
,nEQ=0;];

If[Length[gusciNC]==2 ,
orb1=gusciNC[[1]];orb2=gusciNC[[2]];
num1 = Min[numEperOrb[[orb1]],2*(2LperOrb[[orb1]]+1) - numEperOrb[[orb1]]];
num2 = Min[numEperOrb[[orb2]],2*(2LperOrb[[orb2]]+1) - numEperOrb[[orb2]]];
{Lorb1,Lorb2} = Sort[Flatten[{LperOrb[[orb1]],LperOrb[[orb2]]}]];
];
eneMult = {};
Which[
Length[gusciNC]==1 && Min[nEQ,2(2 LeQ+1)-nEQ]>1,
lista = {nEQeff,LeQ};
mult=Select[listadf,#[[1]][[{1,2}]] == lista&];
(* Warning  in listadf there is the total energy, not for a couple *)
ene=Table[mult[[j,2]].Fk[[orbEQ,orbEQ,All]],{j,1,Length[mult]}];
eneMult=Table[{(ene[[i]]-Min[ene]),mult[[i,1]][[4]],
mult[[i,1]][[5]]},{i,1,Length[mult]}];
eneMult= Sort[eneMult,#1[[1]] < #2[[1]] &];
,
Length[gusciNC]==2 && num1*num2 == 1,
lista = {Lorb1,Lorb2};
mult=Select[listadfdg,#[[1]][[{1,2}]] == lista&];
ene=Table[
mult[[j,2]].Fk[[orb1,orb2,All]]+ mult[[j,3]].Gk[[orb1,orb2,All]]
,{j,1,Length[mult]}];
eneMult=Table[{(ene[[i]]-Min[ene]),mult[[i,1]][[3]],
mult[[i,1]][[4]]},{i,1,Length[mult]}];
eneMult= Sort[eneMult,#1[[1]] < #2[[1]] &];
];
listaMult = Table[
""^ToString[2 eneMult[[j,3]] +1]namesMult[[eneMult[[j,2]]+1]],{j,1,Length[eneMult]}];

If[Length[eneMult] == 0,
Print["No computable structure "];,
Print[Table[{listaMult[[i]],eneMult[[i,1]]*au2Ry*Ry2cm *"cm"^ToString[-1], "=" ,eneMult[[i,1]]"au" } ,{i,1,Length[eneMult]}]//TableForm];
];
Return[eneMult];];

The loop

loopHartree[maxIter_,tolEin_,quantumNum_,grid_,psiIn_,alpha_,stampa_]:=
Module[{psi,psiold,YHartree,Vintvec,vF,Kin,matS,vec8,baseV,au2Ry,
ham,vMat,eigH,eigV,statiordinati,
x,x0,h,rho,rho2,E0,ene,idiag,iPdiag,iMdiag,deltaE,stato,potvec,
Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,tolE,iter},
(* nQ = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}; *)
{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}=quantumNum;
{x,x0,h,rho,rho2} = grid;
{Kin,matS,vec8} = InitMatrici[quantumNum,grid];
psi = psiIn;
psiold = psiIn;
baseV = 1/x[[1]];
E0 =Table[0.,{norbitals}];
ene = E0;
idiag = Table[{i,i},{i,1,Length[x]}];
iPdiag = Table[{i,i+1},{i,1,Length[x]-1}];
iMdiag = Table[{i,i-1},{i,2,Length[x]}];
iter = 1;
deltaE = 1.;
au2Ry = 2.0;
tolE = tolEin/Z^2;
(* --------------------------- *)
While[ iter ≤  maxIter && deltaE > tolE/1.,
YHartree=calcpotHartree[psi,LperOrb,x,x0,h];
Vintvec =Sum[numEperOrb[[k]]*YHartree[[k]],{k,1,norbitals}];
vF = 1./Z*Table[Vintvec - YHartree[[k]],{k,1,norbitals}];
Do[
stato =  numrad[[i]];
potvec  = rho2*(baseV-1./x + 0.5* LperOrb[[i]]*(LperOrb[[i]]+1)/x^2) + vec8;
potvec = potvec + rho2*vF[[i]];
vMat = SparseArray[{idiag-> 10./12.potvec,iPdiag->1./12.Drop[potvec,1],iMdiag->1./12.Drop[potvec,-1]},
{Length[x],Length[x]}];
ham = Kin + vMat;
{eigH,eigV} = Eigensystem[{ham,matS}]; 
statiordinati = Ordering[eigH];
psi[[i]] = eigV[[ statiordinati[[stato]] ]]; 
psi[[i]] = Sign[psi[[i]][[2]]]*psi[[i]]/Sqrt[integ[psi[[i]]^2*rho2,h]];
If[iter > 1 ,
psi[[i]] = (1.-alpha[[i]])*psi[[i]] + alpha[[i]]*psiold[[i]];
psi[[i]] = Sign[psi[[i]][[2]]]*psi[[i]]/Sqrt[integ[psi[[i]]^2*rho2,h]];
,];
psiold[[i]] = psi[[i]];
ene[[i]] = eigH[[statiordinati[[stato]] ]]-baseV;
,
{i,1,norbitals}];
deltaE = Tr[Abs[E0-ene]];
E0 = ene;
If[stampa == True,
Print["Z = ",Z," ene = ",NumberForm[(Z^2)*ene*au2Ry,{10,3},
NumberPadding-> {""," "},NumberSigns-> {"-"," "}]," Ry ", " it.= ",iter];];
iter = iter + 1;];
{ene *(Z^2),psi,iter-1}]

integralePerEpsilon[psi_,Yk_,rho2_,h_,ia_,ib_,fk_,gk_,numEperOrb_,LperOrb_]:=
Module[{fatt0,fattk,ris},
fatt0 = integ[Yk[[ib,ib,1]]*psi[[ia]]*psi[[ib]]*rho2,h]*(
1+ numEperOrb[[ib]]*gk[[ia,ib,1]]);
fattk = Sum[ 
integ[Yk[[ib,ib,k+1]]*psi[[ia]]*psi[[ib]]*rho2,h]*(
(numEperOrb[[ib]]-1)*fk[[ib,k+1]]- numEperOrb[[ib]]*gk[[ia,ib,k+1]]),{k,2,2*LperOrb[[ib]],2}];
ris = (fatt0 - fattk)/(numEperOrb[[ib]] - numEperOrb[[ia]]);
ris];

testEpsilon[ia_,ib_,numEperOrb_,LperOrb_,fk_,gk_]:=
Module[{test0,test1,test2,test},
test0 = {1+ numEperOrb[[ib]]*gk[[ia,ib,1]],
          1+ numEperOrb[[ia]]*gk[[ia,ib,1]]};
test1 = Table[(numEperOrb[[ib]]-1)*fk[[ib,k+1]]- numEperOrb[[ib]]*gk[[ia,ib,k+1]],
{k,2,2*LperOrb[[ib]],2}];
test2 = Table[(numEperOrb[[ia]]-1)*fk[[ia,k+1]]- numEperOrb[[ia]]*gk[[ia,ib,k+1]],
{k,2,2*LperOrb[[ia]],2}];
test = Join[test0,test1,test2];
Max[Abs[test]]==  0];

EpsilonOffD[quantumNum_,grid_,psi_,Yk_,Yorb_,Xorb_,fk_,gk_,
vec8_,idiag_,iPdiag_,iMdiag_,Kin_]:=
Module[
{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
x,x0,h,rho,rho2,indmom,lambda,eps0,eps,dE,ia,ib,fatt0,fattk,
vMat,b1X,potvec,potvec0,lvec,ham,term1,term2},

{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC} = quantumNum;
{x,x0,h,rho,rho2} = grid;
indmom =  Table[Flatten[Position[LperOrb,k]],{k,0,Max[LperOrb]}];
dE = Table[0.,{i,1,norbitals},{j,1,Length[x]}];
eps0 = Table[0.,{i,1,norbitals},{j,1norbitals}];
eps = eps0;
Do[
lvec = indmom[[LperOrb[[gusciNC[[iNC]]]]+1]];
Do[
ib = gusciNC[[iNC]];
ia = lvec[[irun]];
Which[
MemberQ[gusciNC,ia] == False  ,

eps[[ia,ib]] =
integralePerEpsilon[psi,Yk,rho2,h,ia,ib,fk,gk,numEperOrb,LperOrb];
eps[[ia,ib]]  = eps[[ia,ib]] /Z;
eps[[ib,ia]] = eps[[ia,ib]];,

MemberQ[gusciNC,ia] == True && ia ≠ ib &&
testEpsilon[ia,ib,numEperOrb,LperOrb,fk,gk] == False,

potvec0  = rho2*(-1./x + 0.5* LperOrb[[ia]]*(LperOrb[[ia]]+1)/x^2) + vec8;
potvec = potvec0 + rho2*Yorb[[ia]];
vMat = SparseArray[{idiag-> potvec},{Length[x],Length[x]}];
ham = Kin + vMat;
b1X = Xorb[[ia]]*rho2;
term1 = integ[psi[[ib]]*(ham.psi[[ia]] + b1X ),h];

potvec = potvec0 + rho2*Yorb[[ib]];
vMat = SparseArray[{idiag-> potvec},{Length[x],Length[x]}];
ham = Kin + vMat;
b1X = Xorb[[ib]]*rho2;
term2 = integ[psi[[ia]]*(ham.psi[[ib]] + b1X ),h];
eps[[ia,ib]] =  0.5*(term1+term2);
eps[[ib,ia]]=eps[[ia,ib]];,
True,];
,{irun,1,Length[lvec]-1}];
,{iNC,1,Length[gusciNC]}];
lambda = Inner[Times,eps,numEperOrb,List];
 dE = dE+lambda.Transpose[Transpose[psi]*rho2]; 
Return[{dE,lambda}]];

updatepsi[psiold_,psi_,alpha_,rho2_,h_]:=
Module[{norma,psiout,psioldout},
norma = integ[rho2*psi^2,h];
psiout = psi/Sqrt[norma];
psiout = (1.-alpha)*psiout + alpha*psiold;
norma = integ[rho2*psiout^2,h];
psiout=psiout/Sqrt[norma];
psioldout = psiout;
Return[{psioldout,psiout}]];

contazeri[f_]:=
Module[{nf,diff1,diff2,lista,valmax,maxmax,pts,ptmax,numzeri},
nf = Length[f];
diff1 = f[[Range[2,nf]]]- f[[Range[nf-1]]];
diff2 = diff1[[Range[2,nf-1]]];
lista=Flatten[Position[Sign[diff1[[Range[nf-2]]]]*Sign[diff2],_?(# < 0 &)]]  +1;
valmax = Abs[f[[lista]]];
maxmax = Max[valmax];
pts = Transpose[{lista,valmax}];
ptmax = Last[ Select[pts,#1[[2]]== maxmax &]];
numzeri =Length[ Select[pts,
(#1[[1]] ≤ ptmax[[1]] )|| 
(#1[[2]] > 0.1 ptmax[[2]] && #1[[1]] > ptmax[[1]])&]];
Return[numzeri]]

risolviSistemaPsi[L1_,b1X_,rho_]:=
Module[{numzeri,psi},
 psi=LinearSolve[L1,-b1X ];
numzeri = contazeri[psi];
{numzeri,psi}];

calcolaPsi[matS_,L1_,b1X_,enein_,numZ_,rho_]:=
Module[{ene,numzeri,numzerinew,psi,L1New,rapp,dE,iterzeri,segno},
 ene = -Abs[enein]; 
{numzeri,psi} = risolviSistemaPsi[L1,b1X,rho];
If[numzeri == numZ, Return[{numzeri,ene,psi}]];
rapp = 0.5;
segno =  Sign[numzeri-numZ];
dE = rapp*Abs[ene];
iterzeri = 0;
L1New = L1;
While[iterzeri < 20,
If[numzeri > numZ,ene = ene - dE;
L1New = L1New + dE*matS;,
ene = ene + dE;L1New = L1New - dE*matS;];

{numzerinew,psi} = risolviSistemaPsi[L1New,b1X,rho];
Which[
numzerinew == numZ,
Return[{numzerinew,ene,psi}];,

numzerinew > numZ && numzeri > numZ,  ,

numzerinew > numZ && numzeri < numZ, dE = + 1/2 dE;,

numzerinew < numZ && numzeri < numZ, ,

numzerinew < numZ && numzeri > numZ,dE =  +1/2 dE;
];
dE = Min[dE,Abs[ene]/2];
numzeri = numzerinew;

iterzeri = iterzeri +1;
];
Return[{numzeri,ene,psi}]];

myProd[f1_,f2_]:= integ[f1*f2*grid[[5]],grid[[3]]];
myNormalize[vec_,ip_]:=Module[{ipvec} ,
ipvec = ip[vec, vec]; vec/Sqrt[ipvec]];
myProjection[v1_,v2_,ip_]:= Module[{ipv2},
ipv2= ip[v2,v2]; ip[v1, v2] v2/ipv2];
myMultipleProjection[v1_, vecs_, ip_] := 
    Plus @@ Map[myProjection[v1, #, ip]&, vecs];
myGS[vecs_List,ip_]:=Module[{ans},
ans = Fold[Join[#1,
  {#2 - myMultipleProjection[#2, #1, ip]}]&,
                {}, vecs];
 myNormalize[#1,myProd]&/@ans ];

loopHartreeFock[maxIter_,tolEin_,quantumNum_,grid_,psiIn_,beta_,fk_,gk_,dfCoppia_,stampa_,primoE_:1]:=
Module[{psi,psiold,Yorb,Xorb,Kin,matS,vec8,baseV,
ham,vMat,eigH,eigV,statiordinati,ene,deltaE,E0,dE,
b1X,b1Xtemp,indmom,iter,numeroMax,alpha,au2Ry,
Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
x,x0,h,rho,rho2,tolE,idiag,iPdiag,iMdiag,indOrto,Ry,psiOrtog,
Yk,eps,stato,potvec},
(* nQ = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}; *)
{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}=quantumNum;
{x,x0,h,rho,rho2} = grid;
{Kin,matS,vec8} = InitMatrici[quantumNum,grid];
psi = psiIn;
psiold = psiIn;
psiOrtog = psiIn;
alpha = beta;
E0 = Table[1./(2*numPrinc[[i]]^2),{i,1,norbitals}];
numeroMax = Table[0,{norbitals}];
ene = E0;
au2Ry = 2.0;

If[norbitals == 1,
{ene,psi,iter}=loopHartree[maxIter,tolEin,quantumNum,grid,psi,beta,stampa];
Return[{ene,psi,0.*E0,iter}];];

tolE = tolEin/Z^2;

idiag = Table[{i,i},{i,1,Length[x]}];
iPdiag = Table[{i,i+1},{i,1,Length[x]-1}];
iMdiag = Table[{i,i-1},{i,2,Length[x]}];
iter = 1;
deltaE = 1.;
Ry = 2.0 "Ry";
indmom =  Table[Flatten[Position[LperOrb,k]],{k,0,Max[LperOrb]}];
indOrto = Select[indmom,Length[#]>0&];
(* --------------------------- *)
While[ iter ≤  maxIter +1 && ( deltaE > tolE/1. || Max[ene]>0),
Yk=calcpotYk[psi,LperOrb,x,x0,h];
Yorb=calcpotYorb[psi,Yk,fk,gk,dfCoppia,quantumNum];
Xorb=calcpotXorb[psi,Yk,fk,gk,quantumNum];

{dE,eps}=EpsilonOffD[quantumNum,grid,psi,Yk,Yorb,Xorb,fk,gk,
vec8,idiag,iPdiag,iMdiag,Kin];

(* Computation of energy *)

Do[
stato = numrad[[i]];
potvec  = rho2*(-1./x + 0.5* LperOrb[[i]]*(LperOrb[[i]]+1)/x^2) + vec8;
potvec = potvec + rho2*Yorb[[i]];
vMat = SparseArray[{idiag-> 10./12.potvec,iPdiag->1./12.Drop[potvec,1],iMdiag->1./12.Drop[potvec,-1]},
{Length[x],Length[x]}];
ham = Kin + vMat;

b1Xtemp = Xorb[[i]]*rho2 - dE[[i]];
b1X = 10./12.b1Xtemp + 1./12.*( Flatten[{0.,b1Xtemp[[Range[2,Length[x]]]]}] +
Flatten[{b1Xtemp[[Range[1,Length[x]-1]]],0.}]);
ene[[i]] = (psi[[i]].(ham.psi[[i]] + b1X))/(psi[[i]].(matS.psi[[i]]));
,{i,1,norbitals}];

deltaE = Total[Abs[E0-ene]]; 
E0 = ene;
If[stampa == True,
Print["Z = ",Z," ene = ",NumberForm[(Z^2)*ene*au2Ry,{10,5},
NumberPadding-> {""," "},NumberSigns-> {"-"," "}]," Ry ", " it.= ",iter];];

(* loop on orbitals *)

Do[
stato = numrad[[i]];
potvec  = rho2*(-1./x + 0.5* LperOrb[[i]]*(LperOrb[[i]]+1)/x^2) + vec8;
potvec = potvec + rho2*Yorb[[i]];
vMat = SparseArray[{idiag-> 10./12.potvec,iPdiag->1./12.Drop[potvec,1],iMdiag->1./12.Drop[potvec,-1]},
{Length[x],Length[x]}];
ham = Kin + vMat;
L1 = ham - ene[[i]]*matS;
b1Xtemp = Xorb[[i]]*rho2 - dE[[i]];
b1X = 10./12.b1Xtemp + 1./12.*( Flatten[{0.,b1Xtemp[[Range[2,Length[x]]]]}] +
Flatten[{b1Xtemp[[Range[1,Length[x]-1]]],0.}]);
(* psi[[i]] = LinearSolve[L1,-b1X ];  *)
{numeroMax[[i]],ene[[i]],psi[[i]]}= calcolaPsi[matS,L1,b1X,ene[[i]],numrad[[i]],rho];
psi[[i]] = Sign[psi[[i]][[2]]]*psi[[i]]/Sqrt[integ[psi[[i]]^2*rho2,h]];
{psiold[[i]],psi[[i]]}=updatepsi[psiold[[i]],psi[[i]],alpha[[i]],rho2,h]; 
,
{i,primoE,norbitals}];
(* psiOrtog[[indOrto[[k]] ]] = GramSchmidt[psi[[indOrto[[k]]]],
InnerProduct-> ((integ[(#1*#2) *rho2,h])&)];, *)
(* Do[
psiOrtog[[indOrto[[k]] ]] = Orthogonalize[psi[[indOrto[[k]]]],((integ[(#1*#2) *rho2,h])&)];,
{k,1,Length[indOrto]}]; *)
 Do[
psiOrtog[[indOrto[[k]] ]] = myGS[psi[[indOrto[[k]]]],myProd];,
{k,1,Length[indOrto]}]; 
Table[psi[[i]] = Sign[psiOrtog[[i]][[2]]]*psiOrtog[[i]],{i,primoE,norbitals}];
psiold = psi;
iter = iter + 1;];
If[Max[Abs[numeroMax - numrad]]> 0,iter = maxIter+1;Print["Warning: No convergence"];];
{(Z^2)*ene,psi,(Z^2)*eps,iter-1}];

cycleHartreeFock[LSlevel_List:{1},opts___]:=
Module[{ene,psi,eps,cycles,cycles2,eneRy,au2Ry,
Lmultipletto,Smultipletto,fk,gk,dfCoppia,eneHartree,
nomeTermine,Z,enetot,enekin,viriale,alpha,quantumNum,
maxIter,tolE,h,x0,psiIn,primoE,stampa,initHF,hLoc,hDef,gridOld},

{maxIter,tolE,h,x0,psiIn,alpha,initHF,primoE,stampa,hDef,gridOld} = 
{cyclesHF,errHF,hgrid,x0grid,psiInput,alphaHF,inizioHF,firstEHF,printHF,hPerHartreeHF,gridOldHF}/.{opts}/.Options[cycleHartreeFock];

quantumNum = numQ$HF;
au2Ry = 2.0;


(* With a single orbital the Hartree method is used *)
If[numQ$HF[[2]]== 1,h = hDef];


{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];
If[Lmultipletto<0,nomeTermine = "Emed",
nomeTermine=""^ToString[2 Smultipletto +1]namesMult[[Lmultipletto+1]]];Z = quantumNum[[1]];
If[stampa == True,Print["Z = ",Z,"; LSlevel: ",nomeTermine];];


If[Length[alpha]== 0,
alpha=0.5*Initalpha[quantumNum];];

(* here we define the global variable grid *)
(* NOTE for psiIn: usually gridOld is the current grid, if psiIn ha sto be used more than once
the reader must supply gridOld to the function *)

If[initHF == True,
hLoc = Max[hDef,h];
Unprotect[grid];
grid = parametrigrid[quantumNum,hLoc,x0];
psi = InitPsi[quantumNum,grid];
{ene,psi,cycles}=loopHartree[10,tolE,quantumNum,grid,psi,alpha,False];
psi= updategridPsi[x0,h,quantumNum,grid,psi];
,
If[Length[psiIn] ==  0,
Unprotect[grid];
If[Length[gridOld]== 0,
grid = parametrigrid[quantumNum,h,x0];,
grid = gridOld;
h = grid[[3]]; x0 = grid[[2]];];
Protect[grid];
psi = InitPsi[quantumNum,grid];
,
If[Length[gridOld]== 0,gridOld = grid;];
psi= updategridPsi[x0,h,quantumNum,gridOld,psiIn];];
];

{ene,psi,eps,cycles}=loopHartreeFock[maxIter,tolE,quantumNum,grid,psi,alpha,fk,gk,dfCoppia,stampa];

(* If the procedure do not converge an attempt is made with Hartree equation  *)

If[cycles ≥ maxIter -1 && Length[psiIn] ==  0 && inizioHF == False,
hLoc = Max[hDef,h];
Print["switch to Hartree "];
Unprotect[grid];
grid = parametrigrid[quantumNum,hLoc,x0];
Protect[grid];
psi = InitPsi[quantumNum,grid];
{eneHartree,psi,cycles}=loopHartree[10,tolE,quantumNum,grid,psi,alpha,False];
psi= updategridPsi[x0,h,quantumNum,grid,psi];
{ene,psi,eps,cycles2}=loopHartreeFock[maxIter,tolE,quantumNum,grid,psi,alpha,fk,gk,dfCoppia,stampa];
cycles += cycles2;
If[Max[Abs[(ene-eneHartree)/ene]]>0.5,
Print["Warning: Suspect solution"];
Print["change h or use InitHartre with a small h, or change h$Hartree"]; ];
];

{enetot,enekin,viriale}=calcEne[psi,grid,quantumNum,fk,gk,dfCoppia];

If[stampa == False,
Print["Z = ",Z," ene = ",NumberForm[ene*au2Ry,{10,5},
NumberPadding-> {""," "},NumberSigns-> {"-"," "}]," Ry ", " it.= ",cycles];];

Print[namesEl[[Z]]," Z=",Z,";  ",nomeTermine,":  ",
"E = ",NumberForm[enetot "au",{10,7}],"; Kin = ",
NumberForm[enekin "au",{10,7}],"; vir = ",
NumberForm[viriale,{10,7}],";"];
(* global variable *)
Unprotect[epsLSlevel,virialeLSlevel];
epsLSlevel = eps;
virialeLSlevel = viriale;
Protect[epsLSlevel,virialeLSlevel];
Return[{ene,psi,enetot}];
]

Options[cycleHartreeFock] = 
{cyclesHF-> 50,errHF-> (10.^(-5)),hgrid-> 0.1,x0grid-> 0.1,psiInput-> {},alphaHF-> {},inizioHF-> False,firstEHF-> 1, printHF -> True,hPerHartreeHF-> 0.1,gridOldHF-> {}};

calcCore[opts___]:=
Module[{ene,psi,psiTot,cycles,eps,ene1,
fkCore,gkCore,dfCoppiaCore,au2Ry,
beta,gusciC,nQCore,quantumNum,enetot,enekin,viriale,
Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
maxIter,tolE,h,x0,alpha,stampa,siEneHF},

{maxIter,tolE,h,x0,alpha,stampa,siEneHF} = 
{cyclesHF,errHF,hgrid,x0grid,alphaHF,printHF,si$EneHF}/.{opts}/.Options[calcCore];

quantumNum = numQ$HF;
{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}=quantumNum;

(* This is the grid of the full atom *)

Unprotect[grid];
grid = parametrigrid[quantumNum,h,x0];
Protect[grid];

If[Length[alpha]== 0,
alpha=Initalpha[quantumNum];];

 ene=Table[1./(2*numPrinc[[i]]^2),{i,1,norbitals}];
au2Ry = 2.0;

gusciC=Flatten[Position[Positive[2*(2*LperOrb+1) -  numEperOrb],False]];
nQCore = quantumNum;
Table[nQCore[[i]]=quantumNum[[i,gusciC]],{i,3,7}];
nQCore[[8]] = {};nQCore[[2]]=quantumNum[[2]]-Length[numEperOrb[[gusciNC]]];
beta =  0.5*alpha [[gusciC]];

psiTot = InitPsi[quantumNum,grid];
psi = psiTot[[gusciC]];

{fkCore,gkCore}=fgMedi[nQCore[[2]],nQCore[[4]]];
dfCoppiaCore = 0.;

{ene1,psi,eps,cycles}=loopHartreeFock[maxIter,tolE,nQCore,grid,psi,beta,fkCore,gkCore,dfCoppiaCore,stampa];
If[stampa == False,
Print["Z = ",Z," ene = ",NumberForm[ene1*au2Ry,{10,5},
NumberPadding-> {""," "},NumberSigns-> {"-"," "}]," Ry ", " it.= ",cycles];];

psiTot[[gusciC]] = psi;
ene[[gusciC]] = ene1;
enetot = 0.;
If[siEneHF,
{enetot,enekin,viriale}=calcEne[psi,grid,nQCore,fkCore,gkCore,dfCoppiaCore];
Print[namesEl[[Z]]," Z=",Z,";  ",
"E = ",NumberForm[enetot "au",{10,7}],"; Kin = ",
NumberForm[enekin "au",{10,7}],"; vir = ",
NumberForm[viriale,{10,7}],";"];
];
Return[{ene,psiTot,enetot}];
]

Options[calcCore] = 
{cyclesHF-> 20,errHF-> (10.^(-5)),hgrid-> 0.1,x0grid-> 0.1,alphaHF-> {}, printHF -> False,si$EneHF-> False};

Graphics

plotPsi[psi_,zoom_,statiIn_:{}]:=Module[
{Z,numPrinc,sigma,norbitals,lplot,rho,j,quantumNum,x,stati,datapsi},
(* grid = {x,x0,h,rho,rho2} ;
 nQ = {Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC}; *)
quantumNum = numQ$HF;
{x,rho} = grid[[{1,4}]];
{Z,norbitals,numPrinc,sigma} = quantumNum[[{1,2,5,7}]];
If[Length[statiIn]== 0,stati = Range[1,norbitals];,stati = statiIn;];
lplot = Table[0,{i,1,Length[stati]}];
Block[{$DisplayFunction = Identity},
Do[
j = stati[[istat]];
datapsi = Transpose[{x/Z,(Sqrt[Z])(Sqrt[rho])psi[[j]]}];
lplot[[istat]]=ListPlot[datapsi,PlotRange-> {{0,1/zoom*numPrinc[[j]]/sigma[[j]]},All},
PlotStyle-> Hue[1-0.9istat/norbitals],
Joined-> True,PlotLabel-> StringJoin["orbital  ",electronsnames[[j]]]];,
{istat,1,Length[lplot]}]]; lplot];

Data

eneFF=-{0.5`,2.86168`,7.4327269`,14.573023`,24.529061`,37.688619`,54.400934`,74.809398`,99.409349`,128.5471`,161.85891`,199.61463`,241.87671`,288.85436`,340.71878`,397.5049`,459.48207`,526.81751`,599.16479`,676.75818`,759.73572`,848.406`,942.88433`,1043.35637593,1149.8662`,1262.4437`,1381.4146`,1506.8709`,1638.9501`,1777.8481`,1923.261`,2075.3587`,2234.2386`,2399.8676`,2572.4413`,2752.055`,2938.3574`,3131.5457`,3331.6842`,3538.9951`,3753.59772675,3975.54949849,4204.7887`,4441.53948671,4685.88170311,4937.92102287,5197.69847186,5465.13314125,5740.1691`,6022.9317`,6313.4853`,6611.784`,6917.9809`,7232.1384`,7553.9336`,7883.5438`,8221.0638`,8566.9196`,8921.181`,9283.8829`,9655.099`,10034.953`,10423.543`,10820.617`,11226.568`,11641.453`,12065.29`,12498.153`,12940.174`,13391.456`,13851.808`,14321.25`,14799.813`,15287.546`,15784.533`,16290.649`,16806.113`,17330.949`,17865.212`,18408.991`,18961.825`,19524.008`,20095.586`,20676.501`,21266.882`,21866.772`};

eneFFIone = {0.,-2.86167999,-7.23641520,-14.27739481,-24.23757517,
-37.29222375,-53.88800500,-74.37260566,-98.83172017,
-127.81781408,-161.67696256,-199.37180959,-241.67467032,
-288.57313096,-340.34977577,-397.17318261,-459.04859050,
-526.27453408,-599.01757919,-676.57001233,-759.53914376,
-848.20340052,-942.67078342,-1043.13939315,-1149.64938238,
-1262.21301160,-1381.12874921,-1506.59109851,-1638.72824195,
-1777.56754417,-1923.05972188,-2075.08649072,-2233.88833396,
-2399.55857321,-2572.04520991,-2751.56739321,-2938.21993041,
-3131.37377574,-3331.47288120,-3538.80930393,-3753.38951240,
-3975.33370152,-4204.59435845,-4441.32195495,-4685.66417041,
-4937.67592880,-5197.48133281,-5464.87860804,-5739.97839093,
-6022.67832166,-6313.16593964,-6611.50339297,-6917.62727171,
-7231.70894562,-7553.81032698,-7883.39188624,-8220.83156343,
-8566.61237239,-8920.82088049,-9283.51980928,-9654.73624604,
-10034.54319129,-10423.06873095,-10820.41331457,-11225.91478600,
-11641.12828547,-12064.97790897,-12497.85493186,-12939.83196552,
-13391.04630300};

These are the experimental values of energy shifts from the ground state, in atomic units, for elements from Z=1 (H) to Z=24 (Cr). The values have been averaged over the fine structure splitting.  Original data can be found at

expδELS ={{0},{0},{0},{0},{0.000046435131339298695`},{0.00013482582288895526`,0.04644117607740069`,0.09863563667846192`},{0,0.08760898277824589`,0.131400960344475`},{0.0003552787225555736`,0.07229929901653695`,0.15397033720473138`},{0.0006138006284619452`},{0},{0},{0},{0.0003403916565063879`},{0.0006819927668926473`,0.02869965863784401`,0.0701418742743549`},{0,0.051807290569293174`,0.08538379378813932`},{0.0008919304699988349`,0.0420941998731694`,0.10105930631480378`},{0.0013400941367577093`},{0},{0},{0},{0.0004602080858698102`},{0.0010138525048312868`,0.033057825201456094`,0.03894430861737688`,0.055215465789039964`},{0.0014550266086668633`,0.04426671064137332`},{0,0.034596982587850066`,0.09350255949688818`}};

Computation

Computation of the energies with Hartree Fock method

Let us define the conversion between atomic units and Ry.

au2Ry = 2;
Ry2cm = 109737.31568525;

The electron configuration

Here we define the electronics configuration of Carbon, using Mathematica function ElementData

degIoniz = 0;
electrons = ElementData["C","ElectronConfiguration"]

We extract the quantum numbers from the configuration. namesEL is a list of names, defined in the section lists above.

quantumNum = getquantumNum[electrons,degIoniz];

Zatom=quantumNum[[1]];
If[Zatom ≤  Length[namesEl],nameAtom = namesEl[[Zatom]],nameAtom =" ? "];
Print[nameAtom , "   Z = ", Zatom,"  conf.:  ",  ElementData[Zatom,"ElectronConfigurationString"]];

The level

listLS=choiceMult[];

LSlevel = {ground};

Energy computation

{ene,psi,enetot}=cycleHartreeFock[LSlevel,hgrid-> 0.05];//Timing

An equivalent statement is cycleHartreeFock[], as the ground state is the default. We can delete intermediate results from the output by selecting the option printHF, as an instance

{ene,psi,enetot}=cycleHartreeFock[printHF-> False];//Timing

The electrons

The energies of the orbitals with more digits

NumberForm[TableForm[Transpose[{electronsnames,
ene*au2Ry*" Ry"}],TableAlignments->Right],{12,7},
NumberPadding-> {""," "},NumberSigns-> {"-"," "}]

and the plots of the wave functions

zoom = 1;
Show/@plotPsi[psi,zoom]

Show[plotPsi[psi,0.2],PlotLabel-> "orbitals"]

Comparison of LS levels computations

We can get more informations by estimating the energy differences between various terms using ground state wave functions

diffEApprox = approxMultiplets[];

We can also compute independently the energies and compare the results

eneLS=Table[{ene2,psi,enetot2}=cycleHartreeFock[listLS[[k]],printHF-> False];
enetot2,{k,1,Length[listLS]}];

diffE = Sort[eneLS - eneLS[[1]]]

We see an agreement of a few percent:

diffEApprox[[All, 1]]

Computations and experimental data

LS levels

In the vector expδELS (see section DATA) are stored the experimental energy differences between LS levels for atoms up to Z=24.

For the differences in Carbon we have, in atomic units

expδELS[[Zatom]]- Min[expδELS[[Zatom]]]

We see that the approximation is about 20 percent or worse. The configuration mixing and the LS coupling are relevant if one look for accurate results.

Total energy

In Mathematica 6 the function ElementData contains the ionization energy for the last electrons for all elements, in particular up to Z=10 for all electrons. The energies are expressed in kiloJoule/Mole. We use an approximative conversion factor, assuming 1Ry = 0.5 au the ionization energy of H. The ``exact'' convertion factor can be given by the lines commented below if needed:

(* <<Units`
<<PhysicalConstants`
kJM2Ry=Convert[10^3 Joule/AvogadroConstant/Mole,Rydberg]/Rydberg *)

convFact = 0.5/ElementData["H","IonizationEnergies"][[1]]

In atomic units we have

mathEnergies=ElementData[Zatom,"IonizationEnergies"]*convFact

The total energy can be obtained by summing the single ionization energies:

Tr[mathEnergies]

We can compare with above results for the external electrons and for the total energy

{Min[-ene],-enetot}

We see that error on total energy is a fraction of percent.