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

npart=6; (* Number of electrons *)
ndig0=16; (* Accuracy *)
ndig=16; (* Accuracy of FindRoot *)
zmax=10;
zmin=-2;
zstep0=1/5;
ca=1/5; (* If Hessian changes less than ca, then increase Zstep two times,
              else decrease it two times *)
zshift=1/100; (* Shift to imaginary direction to avoid singularities on the real axis *)
ashift=(1+I)/100; (* small initial perturbing *)

If[$VersionNumber==2.2,
   SetDirectory["c:\sergeev\purdue\math\boron"],
   SetDirectory["~/math/boron"] ];
os=OpenWrite["guess"<>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"];
nhalf=Floor[npart/2]; (* Number of bi-particle orbitals *)

(* 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]];
   Do[enlas[n,l]=D[enlas[n],A[l]],{l,n}],{n,npart}];
"---------- Computation of the derivatives completed.">>>watchc;

(* Initial guess at Z=Zmax *)
If[zmax!=10,
  Print["Warning: Initial guess at Z=Zmax may be inaccurate."]];
If[npart==2,(
  a0[1]=1.080;
  a0[2]=0.857;
 )];
If[npart==3,(
  a0[1]=1.027;
  a0[2]=0.882;
  a0[3]=0.228;
 )];
If[npart==4,(
  a0[1]=1.039;
  a0[2]=0.855;
  a0[3]=0.218+0.088 I;
  a0[4]=0.218-0.088 I;
 )];
If[npart==5,(
  a0[1]=1.046;
  a0[2]=0.850;
  a0[3]=0.207+0.082 I;
  a0[4]=0.207-0.082 I;
  a0[5]=a0[6]=0.386
 )];
If[npart==6,(
  a0[1]=1.051;
  a0[2]=0.846;
  a0[3]=0.198+0.076 I;
  a0[4]=0.198-0.076 I;
  a0[5]=0.366
 )];

Do[a0[n]={a0[n],a0[n]+(-1)^n ashift},{n,npart}];

z=N[zmax,ndig0];
zstep=N[zstep0,ndig0];

(* Printing operator for real numbers*)
(* 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} );

cond=True;
en2s=Table[Null,{n,npart},{l,npart}]; (* Matrix of second derivatives *)
d0=1;
While[cond, (* Start loop *)
(* Finding stationary points of the energy functional *)
   la=1/(z+zshift I);
   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)/la^2;
   Do[a[n]=A[n]/.s,{n,npart}];
   Do[en2s[[n,l]]=en2s[[l,n]]=
      (enlas[n,l]/.Append[s,lambda->la]),{n,npart},{l,n}];
   h=Det[en2s];
   d1=h/la^npart; (* scaling *)
   da=Abs[(d0-d1)/d1];
   d0=d1;
   If[da<ca,(cond=(z>=zmin);
         zstep=2 zstep;z=z-zstep;
         Do[a0[n]={a[n],a0[n][[1]]},{n,npart}];
      (* Printing results *)
         line=Flatten[{pri[1/la,8], pri[en0,8],Table[pri[a[n],6],{n,npart}],
               pri[d0,6]}];
         line1=Flatten[{pri[1/la,4], pri[en0,4],Table[pri[a[n],3],{n,npart}],
               pri[d0,3]}];
         Write[ os, Sequence@@line];
         Print[Sequence@@line1] ),
              (zstep=zstep/2;z=z+zstep) ];
]; (* End of loop *)
Close[os];
"---------- Completed.">>>watchc;
Exit[];