Guide to Hartree approximation

Quantum Mechanics
                     K. Konishi, G. Paffuti

In this notebook we present in some details the numerical procedures used in the solution of the Hartree 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];

These are the global variables in this notebook.

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

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

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

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

Use of the notebook

What can be done and how

This notebook can be used for the calculatioin of the ground-state energy and energy of a few first excited levels of an atom in the Hartree approximation. The energies and wave functions are determined. The fine structure (i.e. spin interaction) is neglected in this notebook.
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 requires 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

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

Functions for the basic usage

{ene,psi,enetot}=cycleHartree[];

This single instruction computes the energies of all orbitals, the wave functions and the total energy.

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 calcCoreHartree which compute the same quantities only for the closed shells.

plotPsi[psi,zoom];

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

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

The zoom parameter allow to enlarge or reduce the graph.

Approximate multiplets

One can use Hartree wave function to compute an approximate value of the energy of a given multiplet, i.e. with given (L,S) quantum numbers.

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

elenco=choiceMultiplet[];

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

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

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

The instruction

eneHF=calcEneHFconHartree[psi,LSlevel];

or

eneHF=calcEneHFconHartree[psi];

Compute the energy of LSlevel ( or of the default value). There are some limitations on the computable terms, the reader can find more on the subject in subsequent notebooks.

The instruction

diffE=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.

Options

The function cycleHartree has several options which can be used with the usual Mathematica mechanism  option → val . These are the default values

Options[cycleHartree]

Their meaning is the following

For calcCoreHartree there is a similar set of options

Options[calcCoreHartree]

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

Summary of the theory

The Hartree equation is based on the intuitive description of an electron motion in an effective field generated by the rest of the electrons. It is expected an equation like

-(1/2)(∇^2)(Subscript[ψ, a]) - (Z/r)(Subscript[ψ, a]) + V(r) (Subscript[ψ, a]) = (Subscript[E, a]) (Subscript[ψ, a]);   V(r) = (∑_(b≠a))∫dr'((|(Subscript[ψ, b])(r')((|)^2))/(|r - r'|));

a is the orbital under study. The integral in the definition of V[r] is averaged over directions to give a central field. In this potential there will be a series of levels 1s, 2s, 2p, ... and ψ are the associated eigenfunctions. Equation () in a central field reduces to an equation for the reduced radial wave function. The simplest deduction of the equation is with the variational principle. The trial function is taken as a factorized (and not anti-symmetric)

Ψ((Subscript[r, 1]),(Subscript[r, 2]),...(Subscript[r, N]))= (Subscript[ψ, (Subscript[a, 1])])((Subscript[r, 1])) (Subscript[ψ, (Subscript[a, 2])])((Subscript[r, 2]))...(Subscript[ψ, (Subscript[a, N])])((Subscript[r, N]));
(Subscript[ψ, nlm])(r)= (1/r)(Subscript[P, nl])(r) (Subscript[Y, lm])(θ,φ);         (Subscript[P, nl]) = reduced radial functions

The variational principle imply

(δ/(Subscript[δψ, c]))(∫Ψ(-(1/2)(∑_(a=1)^N)((Subscript[∇, a])^2) - (∑_(a=1)^N)(Z/(Subscript[r, a]))+ (∑_(a< b)^N)(1/(|(Subscript[r, a]) - (Subscript[r, b])|)))Ψ  - (∑_(a=1)^N)(Subscript[E, a])(∫((Subscript[ψ, a])^2)-1)) = 0;

The parameters E appear as Lagrange multipliers to impose the normalization condition.

A multipole expansion of the coulombic term gives (see text):

-(1/2)(Subscript[P, a]^″) + (((L(L+1))/(2(r^2)))-(Z/r))(Subscript[P, a])+((Subscript[q, a]) - 1)(1/r)Y[0,a,a][r] (Subscript[P, a]) + (∑_(b≠a))(1/r)Y[0,b,b](Subscript[P, a]) = (Subscript[E, a])(Subscript[P, a]);

The functions Y[0,b,b]/r describe the electrostatic repulsion (the notation is borrowed from the analogous Hartree Fock quantities). Their explicit form is

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

This expression gives the potential at point r due to the sum of the charge distributed inside the sphere of radius r and external to this sphere.

In the following we will describe in some details the following points

The techniques needed to compute the Y functions are described in detail in notebook [NB8]

Electronic configuration

The input

The notebook accepts two different input 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]

