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

npart=6; (* Number of electrons *)
Off[$ModuleNumber::modnc,General::spell1];
If[$VersionNumber==2.2,
   SetDirectory["c:\sergeev\purdue\math\boron"],
   SetDirectory["~/math/boron"] ];
Clear[m,n,la,lb,f,fn,x,p,alph,beta,
c,nall,pos,coeff];
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],alp fn[6]}]-
    Det[{bet fn[1],alp fn[2],alp fn[3],bet fn[4],alp fn[5],alp fn[6]}]-
    Det[{alp fn[1],bet fn[2],bet fn[3],alp fn[4],alp fn[5],alp fn[6]}]+
    Det[{bet fn[1],alp fn[2],bet fn[3],alp fn[4],alp fn[5],alp fn[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]];
{npart,mpsi}>>watcha;
perm=Permutations[Range[npart]];
lperm=Length[perm];
Do[coeff@perm[[l]]=0,{l,lperm}];

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

os=OpenWrite["cf"<>ToString[npart]<>".dat",PageWidth->Infinity];
Save[os,npart];
Write[os,"Permutations (with their weights) entering matrix elements for an atom with npart electrons"];

(* Normalizing *)
c0=coeff@Range[npart];
Do[pl=perm[[l]];
   Clear[cf];
   Release[cf@@pl]=coeff@pl/c0;
   Save[os,cf],{l,lperm}];
Close[os];