Guide to Hartree-Fock equation

Quantum Mechanics
                     K. Konishi, G. Paffuti

In this notebook we present in some details the numerical procedures used in the solution of the Hartree-Fock equation.

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 is 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 the ground-state energy and energy of the few first excited levels of an atom in the Hartree-Fock approximation. The energies and wave functions are calculated. 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 allows 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.

Lists

The two vectors listadf and listadfdg are defined. These vectors parametrize the difference between the coefficients fk and gk for a given level and their value for the mean energy (see text and below). Some lists with the standard names for multiplets are also defined.

The coefficients are stored in two vectors listadf and listadg with the format

listadf :  {{numElettroniEQ, La, Lb, Ltot, Stot}, δfk};  
 listadfdg: {{La, Lb, Ltot, Stot}, δfk, δgk}

Let us stress that the values stored refer to the total shift in the energy, while the coefficients given above refer to the energy for a single couple.

All these variables are protected and the assignments are done with

Unprotect[var]; var = value; Protect[var];

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

Electronic configuration and multiplets

Electronic configuration

The notebook accepts two different inputs for the electronic configuration. We will store here these data in the variable electrons and take as an example the Carbon (neutral atom). It is necessary to specify the ionization degree (0 for a neutral atom). This information is required to fix the nuclear charge: as an example the Carbon configuration 1(s^2)2(s^2)2(p^2) would correspond to ionized Nitrogen if degIoniz = 1.

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

electrons = ElementData["C","ElectronConfiguration"]

If the reader do not remember the chemical symbol an alternative input is allowed by  ElementData, the atomic number Z

ElementData[6,"ElectronConfiguration"]

The electronic configuration can be printed with

ElementData["C","ElectronConfigurationString"]

From this input configuration the program extracts the quantum numbers that will be used in the rest of the functions, i.e. angular momenta, principal quantum numbers etc. It is convenient to store all these informations in a single variable, quantumNum. We will use the function

quantumNum = getquantumNum[elettrons, degIoniz]

This function calls getquantumNumA, or getquantumNumB for the two different format in input. The first format are different because it contains a string as first component. The variable quantumNum store the information in the form

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

The notations are obvious except for the last two entries.

The parameter sigma in this program is fixed with

sigma = Z - Table[ Total[numEperOrb[[Range[1, i]]]] - ( numEperOrb[[i]] + 1), {i, 1, norbitals}];

From the variable quantum number we can extract all the needed informations, eventually using a vector of indexes. A an example if Z and numEperOrb are required one can use

{Z, numEperOrb} = quantumNum[[{1, 3}]];

In the following we write the code for these functions. For Mathematica functions like Flatten, Inner etc. the reader can use the Mathematica Help.

Code

getquantumNum[electrons_,degIoniz_]:=
Module[{quantumNum,nomi,numPrinc},
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]

Example:

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

quantumNum = getquantumNum[electrons,degIoniz]

We can extract all the values stored in the variable quantumNum

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

LperOrb

or only some of them

{Z,numEperOrb,numrad} = quantumNum[[{1,3,6}]]

The function has defined the global variable numQ$HF,

numQ$HF

The program protect this variable from changes

numQ$HF = {3,5}

the variable keeps its value

numQ$HF

Multiplet

Once an electron configuration is given we have to choose the spectroscopic term to be studied. Each electron configuration gives rise in Russell-Saunders coupling to several (L,S) terms, with an energy E(L,S). This program allows the computations of these energies and the corresponding eigenfunctions within the above limitations on the electronic configuration.

The informations about the term we want to study is contained in a vector variable which will be called LSlevel below. The choice of LSlevel determines the computation to be done

E = (∑_(L.S))E(L,S)(((2L+1)(2S+1))/deg);      deg =  (∑_(L.S))(2L+1)(2S+1);

We write a function choiceMultiplet which allow the choice between different terms

listLS=choiceMultiplet[];

Once the level has been chosen the needed parameter are internally computed by the function parametriMultipletto.

{Lmultipletto, Smultipletto, fk, gk, dfCoppia} = 
parametriMultipletto[quantumNum, LSlevel, listadf, listadfdg];

The output of this function contains the L, S values, and 3 vectors of coefficients, fk, gk, dfCoppia, to be used in the computation of the energy and in the solution of Hartree Fock equation.

These coefficients can be computed directly for the mean energy and the ground state (see below) while the deviations from these values have been computed using notebooks  [NB3, NB4]. The format is

listadf :  {{numElettroniEQ, La, Lb, Ltot, Stot}, δfk};  
 listadfdg: {{La, Lb, Ltot, Stot}, δfk, δgk}

The deviations are identical for n electrons or n vacancies, for  instance Carbon and Oxygen have the same df and dg, as far as the total energy is concerned. We stress that the coefficients in the list refer to the energy, not to the energy per couple.

Details on particular functions are given below for the reader interested in extending the program. For basic use these functions are not needed.

Functions

The general form of the energy in the HF  approximation is (see text)

E = (∑_a)(Subscript[q, a])I[a] + (∑_a)(((Subscript[q, a])((Subscript[q, a])-1))/2)(∑_(k=0))(Subscript[f, k])[a](F^k)[a,a] +
                                          (∑_(a<b))(Subscript[q, a])(Subscript[q, b])(∑_(k=0))((Subscript[f, k])[a,b](F^k)[a,b]+(Subscript[g, k])[a,b](G^k)[a,b]);

The first term embodies the kinetic energy and the electron nucleus interaction.

The functions Fk and Gk are integrals on radial wave functions, q is the number of electrons in the orbital, and fk, gk depend on the level.

For the mean energy:

(Subscript[f, 0])[a]= (Subscript[f, 0])[a,b] =1;                          (Subscript[f, k])[a]= - ((2(Subscript[L, a])+1)/(4(Subscript[L, a])+1))(((({{(Subscript[L, a]), k, (Subscript[L, a])}, {0, 0, 0}})))^2);
 (Subscript[g, k])[a,b]= - (1/2)(((({{(Subscript[L, a]), k, (Subscript[L, b])}, {0, 0, 0}})))^2);            (Subscript[f, k])[a,b] =0;  k > 1;

The variations from these values are taken from the lists.

For the ground state one can consider (Hund rule) the value with the maximum possible value of Sz. If the shell has less than (2L+1) electrons one consider the particles, otherwise in the following for particles we mean the vacancies. Once fixed Sz one can assign the values of the projection Lz for the single particles: L, L-1... up to the total number of particles. The total Lz so constructed gives the L for the state. This construction fix in a unique way the state, i.e. there is a single Slater determinant with these quantum numbers. In this case the matrix element for the Hamiltonian is easily written in term of two particle operators as (see also notebooks [NB3, NB4])

(ij | U | ij) = (∑_k)c[k,(Subscript[L, i]),(Subscript[m, i]),(Subscript[L, i]),(Subscript[m, i])]c[k,(Subscript[L, j]),(Subscript[m, j]),(Subscript[L, j]),(Subscript[m, j])]F[k,(Subscript[L, i]),(Subscript[L, j])];
(ij | U | ji) = (∑_k)((c[k,(Subscript[L, i]),(Subscript[m, i]),(Subscript[L, j]),(Subscript[m, j])])^2)G[k,(Subscript[L, i]),(Subscript[L, j])];

