(* Calculations for arbitrary atom,
                       including p-orbitals *)
(* Step c - minimization of energy functional *)
(*              by equating partial derivatives to zero *)

npart=5; (* Number of electrons *)
ndig=16; (* Accuracy *)

SetDirectory["~/math/boron"];
os=OpenWrite["result"<>ToString[npart]<>".dat",PageWidth->Infinity];
"Minimization of energy functional for an atom with npart electrons">>watchc;
Definition[npart]>>>watchc;
Clear[cf,mat,l,n,enlas,ens,A,a0,a,pri,x];
Get["cf"<>ToString[npart]<>".dat"];
Get["mat"<>ToString[npart]<>".dat"];

(* Summation of contributions from different permutations *)
perm=Permutations[Range[npart]];
lperm=Length[perm];
{mean,meanne,meant,meanee}=
Sum[pl=perm[[l]];
   {l,pl}>>>watchc;
   cf@@pl mat@@pl, {l,lperm}];
Clear[cf,mat];

(* Expressing energy functional and its derivatives through variables a$n=A[n] *)
$ModuleNumber=1;
Do[A[n]=Module[{a},a],{n,npart}];
enla=(meant-meanne+lambda meanee)/mean;
Clear[mean,meanne,meant,meanee];
"---------- Computation of the energy functional completed.">>>watchc;
Do[enlas[n]=D[enla,A[n]],{n,npart}];
"---------- Computation of the derivatives completed.">>>watchc;

(* Initial guess for exponents *)

(* ! Boron (N=Z=5) *)
a0[1]=1.083;
a0[2]=0.769;
a0[3]=0.144+0.046 I;
a0[4]=0.144-0.046 I;
a0[5]=0.248;

Do[a0[n]={a0[n],a0[n]+0.01},{n,npart}];

(* Define range of charges Z *)
zrange=Range[5,4,-1/100];
(*zrange=Flatten[{
 Range[189/100,1/10,-1/100],
 Range[9/100,1/100,-1/1000]
 }];*)
mz=Length[zrange];

(* Printing operator for complex numbers*)
pri[x_,n_]:=(
   xp=SetPrecision[x,n];
   xp=Chop[xp,Abs[xp]/10^n];
   {"  \t", Re[xp], " ", Im[xp]} );

(* Finding stationary points of the energy functional *)

Do[z=zrange[[nz]];
   zp=SetPrecision[z,4]//OutputForm;
   la=1/z;
   en=enla/.lambda->la;
   Do[ens[n]=enlas[n]/.lambda->la,{n,npart}];
   s=FindRoot[Release[Table[ens[n]==0,{n,npart}]],
     Release[Sequence@@Table[{A[n],a0[n]},{n,npart}]],
     WorkingPrecision->ndig,MaxIterations->99];
   en0=en/.s;
   Do[a[n]=A[n]/.s;
      a0[n]={a[n],a0[n][[1]]},{n,npart}];
   line=OutputForm/@Flatten[{zp, pri[en0,12],Table[pri[a[n],6],{n,npart}]}];
   Print[Sequence@@line];
   Write[ os, Sequence@@line]
,{nz,mz}];
Close[os];
"---------- Completed.">>>watchc;