(* Calculations for arbitrary atom,
   including p-orbitals *)
(* Step b - calculation of matrix elements *)

npart=5; (* Number of electrons *)

Off[$ModuleNumber::modnc,General::spell1];
SetDirectory["~/math/boron"];
Clear[f,l,m,n,r,theta,phi,fs,fp,x,la,ma,na,lb,mb,nb,int,A,a,b,intne,int0,
intt,int0ee,intee,l1a,m1a,n1a,l2a,m2a,n2a,l1b,m1b,n1b,l2b,m2b,n2b,a1,a2,b1,b2,
zero,qa,qb];
(* int0[n,a]=Integral[Exp[-a r] r^n dr] *)
int0[n_,a_]=(n+2)!/a^(n+3);
(* int[la,ma,na,lb,mb,nb]=<la ma na|lb mb nb> *)
int[la_,ma_,na_,lb_,mb_,nb_]:=0;
int[l_,m_,na_,l_,m_,nb_]=int0[2 l,A[na]+A[nb]];
(* intne[la,ma,na,lb,mb,nb]=<la ma na|r^-1|lb mb nb> *)
intne[la_,ma_,na_,lb_,mb_,nb_]:=0;
intne[l_,m_,na_,l_,m_,nb_]=int0[2 l-1,A[na]+A[nb]];
(* intt[la,ma,na,lb,mb,nb]=<la ma na|T|lb mb nb>
   where T=L^2/(2 r^2) -1/(2 r^2) d/dr(r^2 d/dr) *)
(*
psia=r^l Exp[-A[na] r];psib=r^l Exp[-A[nb] r];
Unprotect[Re];Re[x_]:=1;
i=Integrate[r^2 psia (-1/(2 r^2) D[r^2 D[psib,r],r]),{r,0,Infinity}];
Simplify[(l(l+1)/2 int0[2 l-2,A[na]+A[nb]]+
Integrate[r^2 psia (-1/(2 r^2) D[r^2 D[psib,r],r]),{r,0,Infinity}])/.{
      Gamma[2 l]->(2 l)!/(2 l),Gamma[2(l+1)]->(2 l)! (2 l+1)}]
*)
intt[la_,ma_,na_,lb_,mb_,nb_]:=0;
intt[l_,m_,na_,l_,m_,nb_]=
   (1 + 3*l + 2*l^2)*A[na]*A[nb]*(A[na] + A[nb])^(-3 - 2*l)*(2*l)!;
(* int0ee[l1a,m1a,l2a,m2a,l1b,m1b,l2b,m2b,a1,a2]=
       <l1a m1a n1a, l2a m2a n2a|r12^-1|l1b m1b n1b, l2b m2b n2b>
   where A[n1a]+A[n1b]=a1, A[n2a]+A[n2b]=a2 *)
(*int0ee[l1a_,m1a_,l2a_,m2a_,l1b_,m1b_,l2b_,m2b_,a1_,a2_]:=0;*)
(* (2,2) *)
int0ee[0,0,0,0, 0,0,0,0, a1_,a2_]=2*(a1^2 + 3*a1*a2 + a2^2)/(a1^2*a2^2*(a1 + a2)^3);
(* 0 *)
int0ee[0,0,0,0, 1,1,0,0, a1_,a2_]=0;
int0ee[0,0,0,0, 0,0,1,1, a1_,a2_]=0;
int0ee[1,1,0,0, 0,0,0,0, a1_,a2_]=0;
int0ee[0,0,1,1, 0,0,0,0, a1_,a2_]=0;
int0ee[0,0,0,0, 1,0,0,0, a1_,a2_]=0;
int0ee[0,0,0,0, 0,0,1,0, a1_,a2_]=0;
int0ee[1,0,0,0, 0,0,0,0, a1_,a2_]=0;
int0ee[0,0,1,0, 0,0,0,0, a1_,a2_]=0;
(* 3*(10,10) *)
int0ee[1,0,0,0, 1,0,0,0, a1_,a2_]=(12*(2*a1^4 + 10*a1^3*a2 + 10*a1^2*a2^2 +
 5*a1*a2^3 + a2^4))/(a1^4*a2^2*(a1 + a2)^5);
