(* Calculations for arbitrary atom,
   including p-orbitals *)
(* Step b - calculation of matrix elements *)
(* Testing *)
Clear[A];
A[1]=1.0872;
A[2]=0.7346;
A[3]=A[4]=0.1326;
A[5]=0.06;

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

os=OpenWrite["test.dat",PageWidth->Infinity];
Write[os,"Matrix <1,2,...,N| {I, Vne, T, Vee} |P(1,2,...,N)>"];
Write[os,"          for an atom with N=npart electrons"];
Save[os,npart];
Print[Definition[npart]];Save["watchb",npart];
qna={lq[#],mq[#],#}&/@Range[npart];
perm=Permutations[Range[npart]];
lperm=Length[perm];
Do[Print["l = ",l];l>>>watchb;
   pl=perm[[l]];
   qnb={lq[#],mq[#],#}&/@pl;
   Clear[mat];
   Release[mat@@pl]=matr[qna,qnb];
   Save[os,mat], {l,lperm}];
Close[os];


(* Calculations for arbitrary atom,
   including p-orbitals *)

Clear[cf,mat,l];
if1="cf"<>ToString[npart]<>".dat";
if2="test.dat";
Get[if1];Get[if2];

perm=Permutations[Range[npart]];
lperm=Length[perm];
{mean,meanne,meant,meanee}=
Sum[Print["l = ",l];
   pl=perm[[l]];
   cf@@pl mat@@pl, {l,lperm}];

lambda=1/npart;
wen=(meant-meanne+lambda meanee)/mean;
Print[Definition[wen]];