====== Reuse every single determinant constituting the eigenstate ====== ;;# asked by [[mailto:ruiwen.xie@tu-darmstadt.de|Ruiwen Xie]] (2021/11/11 12:20) ;;# == == Dear Quanty developers, I am now trying to use quanty to build the model proposed by A.Kotani in PHYSICAL REVIEW B 78, 195115, 2008, in order to get the L23 XAS and XMCD for Ce-based compounds with Ce mix-valency. First I got the eigenstates for Ce 4f^0 and 4f^1 configurations (the script on this forum from BODRY TEGOMO CHIOGO is of great help for this part). The following is the part of script to get eigenstates Verbosity(0) NF = 44 NB = 0 IndexDn_4f = {0,2,4,6,8,10,12} IndexUp_4f = {1,3,5,7,9,11,13} IndexDn_Lf = {14,16,18,20,22,24,26} IndexUp_Lf = {15,17,19,21,23,25,27} IndexDn_2p = {28,30,32} IndexUp_2p = {29,31,33} IndexDn_5d = {34,36,38,40,42} IndexUp_5d = {35,37,39,41,43} NElectrons_4f = 0 V = 0.30001 --hybridization between 4f and ligand Delta_4f_Lf_i = -1.0 -- the 4f energy E_4f1L - E_4f0 NConfigurations = 2 Npsis =15 OppSx_4f = NewOperator("Sx", NF, IndexUp_4f,IndexDn_4f) OppSy_4f = NewOperator("Sy", NF, IndexUp_4f,IndexDn_4f) OppSz_4f = NewOperator("Sz", NF, IndexUp_4f,IndexDn_4f) OppSsqr_4f = NewOperator("Ssqr", NF, IndexUp_4f,IndexDn_4f) OppSplus_4f = NewOperator("Splus", NF, IndexUp_4f,IndexDn_4f) OppSmin_4f = NewOperator("Smin", NF, IndexUp_4f,IndexDn_4f) OppLx_4f = NewOperator("Lx", NF, IndexUp_4f,IndexDn_4f) OppLy_4f = NewOperator("Ly", NF, IndexUp_4f,IndexDn_4f) OppLz_4f = NewOperator("Lz", NF, IndexUp_4f,IndexDn_4f) OppLsqr_4f = NewOperator("Lsqr", NF, IndexUp_4f,IndexDn_4f) OppLplus_4f = NewOperator("Lplus", NF, IndexUp_4f,IndexDn_4f) OppLmin_4f = NewOperator("Lmin", NF, IndexUp_4f,IndexDn_4f) OppJx_4f = NewOperator("Jx", NF, IndexUp_4f,IndexDn_4f) OppJy_4f = NewOperator("Jy", NF, IndexUp_4f,IndexDn_4f) OppJz_4f = NewOperator("Jz", NF, IndexUp_4f,IndexDn_4f) OppJsqr_4f = NewOperator("Jsqr", NF, IndexUp_4f,IndexDn_4f) OppJplus_4f = NewOperator("Jplus", NF, IndexUp_4f,IndexDn_4f) OppJmin_4f = NewOperator("Jmin", NF, IndexUp_4f,IndexDn_4f) Oppldots_4f = NewOperator("ldots", NF, IndexUp_4f,IndexDn_4f) -- Angular momentum operator for the Ligand f-shell OppSx_Lf = NewOperator("Sx", NF, IndexUp_Lf,IndexDn_Lf) OppSy_Lf = NewOperator("Sy", NF, IndexUp_Lf,IndexDn_Lf) OppSz_Lf = NewOperator("Sz", NF, IndexUp_Lf,IndexDn_Lf) OppSsqr_Lf = NewOperator("Ssqr", NF, IndexUp_Lf,IndexDn_Lf) OppSplus_Lf = NewOperator("Splus", NF, IndexUp_Lf,IndexDn_Lf) OppSmin_Lf = NewOperator("Smin", NF, IndexUp_Lf,IndexDn_Lf) OppLx_Lf = NewOperator("Lx", NF, IndexUp_Lf,IndexDn_Lf) OppLy_Lf = NewOperator("Ly", NF, IndexUp_Lf,IndexDn_Lf) OppLz_Lf = NewOperator("Lz", NF, IndexUp_Lf,IndexDn_Lf) OppLsqr_Lf = NewOperator("Lsqr", NF, IndexUp_Lf,IndexDn_Lf) OppLplus_Lf = NewOperator("Lplus", NF, IndexUp_Lf,IndexDn_Lf) OppLmin_Lf = NewOperator("Lmin", NF, IndexUp_Lf,IndexDn_Lf) OppJx_Lf = NewOperator("Jx", NF, IndexUp_Lf,IndexDn_Lf) OppJy_Lf = NewOperator("Jy", NF, IndexUp_Lf,IndexDn_Lf) OppJz_Lf = NewOperator("Jz", NF, IndexUp_Lf,IndexDn_Lf) OppJsqr_Lf = NewOperator("Jsqr", NF, IndexUp_Lf,IndexDn_Lf) OppJplus_Lf = NewOperator("Jplus", NF, IndexUp_Lf,IndexDn_Lf) OppJmin_Lf = NewOperator("Jmin", NF, IndexUp_Lf,IndexDn_Lf) Oppldots_Lf = NewOperator("ldots", NF, IndexUp_4f,IndexDn_4f) -- SUM OF THE OPERATOR OppSx = OppSx_4f + OppSx_Lf OppSy = OppSy_4f + OppSy_Lf OppSz = OppSz_4f + OppSz_Lf OppSsqr = OppSx*OppSx + OppSy*OppSy + OppSz*OppSz OppLx = OppLx_4f + OppLx_Lf OppLy = OppLy_4f + OppLy_Lf OppLz = OppLz_4f + OppLz_Lf OppLsqr = OppLx*OppLx + OppLy*OppLy + OppLz*OppLz OppJx = OppJx_4f + OppJx_Lf OppJy = OppJy_4f + OppJy_Lf OppJz = OppJz_4f + OppJz_Lf OppJsqr = OppJx*OppJx + OppJy*OppJy + OppJz*OppJz OppF0_4f = NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {1, 0, 0, 0}) OppF2_4f = NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {0, 1, 0, 0}) OppF4_4f = NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {0, 0, 1, 0}) OppF6_4f = NewOperator("U", NF, IndexUp_4f, IndexDn_4f, {0, 0, 0, 1}) OppNUp_4f = NewOperator("Number", NF, IndexUp_4f,IndexUp_4f,{1,1,1,1,1,1,1}) OppNDn_4f = NewOperator("Number", NF, IndexDn_4f,IndexDn_4f,{1,1,1,1,1,1,1}) OppN_4f = OppNUp_4f + OppNDn_4f OppNUp_Lf = NewOperator("Number", NF, IndexUp_Lf,IndexUp_Lf,{1,1,1,1,1,1,1}) OppNDn_Lf = NewOperator("Number", NF, IndexDn_Lf,IndexDn_Lf,{1,1,1,1,1,1,1}) OppN_Lf = OppNUp_Lf + OppNDn_Lf zeta_4f_i = 0.087* 0.92 zeta_Lf_i = 0.087* 0.92 U_4f_4f_i = 0 F2_4f_4f_i = 0 * 0.55 F4_4f_4f_i = 0 * 0.55 F6_4f_4f_i = 0 * 0.55 F0_4f_4f_i = U_4f_4f_i + 4 / 195 * F2_4f_4f_i + 2 / 143 * F4_4f_4f_i + 100 / 5577 * F6_4f_4f_i Bz = 0.001 H_i = Chop(F0_4f_4f_i *OppF0_4f + F2_4f_4f_i *OppF2_4f+ F4_4f_4f_i *OppF4_4f + F6_4f_4f_i *OppF6_4f + zeta_4f_i * Oppldots_4f + zeta_Lf_i * Oppldots_Lf + Bz*(2*OppSz_4f + OppLz_4f)) e_4f_i = (28 * Delta_4f_Lf_i - 27 * U_4f_4f_i * NElectrons_4f - U_4f_4f_i * NElectrons_4f^2) / (2 * (14 + NElectrons_4f)) e_Lf_i = NElectrons_4f * (-2 * Delta_4f_Lf_i + U_4f_4f_i * NElectrons_4f + U_4f_4f_i) / (2 * (NElectrons_4f + 14)) H_i = H_i + Chop( e_4f_i * OppN_4f + e_Lf_i * OppN_Lf) -- hybridization Hamiltonian OppV = NewOperator("Number", NF, IndexUp_4f,IndexUp_Lf, {1,1,1,1,1,1,1}) + NewOperator("Number", NF, IndexDn_4f,IndexDn_Lf, {1,1,1,1,1,1,1})+ NewOperator("Number", NF, IndexUp_Lf,IndexUp_4f, {1,1,1,1,1,1,1})+ NewOperator("Number", NF, IndexDn_Lf,IndexDn_4f, {1,1,1,1,1,1,1}) H_i = H_i +Chop(V*OppV) InitialRestrictions = {NF, NB, {" 11111111111111 00000000000000 000000 0000000000",0,1}, {" 00000000000000 11111111111111 000000 0000000000",14 - (NConfigurations - 1), 14}, {" 00000000000000 00000000000000 111111 0000000000",6,6}, {" 00000000000000 00000000000000 000000 1111111111",0,0}, {" 11111111111111 11111111111111 000000 0000000000",14,14}} CalculationRestrictions = {NF, NB, {" 11111111111111 00000000000000 000000 0000000000",0,1}, {" 00000000000000 11111111111111 000000 0000000000",14 - (NConfigurations - 1), 14}, {" 00000000000000 00000000000000 111111 0000000000",6,6}, {" 00000000000000 00000000000000 000000 1111111111",0,0}, {" 11111111111111 11111111111111 000000 0000000000",14,14}} psiList = Eigensystem(H_i, InitialRestrictions, Npsis, {{'restrictions', CalculationRestrictions}}) print(Chop(psiList[1])) -- weight of each configuration Oppf1 = NewOperator ( "Identity" , NF ) Oppf1.Restrictions = { NF , NB , { " 11111111111111 00000000000000 000000 0000000000" , 1 , 1 } } Oppf0 = NewOperator ( "Identity" , NF ) Oppf0.Restrictions = { NF , NB , { " 11111111111111 00000000000000 000000 0000000000" , 0 , 0 } } oppList={H_i, OppSsqr, OppLsqr,OppLz_4f, OppJsqr, OppJz_4f, Oppldots_4f, OppF2_4f, OppF4_4f, OppF6_4f,OppN_4f, OppN_Lf, Oppf1, Oppf0, OppN_5d} print(" # "); for i = 1,#psiList do io.write(string.format("%3i ",i)) for j = 1,#oppList do expectationvalue = Chop(psiList[i]*oppList[j]*psiList[i]) io.write(string.format("%10.4f ",expectationvalue)) end io.write("\n") end Then the output I got is as follows WaveFunction: Wave Function QComplex = 0 (Real==0 or Complex==1) N = 27 (Number of basis functions used to discribe psi) NFermionic modes = 44 (Number of fermions in the one particle basis) NBosonic modes = 0 (Number of bosons in the one particle basis) # pre-factor Determinant 1 -5.799417965437E-02 00000000000100111111111111011111110000000000 2 -7.508426285108E-02 00000000010000111111111101111111110000000000 3 -8.248624204609E-02 00000001000000111111110111111111110000000000 4 -8.272282066785E-02 00000100000000111111011111111111110000000000 5 -7.573216559561E-02 00010000000000111101111111111111110000000000 6 -5.883063271692E-02 01000000000000110111111111111111110000000000 7 5.799417965437E-02 00000000000010111111111110111111110000000000 8 7.508426285107E-02 00000000001000111111111011111111110000000000 9 8.248624204609E-02 00000000100000111111101111111111110000000000 10 8.272282066785E-02 00000010000000111110111111111111110000000000 11 7.573216559561E-02 00001000000000111011111111111111110000000000 12 5.883063271692E-02 00100000000000101111111111111111110000000000 13 2.861976023692E-01 00000000000010111111111111011111110000000000 14 2.629756365689E-01 00000000001000111111111101111111110000000000 15 2.396195327431E-01 00000000100000111111110111111111110000000000 16 2.161284125517E-01 00000010000000111111011111111111110000000000 17 1.925013907583E-01 00001000000000111101111111111111110000000000 18 1.687375751663E-01 00100000000000110111111111111111110000000000 19 -1.437310483771E-01 00000000000001111111111111101111110000000000 20 -1.675216764579E-01 00000000000100111111111110111111110000000000 21 -1.914478021710E-01 00000000010000111111111011111111110000000000 22 -2.155103078794E-01 00000001000000111111101111111111110000000000 23 -2.397100828502E-01 00000100000000111110111111111111110000000000 24 -2.640480233188E-01 00010000000000111011111111111111110000000000 25 -2.885250325532E-01 01000000000000101111111111111111110000000000 26 -5.023566424930E-01 00000000000000111111111111111111110000000000 27 1.448360665541E-01 10000000000000011111111111111111110000000000 # 1 -1.8044 0.1911 0.1911 -0.0082 0.0000 -0.0079 -0.9124 0.0000 0.0000 0.0000 0.7476 13.2524 0.7476 0.2524 2 -1.3223 1.7545 11.7545 -2.8567 12.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 3 -1.3223 1.8567 6.8598 -2.8567 7.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 4 -1.3223 1.1433 12.5732 -2.8567 12.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 5 -1.3223 1.8567 18.2866 -2.8567 15.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 6 -1.3223 1.1433 18.2866 -2.8567 18.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 7 -1.3223 1.8567 24.0000 -2.8567 19.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 8 -1.3223 1.1433 41.1402 -2.8567 42.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 9 -1.3223 1.1433 35.4268 -2.8567 36.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 10 -1.3223 1.8567 29.7134 -2.8567 23.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 11 -1.3223 1.1433 24.0000 -2.8567 24.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 12 -1.3223 1.1433 29.7134 -2.8567 30.0000 -2.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 13 -1.3220 1.3306 13.1170 -2.4719 9.1633 -2.1634 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 14 -1.3214 1.5912 19.5912 -1.7136 20.0000 -1.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 15 -1.3214 1.7136 13.7187 -1.7136 14.0000 -1.5000 -2.0000 0.0000 0.0000 0.0000 1.0000 13.0000 1.0000 0.0000 Then I defined the Hamiltonian for final state 4fLf2p^55d^1 OppUdfG1 = NewOperator("U", NF, IndexUp_5d, IndexDn_5d, IndexUp_4f, IndexDn_4f, {0,0,0}, {1,0,0}) --only considered G1_5d4f -- soc in core 2p Oppcldots=NewOperator("ldots",NF, IndexUp_2p, IndexDn_2p) -- Coulomb interaction between 2p and 4f, omitted multipole part OppUpfF0 = NewOperator("U", NF, IndexUp_2p, IndexDn_2p, IndexUp_4f, IndexDn_4f, {1,0}, {0,0}) OppSz_5d = NewOperator("Sz", NF, IndexUp_5d,IndexDn_5d) OppLz_5d = NewOperator("Lz", NF, IndexUp_5d,IndexDn_5d) zeta_2p = 50 G1df = 1.0 -- taken from Kotani PRB 78 195115 2008 Upf = 10.5 -- taken from Kotani PRB 78 195115 2008 e4ffinal = Delta_4f_Lf_i - 6*Upf H_f = Chop(F0_4f_4f_i *OppF0_4f + F2_4f_4f_i *OppF2_4f+ F4_4f_4f_i *OppF4_4f + F6_4f_4f_i *OppF6_4f + zeta_4f_i * Oppldots_4f + zeta_Lf_i * Oppldots_Lf + Bz*(2*OppSz_4f + OppLz_4f) + Bz*(2*OppSz_5d + OppLz_5d) + zeta_2p * Oppcldots) + Chop(G1df * OppUdfG1 + Upf * OppUpfF0) + Chop(e4ffinal * OppN_4f) -- setting up the transition operators with various polarisations for 2p-->5d excitation -- Dipole X t=math.sqrt(1/2); Akm = {{1,-1,t},{1, 1,-t}} ; TXASx = NewOperator("CF", NF, IndexUp_5d, IndexDn_5d, IndexUp_2p, IndexDn_2p, Akm) --print(OperatorToMatrix(TXASx)) -- Dipole Y Akm = {{1,-1,t*I},{1, 1,t*I}} ; TXASy = NewOperator("CF", NF, IndexUp_5d, IndexDn_5d, IndexUp_2p, IndexDn_2p, Akm) -- Dipole Z Akm = {{1,0,1}} ; TXASz = NewOperator("CF", NF, IndexUp_5d, IndexDn_5d, IndexUp_2p, IndexDn_2p, Akm) --- right circular polarisation TXASr = t*(TXASx - I * TXASy) --- left circular polarisation TXASl =-t*(TXASx + I * TXASy) --get spectrum Spectra_z = CreateSpectra(H_f, TXASz, psiList[1], {{"Emin",-90}, {"Emax",100}, {"NE",5000}, {"Gamma",3.0}}) Spectra_r = CreateSpectra(H_f, TXASr, psiList[1], {{"Emin",-90}, {"Emax",100}, {"NE",5000}, {"Gamma",3.0}}) Spectra_l = CreateSpectra(H_f, TXASl, psiList[1], {{"Emin",-90}, {"Emax",100}, {"NE",5000}, {"Gamma",3.0}}) Spectra_ave = (Spectra_z + Spectra_r + Spectra_l)/3 Spectra_XMCD = Spectra_r - Spectra_l Following the model given by A.Kotani, in order to get proper XMCD, one should rescale the transition operator based on the exchange between 4f and 5d, i.e., the 4f spin for Ce can be fixed according to Hund's rule, and the transition from 2p to 5d thus has preference due to the 4f and 5d exchange. They introduced an enhance parameter alpha and renormalize the transition operator by sqrt(1-alpha*E_df_ex) to consider such effect, where E_df_ex is the exchange interaction energy between 4f and 5d states. As I understand from the Quanty output, the first eigenstate is a combination of several determinants. I want to find a way to easily reuse every determinant, for example, if I define wf1 = NewWavefunction(NF, NB, {{"00000000000100111111111111011111110000000000", 0.05799}}) so that I can first get the corresponding TXASr|wf1>. And then the new wave function TXASr|wf1> will be used to get E_df_ex by using TXASr|wf1> * OppUdfG1 * TXASr|wf1>. With the obtained E_df_ex, I can then renormalize the transition intensity for each single wave function accordingly. So my question is, is there any way to call every determinant in psiList[1] as shown in the output above and reuse them in the script ? Or is this a correct way to renormalize the transition intensity? I am very appreciated if you could tell me your opinions on this regard. Thanks a lot for your help! Best Regards, Ruiwen ~~DISCUSSION|Answers~~