Reuse every single determinant constituting the eigenstate
asked by 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(" # <E> <S^2> <L^2> <L_z^4f> <J^2> <J_z> <l.s> <F[2]> <F[4]> <F[6]> <N^4f> <N^Lf> <Weight_f1> <Weight_f0> ");
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 # <E> <S^2> <L^2> <L_z^4f> <J^2> <J_z> <l.s> <F[2]> <F[4]> <F[6]> <N^4f> <N^Lf> <Weight_f1> <Weight_f0> 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
Answers
Yes there is
https://www.quanty.org/documentation/language_reference/objects/wavefunction/properties/det
However!!!! I highly doubt that this is what you need and want. Why not work with TXASr |psi_1> i.e. act on the full many body state that is a superposition of determinants?
Maurits
Hi Maurits,
Thanks so much for your quick reply! It is of really great help. I was stuck somewhere previously so all my mind was thinking about treating the determinants separately. I have tried your suggestion out and the result seems promising. But I may leave another comment here if I have any further questions!
Thanks again! Ruiwen