-- read the radial wave functions -- they are stored in the file "NiO_Radial" -- order of functions -- r 1S 2S 2P 3S 3P 3D file = io.open("NiO_Radial", "r") Rnl = {} for line in file:lines() do RnlLine={} for i in string.gmatch(line, "%S+") do table.insert(RnlLine,i) end table.insert(Rnl,RnlLine) end -- some constants a0 = 0.52917721092 Rydberg = 13.60569253 Hartree = 2*Rydberg -- dd transitions from 3d (index 7 in Rnl) to -- 3d (index 7 in Rnl) -- define a function that calculates -- function RjRdd (q) Rj0R = 0 Rj2R = 0 Rj4R = 0 dr = Rnl[3][1]-Rnl[2][1] r0 = Rnl[2][1]-2*dr for ir = 2, #Rnl, 1 do r = r0 + ir * dr Rj0R = Rj0R + Rnl[ir][7] * math.SphericalBesselJ(0,q*r) * Rnl[ir][7] * dr Rj2R = Rj2R + Rnl[ir][7] * math.SphericalBesselJ(2,q*r) * Rnl[ir][7] * dr Rj4R = Rj4R + Rnl[ir][7] * math.SphericalBesselJ(4,q*r) * Rnl[ir][7] * dr end return Rj0R, Rj2R, Rj4R end -- the angular part is given as -- C(theta_q, phi_q)^* C(theta_r, phi_r) -- which is a potential expanded on -- spherical harmonics function ExpandOnClm(k,theta,phi,scale) ret={} for m=-k, k, 1 do table.insert(ret,{k,m,scale * math.SphericalHarmonicC(k,m,theta,phi)}) end return ret end -- define nIXS transition operators function TnIXS_dd(q, theta, phi) Rj0R, Rj2R, Rj4R = RjRdd(q) k=0 A0 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj0R) T0 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A0) k=2 A2 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj2R) T2 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A2) k=4 A4 = ExpandOnClm(k, theta, phi, I^k*(2*k+1)*Rj4R) T4 = NewOperator("CF", NF, IndexUp_3d, IndexDn_3d, A4) T = T0+T2+T4 T.Chop() return T end -- q in units per a0 (if you want in units per A -- take 5*a0 to have a q of 5 per A) q=4.5 print("for q=",q," per a0 (",q / a0," per A) The ratio of k=0, k=2 and k=4 transition strength is:",RjRdd(q)) -- define some transition operators qtheta=0 qphi=0 T = TnIXS_dd(q,qtheta,qphi) -- calculate the spectra nIXSSpectra = CreateSpectra(H, T, psi)