# 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

, 2021/11/11 14:52

Yes there is

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

, 2021/11/11 19:22

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