Differences
This shows you the differences between two versions of the page.
documentation:tutorials:nio_ligand_field:temperature [2016/10/09 15:25] – created Maurits W. Haverkort | documentation:tutorials:nio_ligand_field:temperature [2016/10/10 09:41] (current) – external edit 127.0.0.1 | ||
---|---|---|---|
Line 1: | Line 1: | ||
+ | {{indexmenu_n> | ||
+ | ====== Temperature ====== | ||
+ | ### | ||
+ | The effect of temperature can be added by calculating excited states and expectation values of excited states. The temperature dependent expectation value is then created using Boltzmann statistics. | ||
+ | ### | ||
+ | |||
+ | ### | ||
+ | A small example: | ||
+ | <code Quanty Temperature.Quanty> | ||
+ | -- Sofar we calculated eigenstates and expectation values (or spectra) of these | ||
+ | -- eigenstates. At 0 K one would measure the expectation value of the lowest eigenstate | ||
+ | -- at finite temperature one would measure an average over several states weighted by | ||
+ | -- Boltzmann statistics. In this example we calculate the temperature dependent | ||
+ | -- x-ray absorption spectra of NiO. (Ni L23 edge 2p to 3d) within the ligand-field | ||
+ | -- theory approximation | ||
+ | |||
+ | -- The first part is an exact copy of example 41 | ||
+ | |||
+ | 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, | ||
+ | IndexUp_Ld={17, | ||
+ | |||
+ | -- angular momentum operators on the d-shell | ||
+ | |||
+ | OppSx_3d | ||
+ | OppSy_3d | ||
+ | OppSz_3d | ||
+ | OppSsqr_3d =NewOperator(" | ||
+ | OppSplus_3d=NewOperator(" | ||
+ | OppSmin_3d =NewOperator(" | ||
+ | |||
+ | OppLx_3d | ||
+ | OppLy_3d | ||
+ | OppLz_3d | ||
+ | OppLsqr_3d =NewOperator(" | ||
+ | OppLplus_3d=NewOperator(" | ||
+ | OppLmin_3d =NewOperator(" | ||
+ | |||
+ | OppJx_3d | ||
+ | OppJy_3d | ||
+ | OppJz_3d | ||
+ | OppJsqr_3d =NewOperator(" | ||
+ | OppJplus_3d=NewOperator(" | ||
+ | OppJmin_3d =NewOperator(" | ||
+ | |||
+ | Oppldots_3d=NewOperator(" | ||
+ | |||
+ | -- Angular momentum operators on the Ligand shell | ||
+ | |||
+ | OppSx_Ld | ||
+ | OppSy_Ld | ||
+ | OppSz_Ld | ||
+ | OppSsqr_Ld =NewOperator(" | ||
+ | OppSplus_Ld=NewOperator(" | ||
+ | OppSmin_Ld =NewOperator(" | ||
+ | |||
+ | OppLx_Ld | ||
+ | OppLy_Ld | ||
+ | OppLz_Ld | ||
+ | OppLsqr_Ld =NewOperator(" | ||
+ | OppLplus_Ld=NewOperator(" | ||
+ | OppLmin_Ld =NewOperator(" | ||
+ | |||
+ | OppJx_Ld | ||
+ | OppJy_Ld | ||
+ | OppJz_Ld | ||
+ | OppJsqr_Ld =NewOperator(" | ||
+ | OppJplus_Ld=NewOperator(" | ||
+ | OppJmin_Ld =NewOperator(" | ||
+ | |||
+ | -- 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(" | ||
+ | OppF2_3d =NewOperator(" | ||
+ | OppF4_3d =NewOperator(" | ||
+ | |||
+ | -- define onsite energies - crystal field | ||
+ | -- Akm = {{k1, | ||
+ | |||
+ | Akm = PotentialExpandedOnClm(" | ||
+ | OpptenDq_3d = NewOperator(" | ||
+ | OpptenDq_Ld = NewOperator(" | ||
+ | |||
+ | Akm = PotentialExpandedOnClm(" | ||
+ | OppNeg_3d = NewOperator(" | ||
+ | OppNeg_Ld = NewOperator(" | ||
+ | Akm = PotentialExpandedOnClm(" | ||
+ | OppNt2g_3d = NewOperator(" | ||
+ | OppNt2g_Ld = NewOperator(" | ||
+ | |||
+ | OppNUp_2p = NewOperator(" | ||
+ | OppNDn_2p = NewOperator(" | ||
+ | OppN_2p = OppNUp_2p + OppNDn_2p | ||
+ | OppNUp_3d = NewOperator(" | ||
+ | OppNDn_3d = NewOperator(" | ||
+ | OppN_3d = OppNUp_3d + OppNDn_3d | ||
+ | OppNUp_Ld = NewOperator(" | ||
+ | OppNDn_Ld = NewOperator(" | ||
+ | OppN_Ld = OppNUp_Ld + OppNDn_Ld | ||
+ | |||
+ | -- define L-d interaction | ||
+ | |||
+ | Akm = PotentialExpandedOnClm(" | ||
+ | OppVeg | ||
+ | Akm = PotentialExpandedOnClm(" | ||
+ | OppVt2g = NewOperator(" | ||
+ | |||
+ | -- core valence interaction | ||
+ | |||
+ | Oppcldots= NewOperator(" | ||
+ | OppUpdF0 = NewOperator(" | ||
+ | OppUpdF2 = NewOperator(" | ||
+ | OppUpdG1 = NewOperator(" | ||
+ | OppUpdG3 = NewOperator(" | ||
+ | |||
+ | -- dipole transition | ||
+ | |||
+ | t=math.sqrt(1/ | ||
+ | |||
+ | Akm = {{1, | ||
+ | TXASx = NewOperator(" | ||
+ | Akm = {{1, | ||
+ | TXASy = NewOperator(" | ||
+ | Akm = {{1,0,1}} | ||
+ | TXASz = NewOperator(" | ||
+ | |||
+ | TXASr = t*(TXASx - I * TXASy) | ||
+ | TXASl =-t*(TXASx + I * TXASy) | ||
+ | |||
+ | -- 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 | ||
+ | -- 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) | ||
+ | -- 9 eL + (n+1) ed + (n+1)n | ||
+ | -- 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/ | ||
+ | -- eL = nd*((1+nd)*Udd/ | ||
+ | -- | ||
+ | -- For the final state we/they defined | ||
+ | -- The 2p^5 L^10 d^n+1 configuration has an energy 0 | ||
+ | -- The 2p^5 L^9 d^n+2 configuration has an energy Delta + Udd - Upd | ||
+ | -- The 2p^5 L^8 d^n+3 configuration has an energy 2*Delta + 3*Udd - 2*Upd | ||
+ | -- | ||
+ | -- If we relate this to the onsite energy of the p and d orbitals we find | ||
+ | -- 6 ep + 10 eL + n ed + n(n-1) | ||
+ | -- 6 ep + 9 eL + (n+1) ed + (n+1)n | ||
+ | -- 6 ep + 8 eL + (n+2) ed + (n+1)(n+2) Udd/2 + 6 (n+2) Upd == 2*Delta+Udd | ||
+ | -- 5 ep + 10 eL + (n+1) ed + (n+1)(n) | ||
+ | -- 5 ep + 9 eL + (n+2) ed + (n+2)(n+1) Udd/2 + 5 (n+2) Upd == Delta+Udd-Upd | ||
+ | -- 5 ep + 8 eL + (n+3) ed + (n+3)(n+2) Udd/2 + 5 (n+3) Upd == 2*Delta+3*Udd-2*Upd | ||
+ | -- 6 equations with 3 unknowns, but with interdependence yield: | ||
+ | -- epfinal = (10*Delta + (1+nd)*(nd*Udd/ | ||
+ | -- edfinal = (10*Delta - nd*(31+nd)*Udd/ | ||
+ | -- eLfinal = ((1+nd)*(nd*Udd/ | ||
+ | -- | ||
+ | -- | ||
+ | -- | ||
+ | -- note that ed-ep = Delta - nd * U and not Delta | ||
+ | -- note furthermore that ep and ed here are defined for the onsite energy if the system had | ||
+ | -- locally nd electrons in the d-shell. In DFT or Hartree Fock the d occupation is in the end not | ||
+ | -- nd and thus the onsite energy of the Kohn-Sham orbitals is not equal to ep and ed in model | ||
+ | -- calculations. | ||
+ | -- | ||
+ | -- note furthermore that ep and eL actually should be different for most systems. We happily ignore this fact | ||
+ | -- | ||
+ | -- We normally take U and Delta as experimentally determined parameters | ||
+ | |||
+ | -- number of electrons (formal valence) | ||
+ | nd = 8 | ||
+ | -- parameters from experiment (core level PES) | ||
+ | Udd | ||
+ | Upd | ||
+ | Delta | ||
+ | -- parameters obtained from DFT (PRB 85, 165113 (2012)) | ||
+ | F2dd = 11.14 | ||
+ | F4dd = 6.87 | ||
+ | F2pd = 6.67 | ||
+ | G1pd = 4.92 | ||
+ | G3pd = 2.80 | ||
+ | tenDq | ||
+ | tenDqL | ||
+ | Veg | ||
+ | Vt2g = 1.21 | ||
+ | zeta_3d = 0.081 | ||
+ | zeta_2p = 11.51 | ||
+ | Bz = 0.000001 | ||
+ | Hz = 0.120 | ||
+ | |||
+ | ed = (10*Delta-nd*(19+nd)*Udd/ | ||
+ | eL = nd*((1+nd)*Udd/ | ||
+ | |||
+ | epfinal = (10*Delta + (1+nd)*(nd*Udd/ | ||
+ | edfinal = (10*Delta - nd*(31+nd)*Udd/ | ||
+ | eLfinal = ((1+nd)*(nd*Udd/ | ||
+ | |||
+ | F0dd = Udd + (F2dd+F4dd) * 2/63 | ||
+ | F0pd = Upd + (1/15)*G1pd + (3/70)*G3pd | ||
+ | |||
+ | Hamiltonian = F0dd*OppF0_3d + F2dd*OppF2_3d + F4dd*OppF4_3d + zeta_3d*Oppldots_3d + Bz*(2*OppSz_3d + OppLz_3d) + Hz * OppSz_3d + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g + ed * OppN_3d + eL * OppN_Ld | ||
+ | | ||
+ | XASHamiltonian = F0dd*OppF0_3d + F2dd*OppF2_3d + F4dd*OppF4_3d + zeta_3d*Oppldots_3d + Bz*(2*OppSz_3d + OppLz_3d)+ Hz * OppSz_3d + tenDq*OpptenDq_3d + tenDqL*OpptenDq_Ld + Veg * OppVeg + Vt2g * OppVt2g + edfinal * OppN_3d + eLfinal * OppN_Ld + epfinal * OppN_2p + zeta_2p * Oppcldots + F0pd * OppUpdF0 + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3 | ||
+ | |||
+ | -- we now can create the lowest Npsi eigenstates: | ||
+ | Npsi=3 | ||
+ | -- in order to make sure we have a filling of 8 electrons we need to define some restrictions | ||
+ | StartRestrictions = {NF, NB, {" | ||
+ | |||
+ | psiList = Eigensystem(Hamiltonian, | ||
+ | oppList={Hamiltonian, | ||
+ | |||
+ | -- print of some expectation values | ||
+ | |||
+ | print(" | ||
+ | for i = 1,#psiList do | ||
+ | io.write(string.format(" | ||
+ | for j = 1,#oppList do | ||
+ | expectationvalue = Chop(psiList[i]*oppList[j]*psiList[i]) | ||
+ | io.write(string.format(" | ||
+ | end | ||
+ | io.write(" | ||
+ | end | ||
+ | |||
+ | -- We calculate the x-ray absorption spectra for z, right circular and left circular polarized light for the 3 lowest eigen-states. (9 spectra in total) | ||
+ | |||
+ | XASSpectra = CreateSpectra(XASHamiltonian, | ||
+ | |||
+ | -- and put some additional energy broadening on it (0.4 Gaussian and energy dependent lorenzian) | ||
+ | XASSpectra.Broaden(0.4, | ||
+ | |||
+ | -- and now we start to do things different in order to include temperature averaging | ||
+ | -- We have calculated three states. The energy of these states is: | ||
+ | Enlist = {} | ||
+ | for i=1,# | ||
+ | Enlist[i] = psiList[i] * Hamiltonian * psiList[i] | ||
+ | end | ||
+ | -- In order to calculate occupation numbers we need to calculate E^(-Energy/ | ||
+ | -- If the ground-state has an energy of -3.5 eV and we calculate this pre-factor at a temperature | ||
+ | -- of 10 Kelvin we get E^(3.5/ | ||
+ | -- not fit in the computer memory. The solution is simple, we need to make sure the lowest | ||
+ | -- energy is zero. | ||
+ | for i=2,# | ||
+ | Enlist[i] = Enlist[i] - Enlist[1] | ||
+ | end | ||
+ | Enlist[1]=0 | ||
+ | |||
+ | -- Besides the energy I would like to look at the magnetic moment, both spin and angular part | ||
+ | -- so also here we calculate the expectation values in a list | ||
+ | Szlist = {} | ||
+ | Lzlist = {} | ||
+ | for i=1,# | ||
+ | Szlist[i] = psiList[i] * OppSz_3d * psiList[i] | ||
+ | Lzlist[i] = psiList[i] * OppLz_3d * psiList[i] | ||
+ | end | ||
+ | |||
+ | -- We can now calcualte the temperature dependent expectation values: | ||
+ | p={} | ||
+ | print(" | ||
+ | for T=1,1000,10 do | ||
+ | Z=0 | ||
+ | SzT=0 | ||
+ | LzT=0 | ||
+ | ET=0 | ||
+ | for i=1,# | ||
+ | p[i] = exp(-Enlist[i]/ | ||
+ | Z = Z + p[i] | ||
+ | SzT = SzT + p[i] * Szlist[i] | ||
+ | LzT = LzT + p[i] * Lzlist[i] | ||
+ | ET = ET + p[i] * Enlist[i] | ||
+ | end | ||
+ | SzT = SzT / Z | ||
+ | LzT = LzT / Z | ||
+ | MzT = -LzT - 2*SzT | ||
+ | | ||
+ | io.write(string.format(" | ||
+ | io.write(string.format(" | ||
+ | io.write(string.format(" | ||
+ | io.write(string.format(" | ||
+ | io.write(string.format(" | ||
+ | end | ||
+ | |||
+ | -- Note that in order to see the magnetic phase transition you need to make Hz in the Hamiltonian | ||
+ | -- temperature dependent. This can be done using self consistent loops. (The mathematica version | ||
+ | -- has an example of this, will make it here at some point as well) | ||
+ | |||
+ | -- The first column shows you how much energy you need to add to the system in order to heat | ||
+ | -- it. (specific heat) Be aware though that most of the specific heat is due to phonons, not | ||
+ | -- included in this calculation. It does capture nicely the electronic contribution to the | ||
+ | -- specific heat. (Including the change at cross overs for excited states (important in rare- | ||
+ | -- earth systems) and Lambda peaks for phase transitions | ||
+ | |||
+ | -- Now the temperature dependent XAS: | ||
+ | |||
+ | -- The object XASSpectra contains 9 spectra. For 3 different polarizations and 3 different states | ||
+ | -- What we need to do is to sum them according to Boltzmann statistics. | ||
+ | |||
+ | T = 100 | ||
+ | Z=0 | ||
+ | for i=1,# | ||
+ | p[i] = exp(-Enlist[i]/ | ||
+ | Z = Z + p[i] | ||
+ | end | ||
+ | -- we now create the Bolzmann summ for z, right and left polarized absorption (at T=100 Kelvin) | ||
+ | XASSpectraT100 = Spectra.Sum(XASSpectra, | ||
+ | |||
+ | -- and the same at 1000 Kelvin | ||
+ | T = 1000 | ||
+ | Z=0 | ||
+ | for i=1,# | ||
+ | p[i] = exp(-Enlist[i]/ | ||
+ | Z = Z + p[i] | ||
+ | end | ||
+ | -- we now create the Bolzmann summ for z, right and left polarized absorption (at T=1000 Kelvin) (note again that we still have a large magnetization here as the exchange-field is not temperature dependent) | ||
+ | XASSpectraT1000 = Spectra.Sum(XASSpectra, | ||
+ | |||
+ | |||
+ | -- We can print the spectra to file | ||
+ | XASSpectraT100.Print({{" | ||
+ | XASSpectraT1000.Print({{" | ||
+ | |||
+ | |||
+ | gnuplotInput = [[ | ||
+ | set autoscale | ||
+ | set xtic auto | ||
+ | set ytic auto | ||
+ | set style line 1 lt 1 lw 1 lc rgb "# | ||
+ | set style line 2 lt 1 lw 1 lc rgb "# | ||
+ | set style line 3 lt 1 lw 1 lc rgb "# | ||
+ | |||
+ | set xlabel "E (eV)" font " | ||
+ | set ylabel " | ||
+ | |||
+ | set out ' | ||
+ | set size 1.0, 0.6 | ||
+ | set terminal postscript portrait enhanced color " | ||
+ | |||
+ | energyshift=857.6 | ||
+ | intensityscale=48 | ||
+ | set xrange [847:877] | ||
+ | |||
+ | plot " | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | " | ||
+ | |||
+ | |||
+ | |||
+ | ]] | ||
+ | |||
+ | -- write the gnuplot script to a file | ||
+ | file = io.open(" | ||
+ | file: | ||
+ | file: | ||
+ | |||
+ | -- and finally call gnuplot to execute the script | ||
+ | os.execute(" | ||
+ | -- as I like pdf to view and eps to include in the manuel I transform the format | ||
+ | os.execute(" | ||
+ | </ | ||
+ | ### | ||
+ | |||
+ | ### | ||
+ | The output is: | ||
+ | <file Quanty_Output Temperature.out> | ||
+ | # < | ||
+ | 1 | ||
+ | 2 | ||
+ | 3 | ||
+ | Temperature, | ||
+ | | ||
+ | 11 | ||
+ | 21 | ||
+ | 31 | ||
+ | 41 | ||
+ | 51 | ||
+ | 61 | ||
+ | 71 | ||
+ | 81 | ||
+ | 91 | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | | ||
+ | </ | ||
+ | ### | ||
+ | |||
+ | |||
+ | ===== Table of contents ===== | ||
+ | {{indexmenu> |