Grid

The method

We start by rescaling the variables with  Z,  this allows 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 unit the Hartree equation takes the form

-(1/2)(Subscript[Q, a])'' + (((L(L+1))/(2(x^2)))-(1/x))(Subscript[Q, a])+(1/Z)((Subscript[q, a]) - 1)(1/x)Y[0,a,a][x] (Subscript[Q, a]) +(1/Z) (∑_(b≠a))(1/x)Y[0,b,b][x](Subscript[Q, a]) = (Subscript[ϵ, a])(Subscript[Q, a]);  
 (Subscript[E, a]) = (Z^2)(Subscript[ϵ, a])

Let us note that Z appear as a denominator in the interaction energy. Energies in atomic units are (Z^2)times the eigenvalues.

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 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[ξ, j]) = j h;  j = 1, ... N;   
(Subscript[Q, a])[x] = (Sqrt[ρ])(Subscript[ψ, a])[ξ];                      ρ = x0 Exp[ξ];

All points are internal points, for ξ = 0 and in the next to last point (j=N) we impose ψ = 0. 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 appear often in the formulas.

This change of variables transforms the equation in

-(1/2)((d^2)/((dξ)^2))(Subscript[ψ, a]) + (-((ρ^2)/x)+((ρ^2)/(x^2))(((Subscript[L, a])((Subscript[L, a])+1))/2)+(1/8))(Subscript[ψ, a]) + (ρ^2)(1/x)Y[0,a,a](Subscript[ψ, a]) = (ρ^2) (Subscript[ϵ, a]) (Subscript[ψ, a]);

x and ρ are expressed through ξ. 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 scalar product is written as

∫(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 ψ.

The original wave functions are related to ψ by

(Subscript[P, a]) = (Sqrt[Z]) (Sqrt[ρ])(Subscript[ψ, a])

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.

It is convenient to collect all elements of a grid in a single structure, that will be called grid:

grid = {x, x0, h, rho, rho2};       (here    rho2 = rho*rho );

We write a function  parametrigrid[quantumNum,h_,x0_]  to generate the grid. It requires as input the quantum numbers, x0 and h. 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.

Code

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

The potentials

Computational method

More details on this subject are in notebook [NB8]. The computation of the potential functions can be reduced to a solution of a couple of differential equations, or better to a couple of integrals. We will solve a slightly more general problem as these more general functions will be needed in the Hartree-Fock case and for the calculation of the energy. We write r instead of x for clarity

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

The integrals are

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

The more general functions are

(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

and the corresponding relations are

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

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 of Z - functions starts from small r-values, while the computation of Y starts from large r. The discretized integral can be defined by a recursion relation. To write this relation for Z-functions it is necessary to know their values at small r. 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).

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.

These functions are coded below with funZ and funY. Apart from these below we give two more functions. The first one, integ perform the numerical integral in ξ for a function F[ξ] vanishing at the origin (in the cote the value F[0] is omitted). The last function, in  the form

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

just call the functions for Z and Y with k=0 and generate potential terms.

Code

(* 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.)))];
(* qui va 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]];

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

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

Initialization

Definitions

To start the iteration procedure we need a starting expression for the wave functions. We will use the hydrogen wave functions computed with a screening parameter sigma.  To simplify some steps we will initialize also the energy eigenvalues and some matrices that will be used below.

Let us note that in the computation of energies we use the expression 1/n^2  (i.e. with Z=1) as we are using the rescaled energies.

In the function InitEne   the vector alpha is computed which will be used to stabilize the iteration procedure.

The function InitPsi defines psi as stated above. The factor 1/(Sqrt[Z]) is for the normalization (it could be omitted as we will normalize the functions below).

The scale of the functions is (sigma x) . If we call N the screening electrons and use the formula given above for sigma (which had a factor 1/Z):

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

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

