(* Calculations for arbitrary atom,
   including p-orbitals *)
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];
(* ! *) npart=5;
(* 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];
(* ! *)
lq[1]=0;mq[1]=0;
lq[2]=0;mq[2]=0;
lq[3]=0;mq[3]=0;
lq[4]=0;mq[4]=0;
lq[5]=1;mq[5]=0;

Clear[f,fn,x,p,c,nall,n1,n2,n3,n4,n5,n6,n7,n8,n9];
fn[n_]=Array[f[n][x[#]]&,npart];
alp=Array[alph,npart];
bet=Array[beta,npart];
(* ! *)
psi=Det[{alp fn[1],bet fn[2],alp fn[3],bet fn[4],alp fn[5]}]-
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4],alp fn[5]}]-
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4],alp fn[5]}]+
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4],alp fn[5]}] ;
spins=Flatten[Table[{alph[n],beta[n]},{n,npart}]];
psi=Flatten[CoefficientList[psi,spins]];
psi1=Complement[psi,{0}];
(* Test that elements of psi1 enter psi equal number of times *)
If[Length[Union[Map[Count[psi,#]&,psi1]]]==1,psi=psi1,
Print["Warning: Something strange with wavefunction!"]];
If[Length[Complement[Map[Count[psi1,-#]&,psi1],{0}]]>0,
Print["Further removing of elements with an opposite sign is possible (but not done)."]];
mpsi=Length[psi1];
(* ! *)
repl=f[n1_][x[1]] f[n2_][x[2]] f[n3_][x[3]] f[n4_][x[4]] f[n5_][x[5]] -> p[n1,n2,n3,n4,n5];
psi=psi1/.repl;

psipsi=Sum[
   psim=psi[[m]];
   lpsi=Length[psim];
   Sum[{coefa,inda}=psim[[la]]/.c_. p[nall__]->{c,{nall}};
       Clear[pos];
       Do[pos[n]=Position[inda,n][[1,1]],{n,npart}];
       Sum[Print["m, la, lb = ",{m,la,lb}];
          indb=Operate[List&,psim[[lb]]];
          {coefb,indb}=psim[[lb]]/.c_. p[nall__]->{c,{nall}};
          ind=Table[indb[[pos[n]]],{n,npart}];
          coefa coefb p@@ind, {lb,lpsi}]
       ,{la,lpsi}]
   ,{m,mpsi}];
Print[psipsi//Definition];
mpsipsi=Length[psipsi];
qna={lq[#],mq[#],#}&/@Range[npart];

{mean,meanne,meant,meanee}=
Sum[Print["m = ",m];
   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"];
SetDirectory["~/temp"];
Save["temp.m",wen];

la=1/npart;
(* ! *)
Clear[a1,a2,a3,a4,a5,a6,a7,a8,a9];
e=wen/.{lambda->la,A[1]->a1,A[2]->a2,A[3]->a3,A[4]->a4,A[5]->a5};
{a10,a20,a30,a40,a50}={1.,1.,.1,.1,.1};
{a11,a21,a31,a41,a51}={a10,a20,a30,a40,a50}+Table[0.03,{npart}];
s=FindMinimum[ea=e;Print[{a1,a2,ea}];
		ea,{a1,{a10,a11}},{a2,{a20,a21}},{a3,{a30,a31}},{a4,{a40,a41}},
                {a5,{a50,a51}},
	WorkingPrecision->16,MaxIterations->99];
npart^2 s