and this allows the computation of coefficients fk and gk in general form.

The sum over k is bounded by angular momentum conservation laws to  |(Subscript[L, i]) - (Subscript[L, j])| ≤ k ≤ |(Subscript[L, i]) + (Subscript[L, j])|.

The coefficients ck are related to Clebsch - Gordan coefficients with

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

The function fmed computes  fk  for the mean energy case

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]

Hund rule

The function calcdfHund computes the variations from the mean using the coefficients c[k ] given above.

The variable numE selects between electrons and vacancies. The coefficients for the energy are computed and then the mean energy is subtracted. To be consistent with the other coefficients we define the coefficients normalized to a single pair of particles, this is done with the variable fatt.

Let us note that in the subtraction one must use the mean energy per pair times numE, while for the coefficient for a single pair one needs the true number of electrons, numTrue, given as an input in the function.

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

The function calcdFdGHund takes care of a particular case.

For an atom with an incomplete shell and one s-electron one has the same coefficients fk and gk as the complete shell but there is an additional variation in the coefficient gk for the coupling between the incomplete shell and the s electron.

For the instances s-s, s-p, s-d one can use the listed values, but for some cases one can be interested in configurations (L^n)s not covered in the lists. The coefficient depends on the total spin of the system. Here we compute the case in which the incomplete shell is in the fundamental configuration and the total spin is S+1/2 where S is the spin of the incomplete shell.

This function enters in the computation of the ground state energy for atoms or ions which fill the d shell while an s-electron is present, as 3(d^9)4s, so we are talking about a configuration with 2 incomplete shells. The value used for the coefficient can be easily checked by the reader.

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

The function trovaMultipletto simply scans the lists for the wanted couple (L,S) for a given configuration. If the reader uses an LSlevel variable not admitted the program stops and gives an error message.

The lists have the structure

listadf = {{numE,La,Lb,L,S},{dfk, k=0,6}};   listadfdg = {{La,Lb,L,S},{dfk, k=0,6},{dgk, k=0,6}}

and are used only in the following cases:

In the first case obviously the angular momentum of the two particles entering in the interaction are the same, La=Lb.

L,S are those requested by the particular term to be computed.

In the calling sequence for the function it is used a variable numEQ for the number of equivalent electrons and it is constructed a list to be found in listadfg.

In the output the function gives the value read in the lists, or a vector {} if nothing is found. This does not happens if the input is correct (i.e. the LSlevel variable is within those listed by the notebook).

The following functions simply computes again fk and gk for the mean energy or uses trovaMultipletto to write the variations df and dg. For readability they separate the case of one incomplete shell and 2 incomplete shells.

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

The function parametriMultipletto has the task of select between the different cases in the computation of the level parameters. The criterion adopted is the presence on the lists, and is controlled by the Boolean variable StatoLista.

The variable dfCoppia refers to the variations of fk and gk for two electrons both belonging to incomplete shells. As this case is rare in this program the variable is treated apart.

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,
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}];]

Finally the function choiceMult chooses, using the number of incomplete shells and to electronic configuration, between the various cases.

The variable mult gives the term selected from the lists, which must be in agreement with the vector lista built with orbitals of the configuration.

choiceMult[]:=
Module[{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
Smult,Lmult,orb,lista,listanomi,listLS,nomi,mult,nomefond,quantumNum},
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]];

Example:

A typical use is the following.

degIoniz = 0;
electrons = ElementData[6,"ElectronConfiguration"];
quantumNum = getquantumNum[electrons,degIoniz];

listLS=choiceMult[];

LSlevel = listLS[[1]]

Internally the information contained in LSlevel is first used by parametriMultipletto. Let us try a non-existent term

LSlevel = {2,1};

{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];

This does not happen if the level is selected from the list produced by choiceMult[]:

LSlevel = listLS[[3]];
{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];

Here are the quantum numbers (L, S)

{Lmultipletto,Smultipletto}

For the ground state we can write

LSlevel = {fond};
{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];
{Lmultipletto,Smultipletto}

while if we are interested to the mean energy

LSlevel = {};
{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];
{Lmultipletto,Smultipletto}

By convention we assign negative values to L and S for the mean energy.

Grid, equations's structure and initializations

The problem

For a numerical solution we need a grid of points.

The general form of the equations to be solved is

(Subscript[ℒ, a]) (Subscript[P, a])+ (Subscript[V, a]) (Subscript[P, a]) + (Subscript[W, a]) = (Subscript[E, a])(Subscript[P, a]) + (Subscript[δE, ab])(Subscript[P, b])

The first operator is

(Subscript[ℒ, a]) = -(1/2)((d^2)/(dr^2)) + (((Subscript[L, a])((Subscript[L, a])+1))/(2(r^2))) - (Z/r)

We start by rescaling the variables with  Z,  this allow a smoother variation with atomic number

r = (x/Z);    (Subscript[P, a])[r] = (Sqrt[Z]) (Subscript[Q, a])[x];   ∫(P^2)dr = ∫(Q^2)dx

In these variables

(Subscript[ℒ, a]) (Subscript[Q, a])+ (Subscript[𝒴, a]) (Subscript[Q, a]) + (Subscript[𝒳, a]) = (Subscript[ϵ, a])(Subscript[Q, a]) + (Subscript[ϵ, ab])(Subscript[Q, b]);
(Subscript[ℒ, a]) = -(1/2)((d^2)/(d (x^2))) + (((Subscript[L, a])((Subscript[L, a])+1))/(2 (x^2))) - (1/x);    
(Subscript[𝒴, a]) = ((Subscript[V, a])/Z); (Subscript[𝒳, a])= ((Subscript[W, a])/Z);  (Subscript[E, a]) = (Z^2) (Subscript[ϵ, a]); (Subscript[δE, ab]) = (Z^2) (Subscript[ϵ, ab]);

NOTE:  A factor 1/Z appear in the e-e interaction terms. Energies in usual atomic units are obtained from the rescaled results with a multiplication by (Z^2).

The equation could be discretized on a grid with step h in the form  (Subscript[x, j])= j h but different wave functions have quite different characteristic scale and this would require a very large number of points in the grid. In practice the grid must be finer near the origin, to fix correctly the boundary conditions, but can be coarser for large distances where the behavior is quite regular and the wave function is almost vanishing. The simplest solution is to make a change of variables to pass to an exponential grid. We pose

x = x0 (Exp[ξ] - 1);           (Subscript[Q, a])[x] = (Sqrt[ρ])(Subscript[ψ, a])[ξ];        ρ = x0 Exp[ξ];

We take a grid of points in the ξ variable

(Subscript[ξ, j]) = j h;  j = 1, ... N;

h is the grid step, N the numbers of points and x0 is a free parameter which, with h, can be used to vary the grid. For small x the points behave as x0 h. We have introduced the useful variable ρ = x0 Exp[ξ] = x + x0, which appears often in the formulas.

The vector of grid points has an index starting from 1 to follow Mathematica conventions. Let us note the origin do NON belong to the grid.

N is taken even as this is the convention used in the computation of the potential functions (see below).

In ξ variables the equation is

