-- Using operators and wavefunctions as explained in 
-- the Operators and Wavefunctions example
-- and being able to multiply them to get
-- expectation values we can continue and look
-- at eigenstates of operators

-- define the basis
  -- For a p-shell we would like the have 6 
  -- spinorbitals, with the quantum numbers 
  -- spin up ml=-1,ml=0,ml=1 and
  -- spin down with ,ml=-1, ml=0, ml=1
NF=6
NB=0
IndexDn={0,2,4}
IndexUp={1,3,5}

-- Define spin and angular momentum operators.
OppSx   =NewOperator("Sx"   ,NF,IndexUp,IndexDn)
OppSy   =NewOperator("Sy"   ,NF,IndexUp,IndexDn)
OppSz   =NewOperator("Sz"   ,NF,IndexUp,IndexDn)
OppSsqr =NewOperator("Ssqr" ,NF,IndexUp,IndexDn)
OppSplus=NewOperator("Splus",NF,IndexUp,IndexDn)
OppSmin =NewOperator("Smin" ,NF,IndexUp,IndexDn)

OppLx   =NewOperator("Lx"   ,NF,IndexUp,IndexDn)
OppLy   =NewOperator("Ly"   ,NF,IndexUp,IndexDn)
OppLz   =NewOperator("Lz"   ,NF,IndexUp,IndexDn)
OppLsqr =NewOperator("Lsqr" ,NF,IndexUp,IndexDn)
OppLplus=NewOperator("Lplus",NF,IndexUp,IndexDn)
OppLmin =NewOperator("Lmin" ,NF,IndexUp,IndexDn)

OppJx   =NewOperator("Jx"   ,NF,IndexUp,IndexDn)
OppJy   =NewOperator("Jy"   ,NF,IndexUp,IndexDn)
OppJz   =NewOperator("Jz"   ,NF,IndexUp,IndexDn)
OppJsqr =NewOperator("Jsqr" ,NF,IndexUp,IndexDn)
OppJplus=NewOperator("Jplus",NF,IndexUp,IndexDn)
OppJmin =NewOperator("Jmin" ,NF,IndexUp,IndexDn)

Oppldots=NewOperator("ldots",NF,IndexUp,IndexDn)

-- Define the coulomb operator
-- We here define the part depending on F0 
-- separately from the part depending on F2.
-- When summing we can put in the numerical values 
-- of the slater integrals.
OppF0 = NewOperator("U",NF,IndexUp,IndexDn,{1,0})
OppF2 = NewOperator("U",NF,IndexUp,IndexDn,{0,1})
OppU = 5.0 * OppF0 + 4.0 * OppF2

-- Note that the previous definition is the same as
-- OppU = NewOperator("U", NF, IndexUp, IndexDn, 
--                         {5.0,4.0})
 
-- Define the Hamiltonian as a numerical sum of the 
-- previous defined operators.
Hamiltonian = 5.0 * OppF0 + 4.0 * OppF2 + 0.000001 * (2*OppSz + OppLz)
Hamiltonian2 = 6.0 * OppF0 + 4.0 * OppF2 + 0.000001 * (2*OppSz + OppLz)

-- For large systems we do not need to know all
-- eigenstates, but can restrict ourselves to the
-- Npsi lowest states:
Npsi=15

-- In order to make sure we have a filling of 2
-- electrons we need to define some restrictions
StartRestrictions = {NF, NB,  {"111111",2,2}}

-- We now can create the lowest Npsi eigenstates:
psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)

oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, Oppldots, OppF0, OppF2}

-- after we've created the eigen states we can look 
-- at a list of their expectation values
print(" <E>     <S^2>   <L^2>   <J^2>   <S_z>   <L_z>   <l.s>   <F[0]>  <F[2]>");
for key,value in pairs(psiList) do
  expvalue = value * oppList * value
  for k,v in pairs(expvalue) do
    io.write(string.format("%7.4f ",v))
  end;
  io.write("\n")
end