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

npart=2; (* Number of electrons *)
ndig0=16; (* Accuracy *)
ndig=16; (* Accuracy of FindRoot *)
zmin=1/100;
zmax=3;
zstep=1/100;
dz=1/300; (* variation of Z to bracket a stationary point *)

If[$VersionNumber==2.2,
   SetDirectory["c:\sergeev\purdue\math\boron"],
   SetDirectory["~/math/boron"] ];
if="guess"<>ToString[npart]<>".dat";
guess=ReadList[if,Release@Table[Number,{2 npart+6}]];
npoints=Dimensions[guess][[1]];
Do[
  g=Table[{guess[[np,1]],(guess[[np,2 n+3]]+guess[[np,2 n+4]] I)*
(guess[[np,1]]+guess[[np,2]] I)},{np,npoints}];
  guessa[n]=Interpolation[g],{n,npart}];
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,a1,a2,a,pri,x,Z];
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}];
enz=(meant-Z meanne+meanee)/mean;
en1z=-meanne/mean;
Clear[mean,meanne,meant,meanee];
"---------- Computation of the energy functional completed.">>>watchc;
Do[enzs[n]=D[enz,A[n]];en1zs[n]=D[en1z,A[n]];
   Do[enzs[n,l]=D[enzs[n],A[l]],{l,n}],{n,npart}];
"---------- Computation of the derivatives completed.">>>watchc;

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

(* Printing a table of E, E', and E'' to a file resultN.dat *)
Write[os,"Z",OutputForm[","],
 "Re E",OutputForm[","],"Im E",OutputForm[","],
 "Re E'",OutputForm[","],"Im E'",OutputForm[","],
 "Re E''",OutputForm[","],"Im E''",
 Sequence@@Flatten@Table[{OutputForm[","],ToString@OutputForm@Re@A[n]},{n,npart}] ];

A1=Table[Null,{n,npart}];
A2=Table[Null,{m,npart},{n,npart}];

Do[ (* Start cycle *)
(* Finding stationary points of the energy functional *)
   en=enz/.Z->z;
   Do[ens[n]=enzs[n]/.Z->z,{n,npart}];
   Do[a0[n]={guessa[n][z-dz],guessa[n][z+dz]},{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;
   en1=(en1z/.Z->z)/.s;
   Do[a[n]=A[n]/.s,{n,npart}];
   Do[A1[[n]]=(en1zs[n]/.Z->z)/.s,{n,npart}];
   Do[Do[A2[[l,n]]=A2[[n,l]]=(enzs[n,l]/.Z->z)/.s,{l,n}],{n,npart}];
   en2=-A1.MatrixPower[A2,-1].A1;
(* Printing results *)
   line={prire[z,4], pri[en0,6], pri[en1,6], pri[en2,6],Table[prire[a[n],5],{n,npart}]}//Flatten;
   Print[Sequence@@line];
   Write[ os, Sequence@@line];
,{z,zmin,zmax,zstep}]; (* End cycle *)
Close[os];
"---------- Completed.">>>watchc;