-(1/2)((d^2)/((dξ)^2))(Subscript[ψ, a]) + ((1/8) +(ρ^2)((((Subscript[L, a])((Subscript[L, a])+1))/(2 (x^2))) - (1/x))) (Subscript[ψ, a]) + (ρ^2)(Subscript[𝒴, a])(Subscript[ψ, a]) + (ρ^2)(Subscript[𝒳, a])=  (ρ^2)(Subscript[ϵ, a])(Subscript[ψ, a])+(ρ^2)(Subscript[ϵ, ab])(Subscript[ψ, b])

x and ρ are function of ξ.

We used the relations

dx = ρ dξ;    (d/dx)=(1/ρ)(d/(dξ)); ((d^2)/(dx^2))=(1/(ρ^2))(((d^2)/((dξ)^2))-(d/(dξ)));

Let us note that the (integral) scalar product is now written

∫(Subscript[Q, a])(Subscript[Q, b])dx = ∫(Subscript[ψ, a])(Subscript[ψ, b])(ρ^2) dξ

In this formula a factor ρ is due to dx, the second one to the definition of ψ.

Let us note that in the following the ξ variable will appear only as an index for the grid points, F[[i]] will be the value of the function F at the i-th point of the grid, i.e. (Subscript[ξ, i]). This is automatically the value of F in the corresponding x point.

Discretization

We need a discretization for the derivatives and a method to compute integrals.

For the discretization we will use the following approximation, see also [NB7]:

(δ^2)y[j]= (h^2)(1 + (1/12)(δ^2))y''[j];   (δ^2) f[j]= f[j+1]- 2 f[j] + f[j-1];

Substituting for the second derivative the equation we get

-(1/(2(h^2)))(δ^2)ψ +(1 + (1/12)(δ^2))( (1/8)(Subscript[ψ, a]) + (ρ^2)(Subscript[V, a])(Subscript[ψ, a]) + (ρ^2)(Subscript[𝒳, a])) = (1 + (1/12)(δ^2)) (ρ^2)(Subscript[ϵ, a])(Subscript[ψ, a]) + (1 + (1/12)(δ^2))(ρ^2)(Subscript[ϵ, ab])(Subscript[ψ, b])

We have written a single potential term for the coulombic and centrifugal terms.  (δ^2) acts as a matrix on vectors ψ, potential etc.

For the computation of integrals we will use a modification of the Simpson's formula (see  [NB8]):

(∫_(Subscript[x, j-1])^(Subscript[x, j+1]))f[x]dx = 2h(f[j] +(1/6)(δ^2)f[j]-(1/180)(δ^4)f[j])+O((h^7))=
            (h/90)(-f[j+2] + 34 f[j+1] + 114 f[j] + 34 f[j-] - f[j-2])+O((h^7));

For the first interval :

(∫_(Subscript[x, 0])^(Subscript[x, 2]))f[x]dx = (h/3)(f[0]+4 f[1] + f[2])

All functions are exponentially vanishing and there is no need to correct for the last interval.

We will write three functions:

parametrigrid[quantumNum, h, x0]; 
InitMatrici[quantumNum, grid]; 
integ[f, h];

The only one which is useful to the reader is the first one, the last two just define some matrices needed below and the previous formula for the integration

In the function InitMatrici one defines a matrix matS which is in effect (1+(δ^2)/12) rho2, a matrix Kin for the kinetic term in the HF equation, i.e. - 1/2 (δ^2)/h^2 and a vector vec8 which will be used to build a diagonal matrix filled with 1/8 in the HF equation. In Mathematica 6 one can also use the function Band[] for this task.

function parametrigrid

This function requires the choice of the parameters h and x0.

The suggested values are h = 1/10;  x0 = 0.1 to avoid instabilities near the origin due to the singular nature of the potential.
If a finer grid is requested to get more precise results the function updategridPsi can be used, this function will be presented below.

The function parametrigrid gives in output a list of the form

grid = {x, x0, h, rho, rho2};

where  ρ2 = (ρ^2);  x is the vector of points, h and x0 are the grid parameters of the input. The quantum numbers are needed because the largest point on the grid is chosen as  xmax = 40.*Max[numrad]^2 * Z. This appears to be a reasonable value, but can be modified by the user. In  the computation of potential functions by convention the number of points is assumed to be even, this explain the call to the function Mod in the code below.

To get a finer grid on which to extend a computed wave function one can use

{grid,psi} = updategridPsi[x0, h, quantumNum, gridIn, psiIn]

The input needs the new parameters x0, h, the quantum numbers of the configuration, the old grid, in the variable gridIn, and the old wave functions vector, with the variable psiIn.

The wave functions will be stored as a matrix ψ. The i-th component ψ[[i]] is the wave function of the orbital i. In updategridPsi the old function is interpolated on the new grid using the Mathematica function  Interpolation. As result we have the couple {grid,psi} with the new grid and the new ψ.

Hydrogen functions and screening parameters

Finally in the iterative procedure we need a starting expression for wave functions. We will use hydrogen wave functions computed with a screening parameter sigma. The parameter sigma in this program is fixed (and stored in quantumNum ) with

sigma = Z - Table[ Total[numEperOrb[[Range[1, i]]]] - ( numEperOrb[[i]] + 1), {i, 1, norbitals}];

In the rescaled variables the exponential dependence is Exp[-sigma x] and with this choice

sigma x= ((Z - N)/Z) Z r = (Z-N) r

This is the correct dependence for an hydrogen exponential with effective charge Z-N.  In general sigma is the effective charge for the hydrogen function.

The function InitEne generates also a parameter alpha which will be useful below in the solution of the HF equation.

In some problems (see an example at the end of this notebook) it can be useful to have a fixed grid for several configurations. In this case one can obtain the grid with

getgridFixed[maxrad, Z, h, x0]

This function is identical to parametrigrid but Z and the maximum radial number are fixed by the user.

Functions

parametrigrid[quantumNum_,hin_,x0in_]:=Module[
{x,x0,h,Z,nP,xmax,rho,rho2,rho32,grid},
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} = 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},
(* 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}];

These are the reduced wave functions for hydrogen

χ[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)] ;

Initialization of the functions

In the computation of the norm with integ we have to remember that there is factor ρ due to dx = ρ dξ  and a factor ρ due to the rescaling of the wave functions.

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

Examples

The grid

Let us consider again the Carbon and let us choose for LSlevel the ground state

degIoniz = 0;
electrons = ElementData[6,"ElectronConfiguration"];
quantumNum = getquantumNum[electrons,degIoniz];
listLS=choiceMult[];

LSlevel = {ground};

From this point on the function cycleHartreeFock will execute everything.

To understand how it works let us say a few things. There is a call to parametriMultipletto:

{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];
{Lmultipletto,Smultipletto}

After this step the function will use some predefined parameters (which can be changed with the options) to generate a grid of points.

The variable grid is global and  protected. The necessary parameters are x0 and h (there will be two options hgrid and x0grid). To show the mechanism we use a local variable gridLoc

h = 0.1; x0 = 0.1;
gridLoc = parametrigrid[quantumNum,h,x0];

The generated output has the format

{x,x0,h,rho,rho2} = gridLoc;

The total number of points is