int0ee[0,0,1,0, 0,0,1,0, a1_,a2_]=int0ee[1,0,0,0, 1,0,0,0, a2,a1];
(* 3*(10,17) *)
int0ee[1,0,0,0, 0,0,1,0, a1_,a2_]=(8*(a1^2 + 5*a1*a2 + a2^2))/(a1^2*a2^2*(a1 + a2)^5);
int0ee[0,0,1,0, 1,0,0,0, a1_,a2_]=int0ee[1,0,0,0, 0,0,1,0, a1,a2];
(* 3/2*((8,8)-(10,10)) *)
int0ee[1,1,0,0, 1,1,0,0, a1_,a2_]=int0ee[1,0,0,0, 1,0,0,0, a1,a2];
int0ee[0,0,1,1, 0,0,1,1, a1_,a2_]=int0ee[1,0,0,0, 1,0,0,0, a2,a1];
(* 3/2*(20,3) *)
int0ee[1,1,0,0, 0,0,1,1, a1_,a2_]=int0ee[1,0,0,0, 0,0,1,0, a1,a2];
int0ee[0,0,1,1, 1,1,0,0, a1_,a2_]=int0ee[1,0,0,0, 0,0,1,0, a1,a2];
(* 9/2*((22,22)-(24,24)) *)
int0ee[1,0,1,1, 1,0,1,1, a1_,a2_]=1/(5*a1^4*a2^4*(a1 + a2)^7)*
  144*(5*a1^6 + 35*a1^5*a2 + 103*a1^4*a2^2 + 161*a1^3*a2^3 + 
     103*a1^2*a2^4 + 35*a1*a2^5 + 5*a2^6);
int0ee[1,1,1,0, 1,1,1,0, a1_,a2_]=int0ee[1,0,1,1, 1,0,1,1, a1,a2];
(* 9/2*(21,24) *)
int0ee[1,0,1,1, 1,1,1,0, a1_,a2_]=(432*(a1^2 + 7*a1*a2 + a2^2))/
   (5*a1^2*a2^2*(a1 + a2)^7);
int0ee[1,1,1,0, 1,0,1,1, a1_,a2_]=int0ee[1,0,1,1, 1,1,1,0, a1,a2];
intee[l1a_,m1a_,n1a_,l2a_,m2a_,n2a_,l1b_,m1b_,n1b_,l2b_,m2b_,n2b_]=
   int0ee[l1a,m1a,l2a,m2a,l1b,m1b,l2b,m2b,A[n1a]+A[n1b],A[n2a]+A[n2b]];

