(* ::Package:: *)

BeginPackage["OrbitalEntanglement`",{"SqueezeRoutines`"}]

(*

Eddy Ardonne and Nicolas Regnault

Squeeze Routines for polarized and spinful quantum Hall states,
bosonic as well as fermionic. Routines for the spinful case
are based on

The structure of spinful quantum Hall states: a squeezing perspective,
E. Ardonne, N. Regnault,
arXiv:1107.2232

See this paper for references for the polarized case, and for reference
about particle and orbital entanglement spectra.

Routines making use of the additional symmetries in case of L=0 and/or
S=0 states are not implemented in this version. A routine to calculate
only the number of solutions to the consistency condition via an
eigenvalue problem is not included either. Please contact E. Ardonne
about these (and possible other) issues.

*)


oes::usage=
 "oes[la] will calculate the orbital entanglement spectrum, with
  system A consisting of the first la orbitals for spinless bosonic
  and fermionic states."

spinoes::usage=
 "spinoes[la] will calculate the particle entanglement spectrum with
  system A consisting of the first la orbitals, for spinful bosonic
  and fermionic states."


Begin["`Private`"]

(** Some functions, needed to generate the root configurations, which
are subsequently used to generate the Hilbertspaces system A and B **)

lz[vec_]:=With[{norb=Length[vec]},Sum[(-(norb-1)/2-1+i)vec[[i]],{i,1,norb}]];