Length[x]

while the explicit coordinates are

x

integ

Let us consider now an example for the function integ. Define a function fexample

fexample[z_]:= z(1.+1.7 z^2/Pi) Exp[-z];
NIntegrate[fexample[z],{z,0,∞}]

To use integ let us remember that we have the ξ variable and that dx = ρ dξ, then:

integ[rho*fexample[x],h]

InitPsi

We give an example for the initialization. With Carbon configuration

{E0,alpha}=InitEne[quantumNum];
psi = InitPsi[quantumNum,gridLoc];

The energies are those of Hydrogen

E0

and we have 3 wave functions, one for each orbital

Length[psi]

Let us make a graph of one of  them in the usual, non rescaled, variables

Zatom = quantumNum[[1]];
ListPlot[Transpose[{x/Zatom,(Sqrt[Zatom])(Sqrt[rho])psi[[2]]}],Joined-> True,
PlotRange-> {{0,8},{-1,0.5}}]

The potentials

Definitions

More details on this subject are in notebook [NB8].

Various potential can be reduced to the computation of

(Y^k)[a,b,r] = (∫_0^r)((((s/r)))^k)P[a,s]P[b,s]ds + (∫_r^∞)((((r/s)))^(k+1))P[a,s]P[b,s]ds

where P[a,r] are reduced radial wave functions for orbitals a,b.

The computation of these functions can be reduced to a solution of a couple of differential equations, or better to a couple of integrals

(d/dr)(Z^k)[a,b,r]= P[a,r]P[b,r]- (k/r)(Z^k)[a,b,r];
(d/dr)(Y^k)[a,b,r]=(1/r)( (k+1)(Y^k)[a,b,r] - (2k+1)(Z^k)[a,b,r] );
(Z^k)[a,b,0]=0;   (Y^k)[a,b,r] -> (Z^k)[a,b,r]    per    r -> ∞

or:

(d/dr)(r^k) (Z^k)[a,b,r]=(r^k) P[a,r]P[b,r];
(d/dr)(((Y^k)[a,b,r])/(r^(k+1)))= - ((2k+1)/(r^(k+2)))(Z^k)[a,b,r];

To obtain an accurate value for the integrals it is useful to know the small r behavior of the functions (Z^k). In this region the potential is well approximated by the coulombic potential of the nucleus and we can show that (see  [NB8]):

(Z^k)[a,b,r] ∼ (1/((Subscript[L, a])+(Subscript[L, b])+k+3))r P[a,r]P[b,r](1 + ((Z r)/((Subscript[L, a])+(Subscript[L, b])+k+4))((1/((Subscript[L, a])+1))+(1/((Subscript[L, b])+1)))).

The factor Z which appears in the rhs is the atomic number. Let us note that Z[0,a,b][∞] is 1 for a=b and 0 for a≠b (it is the normalization condition for radial wave functions).

We will use an exponential grid.

Let us note first of all that the form of the functions does not change, as the factor 1/Z in dr = dx/Z is reabsorbed in  the change of normalization of the P's:

r = (x/Z);    (Subscript[P, a])[r] = (Sqrt[Z]) (Subscript[Q, a])[x];

With

x = (Subscript[x, 0])*(Exp[ξ]-1);    ρ= (Subscript[x, 0])Exp[ξ]= x+(Subscript[x, 0]);

one obtains

(d/(dξ))(x^k) (Z^k)[a,b,x]=(ρ^2)(x^k) ψ[a,x]ψ[b,x];
(d/(dξ))(((Y^k)[a,b,x])/(x^(k+1)))= - ((2k+1)/(x^(k+2)))ρ (Z^k)[a,b,x];

In the small r behavior again the factor 1/Z produced by r = x/Z cancel with the normalization change in P. The term (Z r) in the expression becomes x and we can write

(Z^k)[a,b,x] ∼ c1 x Pa Pb(1 + c2 x)

In ξ variables, posing for shortness fH = (ρ^2) ψ[a,x]ψ[b,x] for  the k -independent term of the function, we have

(Z^k) ∼ c1 x ρ (Subscript[ψ, a]) (Subscript[ψ, b])(1 + c2 x) = c1 x fH/ρ (1 + c2 x) = c1 x fH/(x+x0) (1 + c2 x)

The numerical computation is performed on the exponential grid ξ.  We will use a modification of the Simpson's formula (see  [NB8]):

(∫_(Subscript[x, j-1])^(Subscript[x, j+1]))f[x]dx = (h/90)(-f[j+2] + 34 f[j+1] + 114 f[j] + 34 f[j-] - f[j-2])+O((h^7));

The computation is reduced to a solution of two recursive relations (see  [NB8]): the computation of Z starts from small r, that of Y from large r.

These functions are heavily used in the program and make very simple operations. It is convenient to use Mathematica Compile instruction (see Help) which produce a faster code.

The function

calcpotHartree[psi, LperOrb, x, x0, h]

has as only task to call the computation of functions Z and Y for k=0 and divide the result by r (i.e. x in the rescaled lengths). This function, the simplest one, is the only function used in Hartree equation.

In Hartree Fock equation the function Y and Z are put together to form the potentials for each orbital. These functions are named below calcpotYorb and calcpotXorb. Their meaning and their use is explained below.

Functions

(* 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},
(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 = 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},
{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},
{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]]

Hartree and Hartree-Fock energy

NOTE: In the following all the energies obtained as results are expressed in atomic units. The internal used quantities are rescaled energies and are multiplied by (Z^2)before the output. The only exception to this rule are intermediate results for single orbitals which, when requested, are printed in Ry. The circumstance will be explicitly signaled.

Hartree energy

Function:  calcEneHartree[psi,grid,quantumNum];

The energy in the Hartree model is written as

E = (∑_a)(Subscript[q, a])∫(Subscript[P, a])(-(1/2)(d/(d (r^2)))+(((Subscript[L, a])((Subscript[L, a])+1))/(2(r^2)))-(Z/r))(Subscript[P, a]) dr +
 (∑_a)(((Subscript[q, a])((Subscript[q, a]) - 1))/2)∫(1/r)Y[0,a,a][r] ((Subscript[P, a])^2) dr+ (∑_(b<a))(Subscript[q, a])(Subscript[q, a])∫(1/r)Y[0,b,b]((Subscript[P, a])^2)dr

The last two terms are due to electrostatic repulsion of the electrons. The first one is for electrons in the same orbital, the second one for electrons in different orbitals. From the expression of Y functions

(1/r)Y[0,a,a][r] = (1/r)(∫_0^r)P[a,s]P[a,s]ds + (1/r)(∫_r^∞)((r/s))P[a,s]P[a,s]ds

We can recognize the electrostatic potential at distance r generated by a charge density  ρ = ((Subscript[P, a])^2).  The first term is due to the internal charge (Gauss theorem) the second one to the external charge (let us note that the r dependence is only in the integration limit).

Let us now consider  exponential coordinates. Remembering all normalizations described above the kinetic part becomes

K = ∫(Subscript[ψ, a])(-(1/2)((d^2)/((dξ)^2))(Subscript[ψ, a])) dξ + ∫(Subscript[ψ, a])(((ρ^2)/(x^2))(((Subscript[L, a])((Subscript[L, a])+1))/2)+(1/8))(Subscript[ψ, a])dξ

