This is an old revision of the document!
Error while Calculating RIXS
asked by Galo Paez Fajardo (2024/03/19 18:40)
Running Quanty on Windows I got this error when: “On entry to DSYEVDSafe minimumPrecisionMIA parameter number 5 had an illegal value” and “On entry to DSYEV Safe minimumPrecisionM parameter number 5 had an illegal value”
System Windows 11, 64-bit operating system, AMD Ryzen 9 4900HS with Radeon Graphics 3.00 GHz, 16.0 GB RAM QuantyWin64.exe summer 2022 version
Quanty input file Verbosity(0)
-- here we calculate the 2p to 3d x-ray absorption of NiO within the Ligand-field theory -- approximation. The first part of the script is very much the same as calculating -- the ground-state with the addition that we now also need a 2p core shell in the basis
-- from the previous example we know that within NiO there are 3 states close to each other -- and then there is an energy gap of about 1 eV. We thus only need to consider the 3 -- lowest states (Npsi=3 later on)
NF=26 NB=0 IndexDn_2p={ 0, 2, 4} IndexUp_2p={ 1, 3, 5} IndexDn_3d={ 6, 8,10,12,14} IndexUp_3d={ 7, 9,11,13,15} IndexDn_Ld={16,18,20,22,24} IndexUp_Ld={17,19,21,23,25}
-- angular momentum operators on the d-shell
OppSx_3d =NewOperator(“Sx” ,NF, IndexUp_3d, IndexDn_3d) OppSy_3d =NewOperator(“Sy” ,NF, IndexUp_3d, IndexDn_3d) OppSz_3d =NewOperator(“Sz” ,NF, IndexUp_3d, IndexDn_3d) OppSsqr_3d =NewOperator(“Ssqr” ,NF, IndexUp_3d, IndexDn_3d) OppSplus_3d=NewOperator(“Splus”,NF, IndexUp_3d, IndexDn_3d) OppSmin_3d =NewOperator(“Smin” ,NF, IndexUp_3d, IndexDn_3d)
OppLx_3d =NewOperator(“Lx” ,NF, IndexUp_3d, IndexDn_3d) OppLy_3d =NewOperator(“Ly” ,NF, IndexUp_3d, IndexDn_3d) OppLz_3d =NewOperator(“Lz” ,NF, IndexUp_3d, IndexDn_3d) OppLsqr_3d =NewOperator(“Lsqr” ,NF, IndexUp_3d, IndexDn_3d) OppLplus_3d=NewOperator(“Lplus”,NF, IndexUp_3d, IndexDn_3d) OppLmin_3d =NewOperator(“Lmin” ,NF, IndexUp_3d, IndexDn_3d)
OppJx_3d =NewOperator(“Jx” ,NF, IndexUp_3d, IndexDn_3d) OppJy_3d =NewOperator(“Jy” ,NF, IndexUp_3d, IndexDn_3d) OppJz_3d =NewOperator(“Jz” ,NF, IndexUp_3d, IndexDn_3d) OppJsqr_3d =NewOperator(“Jsqr” ,NF, IndexUp_3d, IndexDn_3d) OppJplus_3d=NewOperator(“Jplus”,NF, IndexUp_3d, IndexDn_3d) OppJmin_3d =NewOperator(“Jmin” ,NF, IndexUp_3d, IndexDn_3d)
Oppldots_3d=NewOperator(“ldots”,NF, IndexUp_3d, IndexDn_3d)
-- Angular momentum operators on the Ligand shell
OppSx_Ld =NewOperator(“Sx” ,NF, IndexUp_Ld, IndexDn_Ld) OppSy_Ld =NewOperator(“Sy” ,NF, IndexUp_Ld, IndexDn_Ld) OppSz_Ld =NewOperator(“Sz” ,NF, IndexUp_Ld, IndexDn_Ld) OppSsqr_Ld =NewOperator(“Ssqr” ,NF, IndexUp_Ld, IndexDn_Ld) OppSplus_Ld=NewOperator(“Splus”,NF, IndexUp_Ld, IndexDn_Ld) OppSmin_Ld =NewOperator(“Smin” ,NF, IndexUp_Ld, IndexDn_Ld)
OppLx_Ld =NewOperator(“Lx” ,NF, IndexUp_Ld, IndexDn_Ld) OppLy_Ld =NewOperator(“Ly” ,NF, IndexUp_Ld, IndexDn_Ld) OppLz_Ld =NewOperator(“Lz” ,NF, IndexUp_Ld, IndexDn_Ld) OppLsqr_Ld =NewOperator(“Lsqr” ,NF, IndexUp_Ld, IndexDn_Ld) OppLplus_Ld=NewOperator(“Lplus”,NF, IndexUp_Ld, IndexDn_Ld) OppLmin_Ld =NewOperator(“Lmin” ,NF, IndexUp_Ld, IndexDn_Ld)
OppJx_Ld =NewOperator(“Jx” ,NF, IndexUp_Ld, IndexDn_Ld) OppJy_Ld =NewOperator(“Jy” ,NF, IndexUp_Ld, IndexDn_Ld) OppJz_Ld =NewOperator(“Jz” ,NF, IndexUp_Ld, IndexDn_Ld) OppJsqr_Ld =NewOperator(“Jsqr” ,NF, IndexUp_Ld, IndexDn_Ld) OppJplus_Ld=NewOperator(“Jplus”,NF, IndexUp_Ld, IndexDn_Ld) OppJmin_Ld =NewOperator(“Jmin” ,NF, IndexUp_Ld, IndexDn_Ld)
-- total angular momentum OppSx = OppSx_3d + OppSx_Ld OppSy = OppSy_3d + OppSy_Ld OppSz = OppSz_3d + OppSz_Ld OppSsqr = OppSx * OppSx + OppSy * OppSy + OppSz * OppSz OppLx = OppLx_3d + OppLx_Ld OppLy = OppLy_3d + OppLy_Ld OppLz = OppLz_3d + OppLz_Ld OppLsqr = OppLx * OppLx + OppLy * OppLy + OppLz * OppLz OppJx = OppJx_3d + OppJx_Ld OppJy = OppJy_3d + OppJy_Ld OppJz = OppJz_3d + OppJz_Ld OppJsqr = OppJx * OppJx + OppJy * OppJy + OppJz * OppJz
-- define the coulomb operator -- we here define the part depending on F0 seperately from the part depending on F2 -- when summing we can put in the numerical values of the slater integrals
OppF0_3d =NewOperator(“U”, NF, IndexUp_3d, IndexDn_3d, {1,0,0}) OppF2_3d =NewOperator(“U”, NF, IndexUp_3d, IndexDn_3d, {0,1,0}) OppF4_3d =NewOperator(“U”, NF, IndexUp_3d, IndexDn_3d, {0,0,1})
-- define onsite energies - crystal field -- Akm = 2_c_1_-1_-_c_1_1_akm_1_-1_t_1_1_-t TXASx = NewOperator(“CF”, NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm) -- y polarized light is defined as y = Sin[phi]Sin[theta] = sqrt(1/2) I ( C_1^{(-1)} + C_1^{(1)}) Akm = 1_-1_t_i_1_1_t_i TXASy = NewOperator(“CF”, NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm) -- z polarized light is defined as z = Cos[theta] = C_1^{(0)} Akm = 1_0_1 TXASz = NewOperator(“CF”, NF, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
TXASr = t*(TXASx - I * TXASy) TXASl =-t*(TXASx + I * TXASy)
-- the 3d to 2p dipole transition is the conjugate transpose of the 2p to 3d dipole transition TXASx = ConjugateTranspose(TXASx) TXASy = ConjugateTranspose(TXASy) TXASz = ConjugateTranspose(TXASz) TXASl = ConjugateTranspose(TXASl) TXASr = ConjugateTranspose(TXASr)
-- We follow the energy definitions as introduced in the group of G.A. Sawatzky (Groningen) -- J. Zaanen, G.A. Sawatzky, and J.W. Allen PRL 55, 418 (1985) -- for parameters of specific materials see -- A.E. Bockquet et al. PRB 55, 1161 (1996) -- After some initial discussion the energies U and Delta refer to the center of a configuration -- The L^10 d^n configuration has an energy 0 -- The L^9 d^n+1 configuration has an energy Delta -- The L^8 d^n+2 configuration has an energy 2*Delta+Udd -- -- If we relate this to the onsite energy of the L and d orbitals we find -- 10 eL + n ed + n(n-1) U/2 == 0 -- 9 eL + (n+1) ed + (n+1)n U/2 == Delta -- 8 eL + (n+2) ed + (n+1)(n+2) U/2 == 2*Delta+U -- 3 equations with 2 unknowns, but with interdependence yield: -- ed = (10*Delta-nd*(19+nd)*U/2)/(10+nd) -- eL = nd*1) F2dd = 10.622 F4dd = 6.636 F2pd = 6.68 G1pd = 5.066 G3pd = 2.882 tenDq = 0.95 tenDqL = 1.44 Veg = 3.0 Vt2g = 1.74 zeta_3d = 0.091 zeta_2p = 11.21 Bz = 0.000001 --Hz = 0.120
ed = (10*Delta-nd*(19+nd)*Udd/2)/(10+nd) eL = nd*2) / (16+nd) edfinal = (10*Delta - nd*(31+nd)*Udd/2-90*Upd) / (16+nd) eLfinal = 3)
io.write(string.format("%3i ",i))
for j = 1,#oppList do
expectationvalue = Chop(psiList[i]*oppList[j]*psiList[i])
io.write(string.format("%8.3f ",expectationvalue))
end
io.write("\n")
end
-- and we calculate the RIXS spectra. Note that in order to calculate RIXS you need to
-- specify two Hamiltonians.
-- These calculations can be lengthy and one should think that for each incoming energy
-- (E1) one needs to do a new calculation. In this case there are thus 120 calculations
RIXSSpectra = CreateResonantSpectra(XASHamiltonian, Hamiltonian, TXASx, TXASy, psiList[1], emin1_-15_emax1_25_ne1_240_gamma1_0.5_emin2_-0.5_emax2_7.5_ne2_8000_gamma2_0.05)
RIXSSpectra.Print(file_ni3d8_1_rixs.dat)
</WRAP>
~~DISCUSSION|Answers~~
Answers
Dear Galo,
That's an interesting error. Something seems to go wrong with your binding to Lapack. I'm not sure where the error comes form, but will forward it.
Best wishes, Maurits
Dear Maurits,
Thanks for your quick response.
The problem rises because the dipole operators both dag and nondag has the same name. Thus when invoking them in the line
RIXSSpectra = CreateResonantSpectra(XASHamiltonian, Hamiltonian, TXASx, TXASy, psiList[1], {{"Emin1",-15}, {"Emax1",25}, {"NE1",240}, {"Gamma1",0.5}, {"Emin2",-0.5}, {"Emax2",7.5}, {"NE2",8000}, {"Gamma2",0.05}})the software gets confused.
I changed the conjugate transpose operator names to
and now the RIXS line reads
RIXSSpectra = CreateResonantSpectra(XASHamiltonian, Hamiltonian, TXASx, TXASydag, psiList[1], {{"Emin1",-10}, {"Emax1",20}, {"NE1",240}, {"Gamma1",0.5}, {"Emin2",-0.5}, {"Emax2",7.5}, {"NE2",8000}, {"Gamma2",0.05}})Problem solved.
Thanks again. You guys have a fantastic piece of tool that I use a lot in data interpretation.
Sincerely, Galo
Dear Galo,
thanks for finding the error.
In principle your script should have given a zero spectrum, and not this error message, but then at least I do understand why others did not see this error.
Best wishes, Maurits