We will use the function Chop as the input for functions using Compile must be in machine-precision numbers.

Code

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

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

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

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

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

Equations, solutions and energies

The iterative procedure

The equations are

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

For the discretization we shall 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])) = (1 + (1/12)(δ^2)) (ρ^2)(Subscript[ϵ, a])(Subscript[ψ, a])

The operator (1+ (δ^2)/12) acting on a vector V has non zero matrix elements only on the three main diagonals. Its effect,(Subscript[(1 + (1/12)(δ^2)), ij])(Subscript[V, j])(Subscript[ψ, j]), can be written as a tridiagonal matrix (Subscript[V, ij]) built with the vector V, more details and examples are in notebook  [NB7]. After this transformations we have to solve the eigenvalue equations

-(1/(2(h^2)))(Subscript[((δ^2)), ij])(Subscript[ψ, j]) +(Subscript[vMat, ij]) (Subscript[ψ, j])= ϵ  (Subscript[matS, ij]) (Subscript[ψ, j]);

matS is the tridiagonal matrix constructed with rho2, vMat the matrix with the potential (and the constant 1/8). In the function coded below the first term is indicated with Kin.

This equation can be solved with the Mathematica function Eigenvalues.

The function Eigenvalues can be used in two ways: one can either compute all eigenvalues or only a part of them, by selecting an option like Method -> Arnoldii (see Help). The second case is easier to compute if all eigenvalues are positive. To have only positive energies we add a constant to the potential, baseV = 1/x[[1]], i.e. the coulombic nucleus potential changed in sign. This constant is subtracted from the computed eigenvalues. This method will not be used below but we leave the computation of the constant in the function for possible use by the reader.

The matrix matS is the one computed in the function Initmatrici, as well as the matrix Kin and a vector filled with 1/8, called vec8.

After these preliminaries the iterative procedure can be described:

psi = (1 - alpha) psi + alpha psiold

All this is realized in the function loopHartree below. The internal variable stampa controls the printing in each cycle (True or False). The results are printed in Rydberg, we remember that 1 au = 2 Ry.

Energy in the Hartree model

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.

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 done in the function calcEneHartree[psi ,grid, quantumNum] which gives 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.

Functions cycleHartree e calcCoreHartree

To simplify the use of the program we give a unique function, cycleHartree,  which performs all intermediate steps. The template is

cycleHartree[opts___]

The only argument, opts, has 3 blanks, i.e. it stands for a variable number of arguments, these will be the options that we can change.

The default options are fixed, outside the function code, with

Options[cycleHartree] = {cyclesHF-> 20,errHF-> (10.^(-3)),hgrid-> 0.1,x0grid-> 0.1,psiInput-> {},alphaHF-> {},firstEHF-> 1, printHF -> True,gridOldHF-> {}};

Here we briefly explain the various part of the code

Together with cycleHartree we write another function calcCoreHartree almost identical but with less options. This function simply disregard the peripheral electrons and pass to loopHartree only the core. It has to be noted that in the output the whole psi is provided, in such a way that it can be used as input for cycleHartree. The missing orbitals are replaced by hydrogen functions.

For readability we report here the summary presented at the beginning of the notebook for the options

Options

The function cycleHartree has several options which can be used with the usual Mathematica mechanism  option → val . These are the default values

Options[cycleHartree]

Their meaning is the following

For calcCoreHartree there is a similar set of options

Options[calcCoreHartree]

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

Code

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

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

Options[cycleHartree] = 
{cyclesHF-> 20,errHF-> (10.^(-3)),hgrid-> 0.1,x0grid-> 0.1,psiInput-> {},alphaHF-> {},firstEHF-> 1, printHF -> True,gridOldHF-> {}};

cycleHartree[opts___]:=
Module[{ene,psi,eps,cicli,
Lmultipletto,Smultipletto,fk,gk,dfCoppia,
nomeTermine,Z,enetot,enekin,viriale,
alpha,au2Ry,quantumNum,
maxIter,tolE,h,x0,psiIn,primoE,stampa,
gridOld},

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


quantumNum = numQ$HF;

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 
unless the same psiIn is requested more times, in this
case one must provide gridOld to the function *)

