Reuse every single determinant constituting the eigenstate

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 time!

Best Regards, Ruiwen