(* Calculations for arbitrary atom,
   initial computation of configuration *)

npart=6; (* Number of electrons *)
Off[$ModuleNumber::modnc,General::spell1];
SetDirectory["~/math/boron"];
Clear[m,n,la,lb,f,fn,x,p,alph,beta,
c,nall,pos,n1,n2,n3,n4,n5,n6,n7,n8,n9];
fn[n_]=Array[f[n][x[#]]&,npart];
alp=Array[alph,npart];
bet=Array[beta,npart];
(* ! *)
If[npart==2,
psi=Det[{alp fn[1],bet fn[2]}]-
    Det[{bet fn[1],alp fn[2]}] ];
If[npart==3,
psi=Det[{alp fn[1],bet fn[2],alp fn[3]}]-
    Det[{bet fn[1],alp fn[2],alp fn[3]}] ];
If[npart==4,
psi=Det[{alp fn[1],bet fn[2],alp fn[3],bet fn[4]}]-
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4]}]-
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4]}]+
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4]}] ];
If[npart==5,
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]}] ];
If[npart==6,
psi=Det[{alp fn[1],bet fn[2],alp fn[3],bet fn[4],alp fn[5],bet fn[6]}]-
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4],alp fn[5],bet fn[6]}]-
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4],alp fn[5],bet fn[6]}]+
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4],alp fn[5],bet fn[6]}]-
    Det[{alp fn[1],bet fn[2],alp fn[3],bet fn[4],bet fn[5],alp fn[6]}]+
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4],bet fn[5],alp fn[6]}]+
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4],bet fn[5],alp fn[6]}]-
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4],bet fn[5],alp fn[6]}]
];
If[npart==7,
psi=Det[{alp fn[1],bet fn[2],alp fn[3],bet fn[4],alp fn[5],bet fn[6],alp fn[7]}]-
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4],alp fn[5],bet fn[6],alp fn[7]}]-
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4],alp fn[5],bet fn[6],alp fn[7]}]+
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4],alp fn[5],bet fn[6],alp fn[7]}]-
    Det[{alp fn[1],bet fn[2],alp fn[3],bet fn[4],bet fn[5],alp fn[6],alp fn[7]}]+
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4],bet fn[5],alp fn[6],alp fn[7]}]+
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4],bet fn[5],alp fn[6],alp fn[7]}]-
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4],bet fn[5],alp fn[6],alp fn[7]}]
];
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=($ModuleNumber=1;
    Product[Module[{m},f[m_][x[n]]],{n,npart}])->
       ($ModuleNumber=1;
       p@@Table[Module[{m},m],{npart}]);
psi=psi1/.repl;

Print[Definition[mpsi]];Save["watcha",{npart,mpsi}];

psipsi=Sum[
   psim=psi[[m]];
   lpsi=Length[psim];
   Print["m, Length[psim] = ",{m,lpsi}];{m,lpsi}>>>watcha;
   Sum[
       {coefa,inda}=psim[[la]]/.c_. p[nall__]->{c,{nall}};
       Do[pos[n]=Position[inda,n][[1,1]],{n,npart}];
       Sum[
          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}];
(* Normalizing *)
p0=p@@Range[npart];
np0=Position[List@@psipsi,c_. p0][[1,1]];
cp0=psipsi[[np0]]/p0;
psipsi=Expand[psipsi/cp0];
Print[psipsi//Definition];
os=OpenWrite["psipsi"<>ToString[npart]<>".dat"(*,PageWidth->Infinity*)];
Write[os,"Permutations entering matrix elements for an atom with npart electrons"];
Save[os,npart,psipsi];
Close[os];