Verbosity(0)

orbs = {"1s","2s","2p"}
Ind, NF = CreateAtomicIndicesDict(orbs)

psiList = {}
for key1,orb1 in pairs(orbs) do
    for key2, orb2 in pairs(orbs) do
        if key1 == key2 then
            det1 = DeterminantString(NF, Ind[orb1])
            psiList[#psiList + 1] = NewWavefunction(NF, 0, {{det1,1}})
        elseif key1 < key2 then
            det1 = DeterminantString(NF, Ind[orb1])
            det2 = DeterminantString(NF, Ind[orb2])
            psiList[#psiList + 1] = NewWavefunction(NF, 0, {{det1,math.sqrt(1/2)}, {det2,math.sqrt(1/2)}})
        end
    end
end

Hamil = NewOperator("Number",NF,Ind["1s"],Ind["1s"],{0.1,0.2})
Hamil = Hamil + NewOperator("Number",NF,Ind["2s"],Ind["2s"],{0.3,0.4})
Hamil = Hamil + NewOperator("Number",NF,Ind["2p"],Ind["2p"],{0.5,0.6,0.7,0.8,0.9,1.0})
Hamil = Hamil + NewOperator("U",NF,Ind["1s_Up"],Ind["1s_Dn"],Ind["2p_Up"],Ind["2p_Dn"],{0.6},{0.7})
Hamil = Hamil + NewOperator("U",NF,Ind["2s_Up"],Ind["2s_Dn"],Ind["2p_Up"],Ind["2p_Dn"],{0.8},{0.9})
Hamil = Hamil + NewOperator("U",NF,Ind["2p_Up"],Ind["2p_Dn"],{0.5,sqrt(2)})
Hamil = Hamil + 1e-5*NewOperator("Lz",NF,Ind["2p_Up"],Ind["2p_Dn"])
Hamil = Hamil - 4
Hamil = Hamil/2
Hamil = Hamil + ConjugateTranspose(Hamil)
Hamil.Name = "H"

opList = {}
opList[#opList + 1] = NewOperator("Number",NF,Ind["1s"],Ind["1s"],{1,1})
opList[#opList].Name = "N_1s"
opList[#opList + 1] = NewOperator("Number",NF,Ind["2s"],Ind["2s"],{1,1})
opList[#opList].Name = "N_2s"
opList[#opList + 1] = NewOperator("Number",NF,Ind["2p"],Ind["2p"],{1,1,1,1,1,1})
opList[#opList].Name = "N_2p"
opList[#opList + 1] = NewOperator("U",NF,Ind["2p_Up"],Ind["2p_Dn"],{1,0})
opList[#opList].Name = "F0_2p"
opList[#opList + 1] = NewOperator("U",NF,Ind["2p_Up"],Ind["2p_Dn"],{0,1})
opList[#opList].Name = "F2_2p"
opList[#opList + 1] = NewOperator("U",NF,Ind["2p_Up"],Ind["2p_Dn"],{0,1})
opList[#opList].Name = "Operator with a long name"
opList[#opList + 1] = NewOperator("U",NF,Ind["2p_Up"],Ind["2p_Dn"],{0,1})
opList[#opList].Name = "An Operator with small lexicographical value"

StartRestrictions = {NF,0}
StartRestrictions[3] = {DeterminantString(NF, Ind["1s"]),2,2}
StartRestrictions[4] = {DeterminantString(NF, Ind["2s"]),2,2}

Restrictions = {NF,0}
Restrictions[3] = {DeterminantString(NF, Ind["2p"]),0,1}

NPsi = 15
psiList = Eigensystem(Hamil,StartRestrictions,NPsi)

extraList1 = {"Extra1",{}}
extraList2 = {"Extra2",{}}
for i = 1, NPsi do
    extraList1[2][i] = 0.5*i
    extraList2[2][i] = 5/i
end
extraLists = {extraList1, extraList2}

Options = {{"colWidth",10},{"printHeaderEvery",10},{"printSigma",true},{"sortByEnergy",true},{"zero",1E-7}}

print("\n\n\nPrintExpectationValues(psiList,opList):")
PrintExpectationValues(psiList,opList)
print("\n\n\nPrintExpectationValues(psiList,opList,Hamil):")
PrintExpectationValues(psiList,opList,Hamil)
print("\n\n\nPrintExpectationValues(psiList,opList,extraLists):")
PrintExpectationValues(psiList,opList,extraLists)
print("\n\n\nPrintExpectationValues(psiList,opList,Options):")
PrintExpectationValues(psiList,opList,Options)
print("\n\n\nPrintExpectationValues(psiList,opList,Hamil,extraLists):")
PrintExpectationValues(psiList,opList,Hamil,extraLists)
print("\n\n\nPrintExpectationValues(psiList,opList,Hamil,Options):")
PrintExpectationValues(psiList,opList,Hamil,Options)
print("\n\n\nPrintExpectationValues(psiList,opList,Hamil,extraLists,Options):")
PrintExpectationValues(psiList,opList,Hamil,extraLists,Options)
print("\n\n\nPrintExpectationValues(psiList,opList,-Hamil,extraLists,Options):")
PrintExpectationValues(psiList,opList,-Hamil,extraLists,Options)