In the computation we use in a reversed way the approximation for the second derivative with:

((d^2)/((dξ)^2))(Subscript[ψ, a]) ≃ (δ^2)(1 - (1/12)(δ^2))(Subscript[ψ, a])

In ξ variables the rescaled (with the factor 1/Z^2) potential energy takes the form

U = -∫(ρ^2)(Subscript[ψ, a])(1/x)(Subscript[ψ, a])dξ + (1/Z) (∑_a)(((Subscript[q, a])((Subscript[q, a]) - 1))/2)∫(1/x)Y[0,a,a]((Subscript[ψ, a])^2) dξ +(1/Z)(∑_(b<a))(Subscript[q, a])(Subscript[q, b])∫(1/x)Y[0,b,b]((Subscript[ψ, a])^2)dξ

The functions Y/x have been computed by calcpotHartree[], and the integrals are performed with the function integ[].

Everything is assembled in the function

calcEneHartree[psi , grid, quantumNum]

which computes the total energy, the kinetic energy and the virial, i.e. the ratio  U/K . The exact virial would be 2 for the coulombic interactions, then its value can give an idea of the error in the solution.

Hartree-Fock energy

Function:  calcEne[psi,grid,quantumNum,fk,gk,dfCoppia];

In the Hartree-Fock approximation the states are written as superposition of Slater determinants, organized to be eigenstates of L and S. In absence of additional degeneracies these states will be also eigenstates of the Hamiltonian (see text). This kind of coupling is called Russell-Saunders coupling and assume that spin interactions are negligible. Each level with a definite electronic configuration (as computed in Hartree approximation) splits in several terms with (L,S) fixed, each one with a degeneracy (2L+1)(2S+1).

The general form of the energy on these level is (see text)

E = (∑_a)(Subscript[q, a])I[a] + (∑_a)(((Subscript[q, a])((Subscript[q, a])-1))/2)(∑_(k=0))(Subscript[f, k])[a](F^k)[a,a] +(∑_(a<b))(Subscript[q, a])(Subscript[q, b])(∑_(k=0))((Subscript[f, k])[a,b](F^k)[a,b]+(Subscript[g, k])[a,b](G^k)[a,b]);

The first term embodies the kinetic energy and the electron nucleus interaction, and it is identical to the one of Hartree energy.

The last terms describe mutual electron interactions. These terms are linear combinations of radial integrals Fk and Gk. Fk are direct integrals while Gk are exchange integrals. The two forms come from two particle integrals of the form

∫drdr'|(Subscript[ψ, a])(r)((|)^2)(1/(|r - r'|))|(Subscript[ψ, b])(r')((|)^2); 
∫drdr'(Subscript[(OverBar[ψ]), a])(r)(Subscript[(OverBar[ψ]), b])(r)(1/(|r - r'|))(Subscript[ψ, a])(r)(Subscript[ψ, b])(r');

which can be translated in radial integrals using multipole expansion. With the notation

(U^k)(r,s) = (1/r)((((s/r)))^k) θ(r-s) + (1/s)((((r/s)))^k) θ(s-r) = (1/(Subscript[r, >]))(((((Subscript[r, <])/(Subscript[r, >]))))^k)

the radial integrals are:

(F^k)[a,b]=∫(Subscript[P, a])[r](Subscript[P, a])[r](U^k)(r,s) (Subscript[P, b])[s](Subscript[P, b])[s]drds;
 (G^k)[a,b]=∫(Subscript[P, a])[r](Subscript[P, b])[r](U^k)(r,s) (Subscript[P, a])[s](Subscript[P, b])[s]drds;

In these double integrals the first integration, let us say the one in the s-variable, gives functions like Y[k,a,b] already met. It will be computed with the function calcpotYk (almost identical to calcpotHartree ). The last integral is done using the function integ, already introduced.

The coefficients fk, gk depend on the multiplet. It is possible to define a mean energy by averaging eigenvalues of H, this means that each energy level is weighted with its multiplicity (2L+1)(2S+1). For this quantity the coefficient are easily computable (see text):

(Subscript[f, 0])[a]= (Subscript[f, 0])[a,b] =1;  (Subscript[f, k])[a]= - ((2(Subscript[L, a])+1)/(4(Subscript[L, a])+1))(((({{(Subscript[L, a]), k, (Subscript[L, a])}, {0, 0, 0}})))^2);
 (Subscript[g, k])[a,b]= - (1/2)(((({{(Subscript[L, a]), k, (Subscript[L, b])}, {0, 0, 0}})))^2);  (Subscript[f, k])[a,b] =0;  k > 1;

We used 3j-symbols to write the coefficients.

The variations from the mean are computable as they depend only on the uncompleted shells. A systematic way to obtain these deviations is to diagonalize the Hamiltonian for incomplete shells. This has been discussed in notebooks, [NB3, NB4]. We run these notebook for several configurations of incomplete shells and have collected the results in the form of a list, reported below. The list imposes some limitations on the possible calculations, due to the finite number of runs performed. Here we list the limitations and some possibilities:

The coefficients are stored in two vectors listadf and listadg with the format

listadf :  {{numElettroniEQ, La, Lb, Ltot, Stot}, δfk};  
 listadfdg: {{La, Lb, Ltot, Stot}, δfk, δgk}

Let us stress that the values stored refer to the total shift in the energy, while the coefficients given above refer to the energy for a single couple.

In this approach the variation in energy for N electrons gives also the variation for the case of N vacancies, i.e. 2(2L+1)-N electrons (see text).

Multiplets separation

Function: approxMultiplets[];

It can be useful to have an approximation of various possible LS terms using the wave functions computed with the mean energy, or the ground state. Clearly this is only an approximation with respect the full HF calculation as the wave functions are kept fixed.

The energies differ due to the different fk and gk coefficients, the radial integrals have been already computed at this level, then it is sufficient just to read the coefficients fk, gk and perform the sums.

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},

{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},

{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}];

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},

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

Examples

Let us see in some details how calcEne compute the energy by taking the irrealistic case of Carbon with hydrogen functions.

degIoniz = 0;
electrons = ElementData[6,"ElectronConfiguration"];
quantumNum = getquantumNum[electrons,degIoniz];
listLS=choiceMult[];

LSlevel = {ground};

In real computations at this point the function cycleHartreeFock compute the energy of LSlevel and the relative wave functions. For illustrative purposes let us generate some artificial functions.

{Lmultipletto,Smultipletto,fk,gk,dfCoppia}=parametriMultipletto[quantumNum,LSlevel,listadf,listadfdg];

If[Lmultipletto<0,nomeTermine = "Emed",
nomeTermine=""^ToString[2 Smultipletto +1]namesMult[[Lmultipletto+1]]];
Zatom = quantumNum[[1]];
Print["Z = ",Zatom,"; LSlevel: ",nomeTermine]

We chose standard parameter for the grid and use the global variable by unprotecting it.

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

This compute the Hartree energy :

{enetot,enekin,viriale}=calcEneHartree[psi,grid,quantumNum];

