-- 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(" "); 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