matr[qa_,qb_]:=(ndif=0;neq=0;Clear[poseq,posdif];
   Do[If[Take[qa[[n]],2]==Take[qb[[n]],2],
         (neq++;poseq[neq]=n),(ndif++;posdif[ndif]=n)];
      If[ndif>2,Break[]],{n,npart}];
   result=Table[0,{4}]; (* <1>,<1/r>,<T>,<1/rij> *)
   If[ndif>2,Return[result]];
   If[ndif==2,(n1dif=posdif[1];n2dif=posdif[2];
     result[[4]]=intee[
        Sequence@@Flatten[{qa[[n1dif]],qa[[n2dif]],qb[[n1dif]],qb[[n2dif]]}] ]*
      Product[neq=poseq[n];int[ Sequence@@Flatten[{qa[[neq]],qb[[neq]]}] ],{n,npart-2}];
     Return[result])];
   If[ndif==1,(
     matr0=Product[neq=poseq[n]; int[ Sequence@@Flatten[{qa[[neq]],qb[[neq]]}] ],{n,npart-1}];
     ndif=posdif[1];
     result[[4]]=(matr0 Sum[neq=poseq[n];
         intee[ Sequence@@Flatten[{qa[[ndif]],qa[[neq]],qb[[ndif]],qb[[neq]]}] ]/
           int[ Sequence@@Flatten[{qa[[neq]],qb[[neq]]}] ],{n,npart-1}])//Cancel;
     Return[result])];
   If[ndif==0,(
     result[[1]]=matr0=Product[int[ Sequence@@Flatten[{qa[[n]],qb[[n]]}] ],{n,npart}];
     result[[2]]=matr0 Sum[intne[ Sequence@@Flatten[{qa[[n]],qb[[n]]}] ]/
                  int[ Sequence@@Flatten[{qa[[n]],qb[[n]]}] ],{n,npart}];
     result[[3]]=matr0 Sum[intt[ Sequence@@Flatten[{qa[[n]],qb[[n]]}] ]/
                  int[ Sequence@@Flatten[{qa[[n]],qb[[n]]}] ],{n,npart}];
     result[[4]]=matr0 Sum[neq=poseq[n];
      intee[ Sequence@@Flatten[{qa[[n1]],qa[[n2]],qb[[n1]],qb[[n2]]}] ]/
       (int[ Sequence@@Flatten[{qa[[n1]],qb[[n1]]}] ] int[ Sequence@@Flatten[{qa[[n2]],qb[[n2]]}] ]),
         {n1,npart-1},{n2,n1+1,npart}];
     Return[result//Cancel])]
);

Clear[lq,mq];
(* ! *)
Do[lq[n]=mq[n]=0,{n,npart}];
lq[5]=1;mq[5]=0;
lq[6]=1;mq[6]=0;
lq[7]=1;mq[7]=1;
lq[8]=1;mq[8]=1;

Clear[f,fn,x,p,c,nall,n1,n2,n3,n4,n5,n6,n7,n8,n9];
if="psipsi"<>ToString[npart]<>".dat";
Get[if];


mpsipsi=Length[psipsi];
qna={lq[#],mq[#],#}&/@Range[npart];

Print[Definition[mpsipsi]];Save["watchb",{npart,mpsipsi}];
{mean,meanne,meant,meanee}=
Sum[Print["m = ",m];m>>>watchb;
   psipsim=psipsi[[m]];
   {coefb,indb}=psipsim/.c_. p[nall__]->{c,{nall}};
   qnb={lq[#],mq[#],#}&/@indb;
   coefb matr[qna,qnb], {m,mpsipsi}];

wen=(meant-meanne+lambda meanee)/mean;
Print["wen"];"wen">>>watchb;
os=OpenWrite["wen"<>ToString[npart]<>".dat"(*,PageWidth->Infinity*)];
Write[os,"Energy functional for an atom with npart electrons"];
Save[os,npart,wen];
Close[os];

la=1/npart;
(* ! *)
e=wen/.{lambda->la,
($ModuleNumber=1;
   Sequence@@Table[Module[{a},A[n]->a],{n,npart}])};
Clear[a0,a1];
a0=Take[{1.,0.9,.2,.15,.15,.12,.1},npart];
a1=a0+Table[0.03,{npart}];
s=FindMinimum[ea=e;
     Print[{Release[$ModuleNumber=1;
           Table[Module[{a},a],{npart}]],
        ea}];
     {Release[$ModuleNumber=1;
           Table[Module[{a},a],{npart}]],ea}>>>watchb;
     ea,
Release[$ModuleNumber=1;
        Sequence@@Table[Module[{a},{a,{a0[[n]],a1[[n]]}}],{n,npart}]],
	WorkingPrecision->16,MaxIterations->99];
Print[npart^2 s];npart^2 s>>>watchb;

(*
a$1=1.0872;
a$2=0.7346;
a$3=
a$4=0.1326;
wen
*)