Print[namesEl[[Zatom]]," Z=",Zatom,";  ",nomeTermine,":  ",NumberForm[{enetot "au",enekin "au",viriale},{10,7}]]

Let us compute now the HF energy for the selected term and let us approximate the other terms (we repeat that the wave functions are irrealistic):

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

Print[namesEl[[Zatom]]," Z=",Zatom,";  ",nomeTermine,":  ",NumberForm[{enetot "au",enekin "au",viriale},{10,7}]]

These are energy differences between the terms:

diffE = approxMultiplets[];

The Hartree and Hartree-Fock equations

Hartree equation

The Hartree equation can be obtained by varying the previous energy functional

ℰ = (∑_a)(Subscript[q, a])∫(Subscript[P, a])(-(1/2)(d/(d (r^2)))+(((Subscript[L, a])((Subscript[L, a])+1))/(2(r^2)))-(Z/r))(Subscript[P, a]) dr +
                             (∑_a)(((Subscript[q, a])((Subscript[q, a]) - 1))/2)∫(1/r)Y[0,a,a][r] ((Subscript[P, a])^2) dr+ (∑_(b<a))(Subscript[q, a])(Subscript[q, b])∫(1/r)Y[0,b,b]((Subscript[P, a])^2)dr ;

and imposing norm 1 to the radial functions though a set of Lagange multipliers:

(δ/((Subscript[δP, a])[r]))( H - (∑_b) (Subscript[q, b])(Subscript[E, b]) (∫((Subscript[P, b])^2)dr - 1) ) = 0

With scaled coordinates x one easily get:

(-(1/2)(d/(d (x^2)))+(((Subscript[L, a])((Subscript[L, a])+1))/(2(x^2)))-(1/x))(Subscript[Q, a])+ (1/Z)(((Subscript[q, a]) - 1)(1/x)Y[0,a,a][x] + (∑_(b≠a))(Subscript[q, b])(1/x)Y[0,b,b][x])(Subscript[Q, a]) = (Subscript[∈, a])(Subscript[Q, a]) .

In exponential coordinates, naming V the electrostatic repulsion term

-(1/2)((d^2)/((dξ)^2))(Subscript[ψ, a])+(((ρ^2)/(x^2))(((Subscript[L, a])((Subscript[L, a])+1))/2)- ((ρ^2)/x)+(1/8))(Subscript[ψ, a]) + (ρ^2)(Subscript[V, a]) (Subscript[ψ, a]) = (ρ^2)(Subscript[∈, a])(Subscript[ψ, a]) .

This equation is easily discretized by previous methods, see also [NB9]:

-(1/(2(h^2)))(δ^2)ψ +(1 + (1/12)(δ^2))(((ρ^2)/(x^2))(((Subscript[L, a])((Subscript[L, a])+1))/2)- ((ρ^2)/x)+(1/8))(Subscript[ψ, a]) + (1 + (1/12)(δ^2)) (ρ^2)(Subscript[V, a]) (Subscript[ψ, a]) =(1 + (1/12)(δ^2)) (ρ^2)(Subscript[∈, a])(Subscript[ψ, a]).

The equation as the form of an ordinary generalized eigenvalue problem:

M ψ = ϵ S ψ

and can be solved with the function Eigenvalues of Mathematica.

In this notebook we will use the function:

loopHartree[maxIter, tolE, quantumNum, grid, psi, alpha, stampa]

Details on loopHartree

The function starts by using the input function psi to compute potential terms, after this step solves the eigenvalue problem and compute a new psi function. Then compute again the potentials. The loop continues and stops when the difference between new eigenvalues and old eigenvalues is smaller then tolE (a parameter fixed by the input) or when the number of iterations exceed maxIter.

The parameter alpha is used to stabilize the solution. The approach to a fixed point in an iterative process can be plagued by oscillations. The parameter alpha damp them by accepting as the new psi function a linear combination of old and new values

psi = (1-alpha)psi + alpha  psiOld

A quite conservative choice for alpha is done with the function Initalpha. The Boolean parameter stampa (True/False) is used to print the results at each iteration. In the solution the sign of the solution is taken to be positive for small x.

During a cycle the function compute the potential by a call to calcpotHartree. This function gives a list of potential, YHartree in the code. These potentials are multiplied by the occupation number of each orbital then the potential for the orbital to be studied is subtracted, in such a way to reconstruct the combination ((Subscript[q, a])-1). The vector vF used in the function store all electrostatic repulsion for each orbital.

A cycle on all orbitals is now performed, solving the eigenvalues problems. The states are ordered in growing order and the state we are interested in is selected by stato = numrad[[i]]. The radial quantum number is used because this is the number which distinguishes between wave functions for each value of angular momentum L.

The matrix Kin is the kinetic term, while the matrix matV is simply (1+(δ^2)/12) (ρ^2) V.

The same construction works for the other terms in the equation.

The only peculiarity in the code is the addition of a constant baseV to the potential, such that all eigenvalues become positive. This is not really used in the function but could be implemented by the reader, and could give better results for large atoms. We will use in a very limited form the function loopHartree and have not used the parameter baseV.

Hartree - Fock equation

The Hartree-Fock equation can be obtained by varying the previous energy functional

ℰ = (∑_a)(Subscript[q, a])I[a] + (∑_a)(((Subscript[q, a])((Subscript[q, a])-1))/2)(∑_(k=0))(Subscript[f, k])[a](F^k)[a,a] +(∑_(a<b))(Subscript[q, a])(Subscript[q, b])(∑_(k=0))((Subscript[f, k])[a,b](F^k)[a,b]+(Subscript[g, k])[a,b](G^k)[a,b]);

and imposing orthogonality on the wave functions with Lagrange multipliers. Let us note that here we impose explicitly the orthogonality on the wave functions, we do not implicitly assume it as in Hartree case.

We need two Lagrange multipliers, (Subscript[E, a]) for the normalization and (Subscript[E, ab]), with b ≠ a for orthogonality (see also text). The variational principle is written as

(δ/((Subscript[δP, a])[r]))( E - (∑_b) (Subscript[q, b])(Subscript[E, b]) (∫((Subscript[P, b])^2)dr - 1) - (∑_(b<c)) (Subscript[q, b])(Subscript[q, c])(Subscript[E, bc])∫(Subscript[P, b])(Subscript[P, c])dr) = 0

The resulting equation in x coordinates is

(-(1/2)(d/(d (x^2)))+(((Subscript[L, a])((Subscript[L, a])+1))/(2(x^2)))-(1/x))(Subscript[Q, a])+ (1/Z)(((Subscript[q, a]) - 1)(1/x)(∑_(k=0))(Subscript[f, k])[a]Y[k,a,a][x] + (∑_(b≠a))(Subscript[q, b])(1/x)(Subscript[f, k])[a,b]Y[k,b,b][x])(Subscript[Q, a]) +
(1/Z)(1/x)(∑_(b≠a))(∑_(k=0))(Subscript[g, k])[a,b]Y[k,a,b][x](Subscript[Q, b])[x] = (Subscript[∈, a])(Subscript[Q, a])[x] + (∑_(b≠a))(Subscript[q, b])(Subscript[∈, ab])(Subscript[Q, b])[x]

This is an integro-differential equation (like that of Hartree) has the unknown function (Subscript[Q, a]) is both in the derivative and in the potential term Y[k,a,b].