id[num_]:=integerdigits[num,Global`base,Global`norbitals]
fd[rp_]:=fromdigits[rp,Global`base]
sid[num_]:={id[num[[1]]],id[num[[2]]]}
sfd[rp_]:={fd[rp[[1]]],fd[rp[[2]]]}
disknorm[flux_,l_]:=1/Sqrt[2^l l!];
spherenorm[flux_,l_]:=Sqrt[Binomial[flux,l]];
cftnorm[flux_,l_]:=1;


(* default value for the normalization *)

Global`normalization="sphere";


(** Routines to calculate the orbital entanglement spectrum **)

(* Spinless states, both bosonic as well as fermionic *)

oes[na_]:=With[{},

(* Set the normalization *)
Piecewise[{
{onepartnorm[flux_,l_]:=spherenorm[flux,l];,Global`normalization=="sphere"},
{onepartnorm[flux_,l_]:=cftnorm[flux,l];,Global`normalization=="cft"},
{onepartnorm[flux_,l_]:=disknorm[flux,l];,Global`normalization=="disk"}}
,onepartnorm[flux_,l_]:=spherenorm[flux,l];];

base=Global`base;
nparticles=Global`nparticles;
norbitals=Global`norbitals;
hsdim=Global`hsdim;

(* Normalizing the state *)

hsnorm=Table[With[{vec=id[Global`state[i]]},1/Product[onepartnorm[norbitals-1,j-1]^(vec[[j]])//N,
{j,1,norbitals}]],{i,1,hsdim}];

deg=Table[With[{vec=id[Global`state[i]]},Sqrt[nparticles!/Product[id[Global`state[i]][[j]]!//N,{j,1,norbitals}]]],{i,1,hsdim}];

norm=Sqrt[Sum[(hsnorm[[i]]*deg[[i]]* Global`b[i])^2,{i,1,hsdim}]];

(* Initializing angular momentum varialbles *)
lt=(norbitals-1)/2;
lstate=Sum[(lt-j)id[Global`state[1]][[j+1]],{j,0,norbitals-1}];

(* The number of particles in part B *)
nb=norbitals-na;


(* Calculating the reduced density matrices, and their eigenvalues *)

Clear[astate,bstate,astatepresent,bstatepresent,astatepos,bstatepos];
Clear[nlpair,nlpairpresent,nlpairapos,nlpairbpos,npairs,Global`ndata,cmat,svd];
adim=0;
bdim=0;
numnlpairs=0;


Do[
acurstate=base^nb fromdigits[id[Global`state[i]][[1;;na]],base];
bcurstate=fromdigits[id[Global`state[i]][[na+1;;norbitals]],base];
If[astatepresent[acurstate],Null,Null,
astatepresent[acurstate]=True;
adim++;
astate[adim]=acurstate;
astatepos[acurstate]=adim;];
If[bstatepresent[bcurstate],Null,Null,
bstatepresent[bcurstate]=True;
bdim++;
bstate[bdim]=bcurstate;
bstatepos[bcurstate]=bdim;];
currentnlpair={Plus@@id[Global`state[i]][[1;;na]],Sum[(lt-j)id[Global`state[i]][[j+1]],{j,0,na-1}]};
If[nlpairpresent[currentnlpair],
nlpairapos[currentnlpair]=Append[nlpairapos[currentnlpair],astatepos[acurstate]];
nlpairbpos[currentnlpair]=Append[nlpairbpos[currentnlpair],bstatepos[bcurstate]];
,Null,
nlpairpresent[currentnlpair]=True;numnlpairs++;nlpair[numnlpairs]=currentnlpair;
nlpairapos[currentnlpair]={astatepos[acurstate]};
nlpairbpos[currentnlpair]={bstatepos[bcurstate]};]
,{i,1,hsdim}];

Do[nlpairapos[nlpair[i]]=Union[nlpairapos[nlpair[i]]];nlpairbpos[nlpair[i]]=Union[nlpairbpos[nlpair[i]]],{i,1,numnlpairs}];

Do[cmat[nlpair[i]]=Table[With[{hspos=Global`stateposition[astate[i1]+bstate[j1]]},If[NumberQ[hspos],Global`b[hspos]hsnorm[[hspos]]deg[[hspos]]/norm,0]],{i1,nlpairapos[nlpair[i]]},{j1,nlpairbpos[nlpair[i]]}],{i,1,numnlpairs}];

Do[svd[nlpair[i]]=SingularValueList[cmat[nlpair[i]]//N],{i,1,numnlpairs}];

checknorm=Plus@@Flatten[Table[#^2&/@svd[nlpair[i]],{i,1,numnlpairs}]];
Print["Checking the norm gives: ",checknorm];

ententropy=-Plus@@Flatten[Table[2#^2Log[#]&/@DeleteCases[Chop[svd[nlpair[i]]],0],{i,1,numnlpairs}]];
Print["The entanglement entropy is: ",ententropy];


(* Collecting the data *)

nvalues=Union[Table[nlpair[i],{i,1,numnlpairs}][[All,1]]];

Global`nvalues=nvalues;
Do[npairs[i]=Cases[Table[nlpair[i],{i,1,numnlpairs}],x_/;x[[1]]==i],{i,nvalues}];

Clear[Global`ndata];
Do[Global`ndata[i]=DeleteCases[Flatten[Table[{j[[2]],-2Log[v]},{j,npairs[i]},{v,svd[{i,j[[2]]}]}],1],x_/;x[[2]]===Indeterminate],{i,nvalues}];

lvalues=Union[Table[nlpair[i],{i,1,numnlpairs}][[All,2]]];

Global`nlzstates=Table[If[n>=Min[nvalues],Count[Global`ndata[n],x_/;x[[1]]==l],0],{n,0,Max[nvalues]},{l,0,Max[lvalues],1/Length[Union[Mod[#,1]&/@lvalues]]}];

Global`lzstates=Plus@@Global`nlzstates;


];


(* Spinful states, both bosonic as well as fermionic *)

spinoes[na_]:=With[{},
(* Set the normalization *)
Piecewise[{
{onepartnorm[flux_,l_]:=spherenorm[flux,l];,Global`normalization=="sphere"},
{onepartnorm[flux_,l_]:=cftnorm[flux,l];,Global`normalization=="cft"},
{onepartnorm[flux_,l_]:=disknorm[flux,l];,Global`normalization=="disk"}}
,onepartnorm[flux_,l_]:=spherenorm[flux,l];];

base=Global`base;
nparticles=Global`nparticles;
norbitals=Global`norbitals;
hsdim=Global`hsdim;

(* Normalizing the state *)

hsnorm=Table[With[{vec=sid[Global`state[i]]},1/Product[onepartnorm[norbitals-1,j-1]^(vec[[1,j]]+vec[[2,j]])//N,
{j,1,norbitals}]],{i,1,hsdim}];

deg=Table[With[{vec=sid[Global`state[i]]},Sqrt[nparticles!/Product[vec[[1,j]]!*vec[[2,j]]!,{j,1,norbitals}]]//N],{i,1,hsdim}];

norm=Sqrt[Sum[(hsnorm[[i]]*deg[[i]]* Global`b[i])^2,{i,1,hsdim}]];

(* Initializing angular momentum varialbles *)
lt=(norbitals-1)/2;
lstate=Sum[(lt-j)*(sid[Global`state[1]][[1,j+1]]+sid[Global`state[1]][[2,j+1]]),{j,0,norbitals-1}];


(* The number of particles in part B *)
nb=norbitals-na;


(* Calculating the reduced density matrices, and their singular values *)

Clear[astate,bstate,astatepresent,bstatepresent,astatepos,bstatepos];
Clear[nlspair,nlspairpresent,nlspairapos,nlspairbpos];
adim=0;
bdim=0;
numnlspairs=0;

Do[
acurstate={base^nb fromdigits[sid[Global`state[i]][[1,1;;na]],base],base^nb fromdigits[sid[Global`state[i]][[2,1;;na]],base]};
bcurstate={fromdigits[sid[Global`state[i]][[1,na+1;;norbitals]],base],fromdigits[sid[Global`state[i]][[2,na+1;;norbitals]],base]};
If[astatepresent[acurstate],Null,Null,
astatepresent[acurstate]=True;
adim++;
astate[adim]=acurstate;
astatepos[acurstate]=adim;];
If[bstatepresent[bcurstate],Null,Null,
bstatepresent[bcurstate]=True;
bdim++;
bstate[bdim]=bcurstate;
bstatepos[bcurstate]=bdim;];
currentnlspair={Plus@@sid[Global`state[i]][[1,1;;na]]+Plus@@sid[Global`state[i]][[2,1;;na]],Sum[(lt-j)sid[Global`state[i]][[1,j+1]]+(lt-j)sid[Global`state[i]][[2,j+1]],{j,0,na-1}],1/2(Plus@@sid[Global`state[i]][[1,1;;na]]-Plus@@sid[Global`state[i]][[2,1;;na]])};
If[nlspairpresent[currentnlspair],
nlspairapos[currentnlspair]=Append[nlspairapos[currentnlspair],astatepos[acurstate]];
nlspairbpos[currentnlspair]=Append[nlspairbpos[currentnlspair],bstatepos[bcurstate]];
,Null,
nlspairpresent[currentnlspair]=True;numnlspairs++;nlspair[numnlspairs]=currentnlspair;
nlspairapos[currentnlspair]={astatepos[acurstate]};
nlspairbpos[currentnlspair]={bstatepos[bcurstate]};]
,{i,1,hsdim}];

Do[nlspairapos[nlspair[i]]=Union[nlspairapos[nlspair[i]]];nlspairbpos[nlspair[i]]=Union[nlspairbpos[nlspair[i]]],{i,1,numnlspairs}];

Do[cmat[nlspair[i]]=Table[With[{hspos=Global`stateposition[astate[i1]+bstate[j1]]},If[NumberQ[hspos],Global`b[hspos]hsnorm[[hspos]]deg[[hspos]]/norm,0]],{i1,nlspairapos[nlspair[i]]},{j1,nlspairbpos[nlspair[i]]}],{i,1,numnlspairs}];

Do[svd[nlspair[i]]=SingularValueList[cmat[nlspair[i]]//N],{i,1,numnlspairs}];

checknorm=Plus@@Flatten[Table[#^2&/@svd[nlspair[i]],{i,1,numnlspairs}]];
Print["Checking the norm gives: ",checknorm];

ententropy=-Plus@@Flatten[Table[2#^2Log[#]&/@DeleteCases[Chop[svd[nlspair[i]]],0],{i,1,numnlspairs}]];
Print["The entanglement entropy is: ",ententropy];

(* Collecting the data *)

(*
nvalues=Union[Table[nlspair[i],{i,1,numnlspairs}][[All,1]]];
Global`nvalues=nvalues;
svalues=Union[Table[nlspair[i],{i,1,numnlspairs}][[All,3]]];
Global`svalues=svalues;
Do[npairs[i]=Cases[Table[nlspair[i],{i,1,numnlspairs}],x_/;x[[1]]==i],{i,nvalues}];

nsvalues=Union[Table[nlspair[i],{i,1,numnlspairs}][[All,{1,3}]]];
Do[nspairs[i]=Cases[Table[nlspair[i],{i,1,numnlspairs}],x_/;x[[1]]==i[[1]]&&x[[3]]==i[[2]]],{i,nsvalues}];
*)

nlspairlist=Table[nlspair[i],{i,1,numnlspairs}];
Global`nlspairlist=nlspairlist;
nspairlist=Union[nlspairlist[[All,{1,3}]]];

Clear[Global`nsdata]
Do[Global`nsdata[i]=Flatten[Table[DeleteCases[Table[{j[[2]],-2Log[v]},{v,svd[j]}],x_/;x[[2]]===Indeterminate],{j,Cases[nlspairlist,x_/;x[[1]]==i[[1]]&&x[[3]]==i[[2]]]}],1],{i,nspairlist}];
(*
Clear[Global`nsdata];
Do[Global`nsdata[i]=DeleteCases[Table[{i[[2]],-2Log[v]},{v,svd[i]}],x_/;x[[2]]===Indeterminate],{i,nlspairlist}];
*)

nspairs=Cases[nspairlist,x_/;Global`nsdata[x]!={}];
Global`nspairs=nspairs;
Global`nvalues=Union[nspairs[[All,1]]];
Global`svalues=Union[nspairs[[All,2]]];
Global`lvalues=Union[Flatten[Table[Cases[nlspairlist,x_/;x[[1]]==i[[1]]&&x[[3]]==i[[2]]],{i,nspairs}],1][[All,2]]];

Global`nlzszstates=Table[If[Head[Global`nsdata[{n,s}]]===List,Count[Global`nsdata[{n,s}],x_/;x[[1]]==l],0],{n,0,Max[Global`nvalues]},{s,0,Max[Global`svalues],1/2},{l,0,Max[Global`lvalues],1/Length[Union[Mod[#,1]&/@Global`lvalues]]}];


];


End[ ]

EndPackage[ ]

