print("------------------------------ Start bound states ------------------------------")
TimeStart("U")
U2 = NewOperator("AU",NF,indices,G,F,core_orbitals.kappas,valence_orbitals.kappas,{{"conserve",true},{"1P-Operator",true}}) -- Coulomb interaction between core and valence orbitals
U3 = NewOperator("AU",NF,valence_orbitals,G_vale,F_vale) -- Coulomb interaction between valence orbitals
TimeEnd("U")

H0 = NewOperator("Dirac",NF,Z,valence_orbitals,G_vale,F_vale)
U = U2 + U3
Jops = CreateAngularMomentumOperators(orbital_labels, indices,NF)
Nops = CreateNumberOperators(orbital_labels, indices,NF)
H = Chop( (H0 + U) * EnergyUnits.Ha.value + 6 * EnergyUnits.Tesla.value * (Jops.Lz + 2 * Jops.Sz) )

YtJ = YtojjzMatrix({1,2})
--Id1  = Matrix.Identity(14) --TODO: change for different valence
Id1  = Matrix.Identity(12)
ULtoJ   = Matrix.Flatten({{Id1,0},{0,YtJ}}) 

-- Crystal field Hamiltonian: -- TODO: change for different valence
IndexUp_3d = {19,21,23,25,27}
IndexDn_3d = {18,20,22,24,26}
IndexUp_2p = {13,15,17}
IndexDn_2p = {12,14,16}

Akm = {{1,-1,sqrt(1/2)},{1, 1,-sqrt(1/2)}}
TXASx = Rotate(NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm),ULtoJ)
Akm = {{1,-1,sqrt(1/2)*I},{1, 1,sqrt(1/2)*I}}
TXASy = Rotate(NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm),ULtoJ)
Akm = {{1,0,1}}
TXASz = Rotate(NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm),ULtoJ)

Akm = PotentialExpandedOnClm("Oh", 2, {0.6,-0.4})
OpptenDq = Rotate(NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, Akm),ULtoJ)
tenDq   = 2.6

H = H + tenDq * OpptenDq

H = (H + ConjugateTranspose(H))/2



psi0 = NewWavefunction(NF,NB, {{StartDeterminantString(NF,Z-Valence),1}})
shift = Braket(psi0,H,psi0) + 1
H = H - shift


TimeStart("Eigensystem")

ClosedOrbitalLabels = {"1s1/2","2s1/2","3s1/2","3p1/2","3p3/2"} 
OpenOrbitalLabels   = {"2p1/2","2p3/2"}
ExcitedOrbitalLabels= {"3d3/2","3d5/2"}
EmptyOrbitalLabels  = {}

Indices = GenerateIndices(ClosedOrbitalLabels,OpenOrbitalLabels,ExcitedOrbitalLabels,EmptyOrbitalLabels,Z,Valence)

Restrictions = {Indices.NF, NB, {DeterminantString(Indices.NF,Flatten(Indices.All)),Indices.NElectrons,Indices.NElectrons},
    {DeterminantString(Indices.NF,Flatten(Indices.Closed)),Indices.NClosedElectrons,Indices.NClosedElectrons},
    {DeterminantString(Indices.NF,Flatten(Indices.Open)),Indices.NOpenElectrons-1,Indices.NOpenElectrons},
    {DeterminantString(Indices.NF,Flatten(Indices.Excited)),0,1}}

print("Restrictions:\n")
print(Restrictions,"\n")

phi0 = Chop(Eigensystem(H,Restrictions,61,{{"Restrictions",Restrictions}})) -- ,{{"Restrictions",StartRestrictions},{"DenseBorder",2601}}

--ExpectationValues({phi0[1],phi0[2],phi0[3],phi0[4]},H,shift,Jops,Nops,Indices.Labels)--{phi0[1],phi0[2],phi0[3],phi0[4]}

-- now we select all states that can be excited from the ground-state by X-ray absorption
olap = Chop(Braket(TXASx*phi0[1],phi0),1E-6)
phiX0 = {}
idx = 1
for i = 1, #olap do
    if olap[i] ~= 0 then
        phiX0[idx] = phi0[i]
        idx = idx + 1
    end
end 

-- from these we select the states that are in the energy region of interest:
E0 = BraketDiagonal(phiX0,H,phiX0)

phiX = {}
idx = 1
for i = 1, #phiX0 do
   if E0[i] >= 450 and E0[i] <= 500 then
       phiX[idx] = phiX0[i]
       idx = idx + 1
   end
end

olapX = Chop(Braket(TXASx*phi0[1],phiX),1E-6)

--PrintExpectationValues(Join({phi0[1]},phiX),{Jops.Jsqr,Jops.Lsqr,Jops.Ssqr,Jops.Jz},H,{{"residue",Join({0},olapX)}})
poles = BraketDiagonal(phiX,H,phiX)
residues = Chop(Transpose({Braket(TXASx*phi0[1],phiX)}),1E-6)

ClosedOrbitalLabels = {"1s1/2","2s1/2"}
OpenOrbitalLabels   = {"3s1/2","3p1/2","3p3/2","2p1/2","2p3/2"}
ExcitedOrbitalLabels= {"3d3/2","3d5/2"}
EmptyOrbitalLabels  = {} 

Indices = GenerateIndices(ClosedOrbitalLabels,OpenOrbitalLabels,ExcitedOrbitalLabels,EmptyOrbitalLabels,Z,Valence)
Indices.Closed.kappas = Take(Indices.All.kappas,1,#Indices.Closed)
Indices.Open.kappas   = Take(Indices.All.kappas,#Indices.Closed+1,#Indices.All)


StartRestrictions = {Indices.NF, NB, {DeterminantString(Indices.NF,Flatten(Indices.All)),Indices.NElectrons-1,Indices.NElectrons-1},
                     {DeterminantString(Indices.NF,Flatten(Indices.Closed)),Indices.NClosedElectrons,Indices.NClosedElectrons},
                     {DeterminantString(Indices.NF,Flatten(Indices.Open)),Indices.NOpenElectrons-2,Indices.NOpenElectrons-1}}

print("StartRestrictions:\n")
print(StartRestrictions,"\n")

print("---------------------------------------------------------------------------------------------------------------------\n")

psi1 = Chop(Eigensystem(H,StartRestrictions,924,{{"Restrictions",StartRestrictions},{"DenseBorder",2601}})) 

E1 = BraketDiagonal(psi1,H,psi1)

psi = {}
EAuger = {}
i = 1
while i <= #E1 and E1[i] <= 500 do
    psi[i] = psi1[i]
    EAuger[i] = E1[i]
    i = i +1
end

ExpectationValues(Join({phi0[1]},psi),H,shift,Jops,Nops,Indices.Labels)

TimeEnd("Eigensystem")
print("------------------------------   End bound states ------------------------------")