The terms with an explicit dependence on (Subscript[Q, a])are named direct potential terms, while the other are due to exchange potentials. Let us name 𝒴 the direct terms and 𝒳 the exchange terms. The equation takes the form

-(1/2)(d/(d (x^2)))(Subscript[Q, a]) + (Subscript[𝒴, a]) (Subscript[Q, a]) +  (Subscript[𝒳, a]) = (Subscript[∈, a])(Subscript[Q, a]) + (∑_(b≠a))(Subscript[q, b])(Subscript[∈, ab])(Subscript[Q, b])

This equation transforms in the ξ variable exactly in the way of Hartree equation.

In the iterative process the  𝒳  term is considered as an Inhomogeneous term, i.e. is computed with the current value of psi.

The main steps in solving for the orbital a are

L (Subscript[ψ, a]) + (Subscript[X, a]) = (Subscript[ϵ, a]) S (Subscript[ψ, a]) + dE

((((Subscript[ψ, a]),L (Subscript[ψ, a])) +( (Subscript[ψ, a]),(Subscript[X, a])))/(((Subscript[ψ, a]),S (Subscript[ψ, a])))) = (Subscript[ϵ, a])

ψ = LSolve[ L-ϵ S,- X+dE].

Some finer details

Nothing in this process implies that when we start with a wrong function the process converges toward the solution we are looking for instead of another eigenstate. To remedy instead of a direct call to LSolve the program make a call to a function calcolaPsi which itself call risolviSistemaPsi. This last function call LSolve and at the same time with a function contazeri computes the number of nodes of the function. As the function oscillates the function compute in effect the number of maxima and minima. A cut is used in such a way to accept only extrema which are 1/10 of the maximum, this is due to the wild (and small) oscillations at large distances (mainly due to numerical artifacts). The number of zeros must be equal to the radial number.

If this is true the function calcolaPsi accept the function otherwise change the value of the energy, around the initial value, to find a correct solution. It appears that at most within 20 iterations the function do its job.

The last value of ψ is returned to the principal cycle and is assumed as the current value of ψ, the same is true for the energy.

After a cycle on the orbitals the loop starts again up to a maximum of  maxIter cycles.

To the end of each cycle on the orbital the function are made orthogonal by using the Mathematica function Orthogonalize (in Mathematica 6, in previous versions the function GramSchmidt was used). Orbitals with different L are automatically orthogonal by their angular content, then the orthogonality must be imposed to functions with the same L. The vector  indOrto groups together orbitals with the same L. Take for instance a configuration sspspd. The vector of angular momentum would be LperOrb = {0, 0, 1, 0, 1, 2}. The following instruction make the necessary partition

LperOrb = {0,0,1,0,1,2};
indmom =  Table[Flatten[Position[LperOrb,k]],{k,0,Max[LperOrb]}]

indOrto = Select[indmom,Length[#]>0&]

We distinguish between indmom and indOrto vectors as there can be some missing orbitals, as in a configuration s-d.

Off diagonal coefficients

Let us now consider the problem of off diagonal Lagrange multipliers.

As explained in the text the HF equation admits in general orthogonal transformations which change the ψ and the Lagrange multipliers but leave invariant the minimum energy. Then the solution of the equation cannot be unique if one does not fix the free parameters. This problem is the analogous to the problem of phase-fixing in the usual Schrödinger equation. The problem is solved by Koopman theorem (see text) which states that when these rotations are admissible then one can use them to put (Subscript[ϵ, ab])= 0. At the same time the diagonal parameters (Subscript[ϵ, a])can be interpreted as the ionization energy, in the sense that if one imagine the ion built with the same wave functions as the full atom then (Subscript[ϵ, a])is exactly the ionization energy for that orbital.

In the text the problem has been discussed and the following cases have been classified

(Subscript[q, b])(Subscript[ϵ, ab])-(Subscript[q, a])(Subscript[ϵ, ba]) =(∑_(k=0))(∑_(c≠a,b))(Subscript[q, c])(R^k)[ac,bc]((Subscript[g, k])[a,c]-(Subscript[g, k])[b,c])+
                  (R^0)[bb,ab](1+(Subscript[q, b])(Subscript[g, 0])[a,b]) - (R^0)[aa,ab](1+(Subscript[q, a])(Subscript[g, 0])[a,b])+
      (∑_(k>0))(R^k)[aa,ab](((Subscript[q, a])-1)(Subscript[f, k])[a] - (Subscript[q, a])(Subscript[g, k])[a,b]) -(∑_(k>0))(R^k)[bb,ab](((Subscript[q, b])-1)(Subscript[f, k])[b] - (Subscript[q, b])(Subscript[g, k])[a,b]);

with

(R^k)[ab,cd]= ∫∫(Subscript[P, a])[s](Subscript[P, b])[s](U^k)[s,t] (Subscript[P, c])[t](Subscript[P, b])[t]dsdt;         (U^k)[s,t] = (((((Min[s,t])/(Max[s,t]))))^k);

There are now the following possibilities

In the text the following case are analyzed (for states with the same angular momentum)

((Subscript[q, b])-(Subscript[q, a]))(Subscript[ϵ, ab]) = (R^0)[bb,ab](1+(Subscript[q, b])(Subscript[g, 0])[a,b])- (∑_(k>0))(R^k)[bb,ab](((Subscript[q, b])-1)(Subscript[f, k])[b] - (Subscript[q, b])(Subscript[g, k])[a,b]);

2 (Subscript[q, a])(Subscript[ϵ, ab]) = ((Subscript[P, b]),(Subscript[ℱ, a]) (Subscript[P, a])) + ((Subscript[P, a]),(Subscript[ℱ, b]) (Subscript[P, b]));

All these cases are analyzed in the functions integralePerEpsilon, testEpsilon, EpsilonOffD. The last function gives as output the coefficients times the occupation number of the orbitals and the wave functions, in the form of a vector dE which can be inserted directly in the HF equation. It gives also the values (Subscript[q, b])(Subscript[ϵ, ab]). These parameters are energies and at the end of the loop they are multiplied by (Z^2) as all energies and put in the output.

Some more comments on this problem are added after the code of the function EpsilonOffD.

Warning: Mathematica 6

After each iteration we orthogonalize the eigenvectors. We found some problems in using the built-in function Orthogonalize in Mathematica 6 (a similar code was running on Mathematica 5). To avoid problems we write a pedestrian version of Gram-Schmidt orthogonalization procedure adapted to the scalar product we use (see function myProd etc. )

loopHartreeFock

All these steps are contained in a single function, loopHartreeFock, similar to the Hartree equation

loopHartreeFock[maxIter, tolE, quantumNum, grid, psi, alpha, fk, gk, dfCoppia, stampa, primoE_: 1]

All variables should be familiar to the reader at this point. As output the function gives

{ene, psi, eps, cycles}

i.e. the diagonal parameters (the energies) for each orbital, the wave functions, the off diagonal parameters times the occupation number and the number of cycles needed to find the result.

Let us note that with a single orbital there is not exchange term and the procedure must be a little modified. In this case the program transfer the control directly to the function loopHartree. This is the case for the ground state of Helium.

The loop stops when the convergence or maxIter are reached. loopHartrreFock can be used directly be the reader or by a driver (this is the usual case, see below). After this step the driver must make a choice on what to do. It can happen that a convergence is reached on a wrong function (i.e. with a wrong number of zeros). In this case the function gives as cycles the number maxIter and write an error message for the user.

The functions cycleHartreeFock and calcCore

To allow a simpler use of the programs we write a unique function cycleHartreeFock which do all the intermediate work and some checks on the convergence of the procedure. The call sequence is

cycleHartreeFock[LSlevel_List:{1},opts___]

The input is an optional list (LSlevel) and zero or more optional arguments.

If nothing is written the variable LSlevel take the value {1}, i.e. a string of length one and this trigger the computation of the ground state of the atom.

The second argument has 3 blanks, then it represent a variable number of arguments, even zero: this argument will contain the options we can fix. The default options are fixed, outside the function with

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

Let us give some details on the internal structure of cycleHartreeFock.

Warning: it can happen in some cases that for particular combinations of parameters the result converge toward a solution which is probably wrong. A signal of this situation is a big difference between eigenvalues and the corresponding eigenvalues in the Hartree equations. In this case the function print a Warning message and the reader must try the computation with different parameters.

This situation appears as an example in the excited configuration 4d for Li and for a particular combination of the parameters h and x0. In the usual cases (atoms or ions in low excitation levels) we have not found any strange behavior. In the case of Li with a change in the grid the fake solution disappears, confirming the nature of numerical instability.

Together with cycleHartreeFock we write the function calcCore, almost identical but with less options. This function simply disregard the peripheral electrons and solve the HF equation for the closed shells. Let us note that in the output all orbitals are present: the peripheral electrons are written as hydrogen functions.

We give a small summary of the options, reporting what has been shown at the beginning of the file

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 for 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,
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];

Comments on the function EpsilonOffD.

The function compute off diagonal elements of the matrix (Subscript[ϵ, ab]). Let us consider as an example the excited level of Mg with an electronic configuration

1 s^2 2 s^2 2 p^6 3 s 3 p

There are 5 complete orbitals. The angular momenta for closed and non closed shells are (pgusciNC list the non closed shells)

pLorb = {0,0,1,0,1}; pgusciNC = {4,5};

indmom make a partition of the orbitals using the angular momentum

pindmom=Table[Flatten[Position[pLorb,k]],{k,0,Max[pLorb]}]

In the function there are two loops. The external one, (with index iNC) runs on incomplete shells, while the internal one (variable irun) runs on the complete shells for all angular momenta except the last one. The instruction Which inside the loop determines how to compute eps. For Mg the variable iNC vary from 1 to 2 (number of incomplete shells). The variable ib takes the value gusciNC[iNC], then has the values 4 and 5 respectively. For each value of ib there is a value for the index ia (see the code): in the first case (electron s) ia = 1,2; in the second case ia = 3 (electron p).

Both incomplete shells could have the same angular momentum, take for instance He in the configuration 1s2s. In this case in the loop at the first passage ib = 1, but ia is absent. In the second passage ib=2 and ia =1. In a configuration 1s^2 2s 3s: in the first passage ib = 2, ia = 1; in the second ib = 3, ia = 1,2.

In the function code a couple of times one has to consider products like f[a,b] g[b]  for all b but not summed. As an instance as the radial function in HF have to be multiplied in each point by rho[x] one has to make the product psi[[k]][[x]]*rho[[x]]. This product can be done in two different ways. Let us take as an example a couple of functions pϕ and pω (analogous to psi and rho) on a grid of 5 points:

pϕ = Table[ϕ[i,j],{i,1,2},{j,1,5}];
pω = Table[ω[i],{i,1,5}];
pϕ[[1]]

Il first method is

Transpose[Transpose[pϕ]*pω]

The second (see the Help of Mathematica  for Inner ):

Inner[Times,pϕ,pω,List]

The second method is the one used in constructing the matrix λ[a,b] = ϵ[a,b] q[b].

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},