If[Length[psiIn] ==  0,
Unprotect[grid];
grid = parametrigrid[quantumNum,h,x0];
Protect[grid];
psi = InitPsi[quantumNum,grid];
,
If[Length[gridOld]== 0,gridOld = grid;];
psi= updategridPsi[x0,h,quantumNum,gridOld,psiIn];
];

au2Ry = 2.0;

{ene,psi,cicli}=loopHartree[maxIter,tolE,quantumNum,grid,psi,alpha,stampa];

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

  Z = quantumNum[[1]];

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

Print[namesEl[[Z]]," Z=",Z,";  ",
"E_H = ",NumberForm[enetot "au",{10,7}],"; Kin = ",
NumberForm[enekin "au",{10,7}],"; vir = ",
NumberForm[viriale,{10,7}],"; "]; 
Return[{ene,psi,enetot}];
];

calcCoreHartree[opts___]:=
Module[{ene,psi,psiTot,cicli,eps,ene1,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;

(* Questa grid è quella dell'atomo completo *)

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

{ene1,psi,cicli}=loopHartree[maxIter,tolE,nQCore,grid,psi,beta,stampa];

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

psiTot[[gusciC]] = psi;
ene[[gusciC]] = ene1;
enetot = 0.;
If[siEneHF,
{enetot,enekin,viriale}= calcEneHartree[psi,grid,nQCore];
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[calcCoreHartree] = 
{cyclesHF-> 20,errHF-> (10.^(-5)),hgrid-> 0.1,x0grid-> 0.1,alphaHF-> {}, printHF -> False,si$EneHF-> False};

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

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

Graphics

plotPsi[psi_,zoom_,statiIn_:{}]:=Module[{Z,numPrinc,sigma,norbitals,lplot,rho,j,quantumNum},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,((numPrinc[[j]])/(zoom sigma[[j]]))},All},PlotStyle->Hue[1-((0.9` istat)/norbitals)],Joined->True,PlotLabel->"orbital  "<>electronsnames[[j]]],{istat,1,Length[lplot]}]];lplot];

List of names

We define for an easy use a list with the names of elements and another one with the standard names of multiplets.

Unprotect[namesEl,namesMult];
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"};

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

Example

Let us define the electronic configuration with the Mathematica  function ElementData:

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

Here we compute the quantum numbers

quantumNum = getquantumNum[electrons,degIoniz]

and print some information

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

This is the default computation with Hartree equation

{ene,psi,enetot}=cycleHartree[];

Hartree-Fock energy

Analytic form of the energy

Functions computed in Hartree approximation can be used to compute the energy in the Hartree-Fock approximation, i.e. Slater determinants are written in terms of Hartree functions. In this way it is possible to have an approximation for the energy of the different spectral terms ( (L,S) levels).

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 Russel-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[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).

It is also simple, for each configuration, to write the term corresponding to Hund rule, i.e. with the maximum value for S and with the maximum value of L compatible with it. If n is the effective number of electrons (electrons or vacancies depending if the shell is less or more than half filled) then S = n/2 and L can be computed by writing Lz= L,,L-1,...  (see text).

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

elenco=choiceMultiplet[];

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

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

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

The instruction

eneHF=calcEneHFconHartree[psi,LSlevel];

or

eneHF=calcEneHFconHartree[psi];

compute the enrgy of LSlevel ( or of the default value).

The function

diffE=approxMultiplets[];

compute and print energies differences between computable energies.

NOTE: For more details on the computations of HF energies see the notebook on HF equation.

Code

Lists

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

Functions for the electronic configuration

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

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

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

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

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

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

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

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

parametriMultipletto[quantumNum_,LSlevelIn_,listadf_,listadfdg_]:=
Module[{gk,fk,la,lb,df,dg,dfTot,gusciNC,orbEQ,momL,numEeff,fatt,orbL,orbS,
Linput,Sinput,LMult,SMult,dfCoppia,
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}];];

(* Il LSlevel è ora di lungh. 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}];]

choiceMultiplet[]:=
Module[{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC,
Smult,Lmult,orb,lista,listanomi,elenco,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:   mean E: {}","  ground.: {ground} "];
elenco = {};
listanomi = {};
nomefond = "";
Which[
Length[gusciNC] == 0,

nomefond = ""^ToString[1]nomi[[1]];
listanomi = ""^ToString[1]nomi[[1]];
elenco = {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]}];
elenco = 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]}];
elenco = Table[{mult[[j]][[1,3]],mult[[j]][[1,4]]},{j,1,Length[mult]}];
];
Print["or {L,S} or elenco[[k]] from: ",listanomi," groundstate = ",nomefond];
Return[elenco]];

Energy

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}];
(* Integrali 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;
(* qui vengono definite le variabili globali *)
Unprotect[virialeLSlevel,FkHF,GkHF];
virialeLSlevel = viriale;
FkHF = (Z^2)*Fk; GkHF = (Z^2)*Gk;
Protect[virialeLSlevel,FkHF,GkHF];

{enetot,enekin,viriale}];

calcEneHFconHartree[psi_,LSlevel_:{fond}]:=
Module[{Lmultipletto,Smultipletto,fk,gk,dfCoppia,
nomeTermine,Z,enetot,enekin,viriale,
quantumNum},

quantumNum = numQ$HF;
Z = quantumNum[[1]];

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

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

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

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

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

{Z,norbitals,numEperOrb,LperOrb,numPrinc,numrad,sigma,gusciNC} = quantumNum;
If[Length[gusciNC]==1 ,
orbEQ = gusciNC[[1]];
nEQ = numEperOrb[[orbEQ]];
nEQeff = Min[nEQ,2(2LeQ+1)-nEQ];
 LeQ =LperOrb[[orbEQ]];
,nEQ=0;];
eneMult = {};
Which[
Length[gusciNC]==1 && Min[nEQ,2(2 LeQ+1)-nEQ]>1,
lista = {nEQeff,LeQ};
mult=Select[listadf,#[[1]][[{1,2}]] == lista&];
(* Warning  in listadf is stored the total energy, not the energy per couple *)
ene=Table[mult[[j,2]].Fk[[orbEQ,orbEQ,All]],{j,1,Length[mult]}];
eneMult=Table[{(ene[[i]]-Min[ene]),mult[[i,1]][[4]],
mult[[i,1]][[5]]},{i,1,Length[mult]}];
eneMult= Sort[eneMult,#1[[1]] < #2[[1]] &];
];
listaMult = Table[
""^ToString[2 eneMult[[j,3]] +1]namesMult[[eneMult[[j,2]]+1]],{j,1,Length[eneMult]}];

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

Examples

Example 1

We give again the computation for Carbon

au2Ry = 2;
Ry2cm = 109737.31568525;

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

quantumNum = getquantumNum[electrons,degIoniz];

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

{ene,psi,enetot}=cycleHartree[];

Plots

psigraf = psi;

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

All the orbitals plotted for a wider range of radial variable

Show[plotPsi[psigraf,0.1],PlotLabel->"orbitals"]

Spectral terms

We use the Hartree functions to compute the (L,S) levels.  Let us choose an LS level:

elenco = choiceMultiplet[];

LSlevel = {ground};

The ground LSlevel is the default then in effect previous instruction could be omitted.

The Hartree - Fock energy is

eneHF = calcEneHFconHartree[psi,LSlevel];

or (using the default value for LSlevel)

eneHF = calcEneHFconHartree[psi];

We can compute differences between spectral terms:

Warning: Before computing the differences one has to call at least once calcEneHFconHartree

diffE=approxMultiplets[];

Example 2

In this example we perform the same calculation but we show after this how the result can be improved by choosing a finer grid. We do not print intermediate results.

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

quantumNum = getquantumNum[electrons,degIoniz];

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

Note the choice of the options for the print

{ene,psi,enetot}=cycleHartree[printHF -> False];

Let us select now a smaller value for h and use the variable psiInput in the options to pass to the function the computed value of psi. The function automatically adapt the function to the new grid.

{ene,psi,enetot}=cycleHartree[errHF-> (10.^(-4)),hgrid-> 1./20,psiInput-> psi];

We see that the virial is closer to its expected value 2.