{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},
 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},
(* 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-> 30,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,
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},
(* 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];

Examples

Example 1

Let us define some conversion constant for au to Rydberg and from Rydberg to cm^-1:

au2Ry = 2;
Ry2cm = 109737.31568525;

We have to choose an element and compute the quantum numbers of the configuration

degIoniz = 0;
electrons = ElementData[16,"ElectronConfiguration"];
quantumNum = getquantumNum[electrons,degIoniz];
listLS=choiceMult[];

We choose the LSlevel we want to study

LSlevel = {ground};

The standard choice is now to use cycleHartreeFock without options

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

We can get more informations by estimating the energy differences between various terms

diffEne = approxMultiplets[];

Or we can print with more digits the energies of the orbitals

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

We can plot each orbital separately

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

or all together

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

or only some of them

Show[plotPsi[psi,0.2,{1,2}],PlotLabel-> "orbitals"]

Example 2

In this second example we redo the previous calculation but at the end we will look for a more accurate solution.

degIoniz = 0;
electrons = electrons = ElementData[16,"ElectronConfiguration"];
quantumNum = getquantumNum[electrons,degIoniz];
listLS=choiceMult[];

Choice of LSlevel:

LSlevel = {fond};

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

We choose different values for  h, x0 and tolE and we use the function just computed as input

{ene,psi,enetot}=cycleHartreeFock[LSlevel,psiInput-> psi,hgrid-> 1./50,errHF-> 10^-7];//Timing

The virial is much closer to 2.

Example 3

As third example let us make a computation with a fixed grid.

The purpose here is just to have functions of different configurations on the same grid to make simple graphs, but the procedure can be useful in other cases.

Let us suppose for instance that we want to compare the orbital 1s of Li and the same orbital for an excited configuration 1s^2 3s. For the second configuration by default a different grid is produced and to compare the results we could save the grid of the first computation and transfer the grid with updatePsi, but it is easier to work with the same grid in both cases.

The default grid in the first case is

degIoniz = 0;
electrons = {{"1s",2},{"2s",1}};
quantumNum = getquantumNum[electrons,degIoniz];
LSlevel = {ground};

{ene,psi2s,enetot} = cycleHartreeFock[LSlevel,printHF-> False];

Length[grid[[1]]]

Let us fix instead:

maxrad = 3; Zatom = 3;
hF = 0.1; x0F = 0.1;
gridF = getgridFixed[maxrad,Zatom,hF,x0F];

{ene,psi2s,enetot} = cycleHartreeFock[LSlevel,printHF-> False, gridOldHF-> gridF];

Length[grid[[1]]]

We can do the same computation for the second configuration

electrons = {{"1s",2},{"3s",1}};
quantumNum = getquantumNum[electrons,degIoniz];
{ene,psi3s,enetot} = cycleHartreeFock[LSlevel,printHF-> False,gridOldHF-> gridF];

Length[grid[[1]]]

Show/@{plotPsi[psi3s,0.5,{1}],plotPsi[psi2s,0.5,{1}]}

The two functions coincide, in agreement with the idea that core functions in alkali are almost independent on the peripheral electron.

Show[plotPsi[psi3s,0.5,{1}],plotPsi[psi2s,0.5,{1}]]