Differences

This shows you the differences between two versions of the page.

Link to this comparison view

Next revision		Previous revision
documentation:tutorials:small_programs_a_quick_start:xas [2016/10/07 19:54] – created Maurits W. Haverkortdocumentation:tutorials:small_programs_a_quick_start:xas [2017/01/15 15:48] (current) – Marius Retegan
Line 1: Line 1:
 +{{indexmenu_n>3}}
 +====== XAS ======
  
 +###
 +Example three calculates the $2p$ to $3d$ x-ray absorption spectra of transition metal compounds in cubic symmetry. It does so on a crystal field theory level. (Not the most accurate, but very informative. Ligand field calculations can be found later in the tutorials) As input it requires the ion name, the scaling factor for the Slater integrals, the size of the crystal field, a possible magnetic field and a temperature. The output contains the x-ray absorption spectra for all possible polarizations. Besides plots of a few standard polarizations the conductivity tensor is given, which allows one to create any experimental geometry possible.
 +###
 +
 +###
 +This is a nice little script to play with and learn what x-ray absorption spectroscopy measures. For example have a look how the linear and circular dichroism change in NiO (Ni$^{2+}$ with 10Dq=1.1, red=0.9 or so) as a function of temperature. You can see these in the conductivity tensor.
 +###
 +
 +###
 +The configuration file is:
 +<code Quanty conf.in>
 +-- 3d ion (e.g. Sc ... Cu) :
 +ion = "Ni"
 +-- Universal scaling parameter for the Slater integrals (beta):
 +beta = 0.90
 +-- Cubic Crystal Field Strength (10Dq) [eV]
 +tenDq = 1.1
 +-- External magnetic field Bext [Tesla]
 +Bx = 5 ; By = 0 ; Bz = 10
 +-- Temperature [Kelvin]
 +T = 1
 +</code>
 +###
 +
 +###
 +{{:documentation:tutorials:small_programs_a_quick_start:xasspec.png?nolink}}
 +
 +$L_{2,3}$ x-ray absorption spectra of Ni$^{2+}$ including magnetic circular dichroism.
 +###
 +
 +###
 +Besides some text this program also creates plots of the x-ray absorption spectra. In the configuration file the magnetic field is set in the $(102)$ direction. The temperature is $1$ Kelvin and as the field size in the $z$ direction is 10 Tesla the Ni ion will be fully magnetically polarized. One can calculate the magnetic circular dichroic spectra by calculating the spectra for left and right circular polarized light and calculating the difference. This is done in this example and in the figure above one can see in the top pannel the XAS for $z$, $R$ (right circular) and $L$ (left circular) polarized light. The bottom pannel shows the isotropic spectra and the magnetic circular dichroic spectrum.
 +###
 +
 +###
 +{{:documentation:tutorials:small_programs_a_quick_start:sigmatensor.png?nolink}}
 +
 +The conductivity tensor for excitations from $2p$ to $3d$ in Ni$^{2+}$
 +###
 +
 +###
 +In principle one can calculate the spectra for any magnetic field direction, but also for any polarization direction. The Poyinting vector must not be in the $z$ direction and there are thus infinite possibilities to define left circular polarized light. In optics this is solved by calculating the optical conductivity tensor. This is a three by three matrix that describes the optical properties of the material for any posible polarization. If $\sigma(\omega)$ is the energy dependent conductivity tensor (three by three matrix) and $\varepsilon$ the polarization (a vector of length three) then the absorption is give as: $I_{XAS} = -\mathrm{Im}[\varepsilon^* \cdot \sigma(\omega) \cdot \varepsilon$]. The conductivity tensor for Ni$^{2+}$ with a field in the $(102)$ direction is shown in the figure above. 
 +###
 +
 +###
 +The figure above shows several interesting features. The diagonal elements show magnetic linear dichroism, but also some of the off diagonal elements are influenced by magnetic linear dichroism. The $yz$ component for example is symmetric (equal to the $zy$ component) and also given by a magnetic linear dichroic effect. (See [[http://journals.aps.org/prb/abstract/10.1103/PhysRevB.82.094403|Phys. Rev. B 82 094402 (2010)]]). The magnetic circular dichroism due to the field in the $z$ direction can be found as the antisymmetric part of $xy$ and $yx$. The magnetic circular dichroism due to the field in the $x$ direction can be found as the antisymmetric part of $yz$ and $zy$.
 +###
 +
 +###
 +Information on the ground state is given in the text output, which is:
 +<file Quanty_Output XAS.out>
 +--> This program is the spectra calculator.
 +--> Assuming a specific ionic configuration for Me^{2+} , the L-edge spectra are produced and plotted in a separate file.
 +
 +--> Chosen parameters are as follows:
 +  Ion   beta       10Dq       Bext(x)      Bext(y)      Bext(z)       Temperature
 +  Ni    0.9000000  1.1000000  5.0000000  0.0000000   10.0000000       1.00 
 +
 +---> Setting up the Slater Integrals ...
 +---> Initial values are taken from M.W.Haverkort's PhD Thesis (p.156)
 +
 +
 +--> Assumed initial configuration: 2p6d8
 +--> INITIAL STATE PARAMETERS
 +  F0dd    F2dd    F4dd    zeta_3d
 + 0.5666 11.0097  6.8373  0.0830 
 +
 +--> Assumed final configuration: 2p5d9
 +--> FINAL STATE PARAMETERS
 +  F0dd    F2dd    F4dd    zeta_3d zeta_2p F2pd    G1pd    G3pd 
 + 0.6025 11.7045  7.2756  0.1020 11.5070  6.9480  5.2047  2.9610 
 +
 +--> Calculating the lowest 45 eigenstates...
 +
 +--> Ground state properties
 +
 +  <E>    <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <l.s>  <F[2]> <F[4]> <Neg> <Nt2g>
 +-2.706  1.999 12.000 15.154 -0.889 -0.266 -0.334 -1.020 -0.878  2.011  5.989 
 +-2.705  1.999 12.000 15.146 -0.000 -0.004 -0.332 -1.020 -0.878  2.011  5.989 
 +-2.703  1.999 12.000 15.138  0.889  0.257 -0.330 -1.020 -0.878  2.011  5.989 
 +-1.640  1.988 11.988 16.790 -0.167 -0.152 -0.873 -1.019 -0.877  3.008  4.992 
 +-1.640  1.988 11.988 16.795 -0.012  0.001 -0.875 -1.019 -0.877  3.008  4.992 
 +-1.621  1.999 12.000 15.825 -0.442 -0.029 -0.484 -1.020 -0.878  3.007  4.993 
 +-1.621  1.999 12.000 15.835  0.157  0.103 -0.487 -1.020 -0.878  3.007  4.993 
 +-1.620  1.999 12.000 15.840  0.441  0.078 -0.488 -1.020 -0.878  3.008  4.992 
 +-1.570  1.997 11.989 12.550 -0.436 -0.134  0.265 -1.020 -0.877  2.996  5.004 
 +-1.569  1.997 11.989 12.553 -0.058  0.005  0.265 -1.020 -0.877  2.996  5.004 
 +-1.568  1.997 11.989 12.557  0.450  0.122  0.264 -1.020 -0.877  2.996  5.004 
 +-1.549  2.000 12.000 12.001  0.068 -0.014  0.500 -1.020 -0.878  3.000  5.000 
 +-0.944  1.995 11.483 18.926 -0.035 -0.036 -1.474 -1.003 -0.887  3.584  4.416 
 +-0.876  1.999 11.421 15.293 -0.448  0.776 -0.506 -1.002 -0.887  3.568  4.432 
 +-0.876  1.999 11.421 15.303  0.035  0.014 -0.509 -1.002 -0.887  3.569  4.431 
 +-0.876  1.999 11.422 15.309  0.419 -0.769 -0.510 -1.002 -0.887  3.569  4.431 
 +-0.797  1.726 10.766  9.049 -0.012  0.023 -0.081 -0.981 -0.878  3.334  4.666 
 +-0.797  1.726 10.766  9.050  0.011  0.010 -0.081 -0.981 -0.878  3.334  4.666 
 +-0.766  1.997 11.279  9.961 -0.427  0.434  0.788 -0.998 -0.890  3.544  4.456 
 +-0.766  1.997 11.280  9.957  0.006  0.001  0.789 -0.998 -0.890  3.544  4.456 
 +-0.766  1.997 11.280  9.952  0.450 -0.423  0.790 -0.998 -0.890  3.544  4.456 
 +-0.489  0.290  9.837  9.162  0.000 -0.004  1.144 -0.902 -0.808  2.295  5.705 
 +-0.489  0.290  9.837  9.162  0.000 -0.001  1.144 -0.902 -0.808  2.295  5.705 
 + 0.459  0.151  7.374  7.664 -0.032 -0.490 -0.716 -0.887 -0.805  3.132  4.868 
 + 0.460  0.150  7.379  7.668 -0.000  0.003 -0.715 -0.887 -0.805  3.132  4.868 
 + 0.460  0.150  7.384  7.672  0.032  0.489 -0.714 -0.887 -0.805  3.132  4.868 
 + 0.804  0.030 18.827 18.853 -0.000 -0.016  0.042 -0.767 -0.843  2.420  5.580 
 + 0.923  1.992  2.783  6.359 -0.098 -0.081 -0.304 -0.737 -1.034  3.446  4.554 
 + 0.923  1.992  2.782  6.358 -0.008  0.005 -0.306 -0.737 -1.034  3.446  4.554 
 + 0.953  1.853  3.252  6.569 -0.417 -0.343  0.404 -0.749 -1.017  3.429  4.571 
 + 0.953  1.853  3.250  6.569  0.085  0.078  0.402 -0.749 -1.017  3.429  4.571 
 + 0.954  1.853  3.250  6.569  0.414  0.364  0.402 -0.749 -1.017  3.428  4.572 
 + 0.987  1.994  2.636  2.910 -0.437 -0.386  0.409 -0.732 -1.037  3.421  4.579 
 + 0.988  1.994  2.636  2.909 -0.020  0.017  0.409 -0.732 -1.037  3.420  4.580 
 + 0.989  1.994  2.636  2.909  0.454  0.372  0.409 -0.732 -1.037  3.420  4.580 
 + 1.019  1.972  2.755  1.287  0.028 -0.024  0.736 -0.731 -1.036  3.402  4.598 
 + 1.264  0.009 19.943 19.958 -0.002 -0.448  0.084 -0.776 -0.855  3.003  4.997 
 + 1.264  0.009 19.943 19.958  0.000  0.008  0.084 -0.776 -0.855  3.003  4.997 
 + 1.265  0.009 19.943 19.958  0.002  0.444  0.085 -0.776 -0.855  3.003  4.997 
 + 2.030  0.004 16.629 16.635 -0.000 -0.060  0.119 -0.810 -0.836  3.917  4.083 
 + 2.030  0.004 16.628 16.635  0.000  0.010  0.119 -0.810 -0.836  3.917  4.083 
 + 2.104  0.004 18.097 18.103  0.000 -1.807  0.090 -0.795 -0.845  3.888  4.112 
 + 2.105  0.004 18.101 18.107  0.000  0.054  0.090 -0.795 -0.845  3.888  4.112 
 + 2.106  0.004 18.106 18.112 -0.000  1.820  0.090 -0.795 -0.845  3.888  4.112 
 + 5.913  0.003  0.936  0.933  0.000  0.000  0.196 -0.581 -0.585  3.594  4.406 
 +
 +Now the finite temperature properties will be calculated.
 +Boltzmann averaging of the M=2Ms+Ml operator and of the X-ray absorption spectra.
 +
 +Partition function is sufficiently converged after Npsi= 2
 +Z,dZ 1.0000000364703 1.3291678830223e-15
 +
 +For a magnetic field of ( 5.000000 0.000000 10.000000 ) Tesla
 +At temperature T=1.000000 Kelvin, the magnetic moment is:
 +     Spin      Orbital   Total
 +x    0.44461   0.13279   1.02202
 +y    0.00000   0.00000   0.00000
 +z    0.88920   0.26553   2.04393
 +Printing the components of the conductivity tensor (sigma[i,j])
 +Printing the XAS spectra
 +Prepare gnuplot-file for XAS spectra
 +Prepare gnuplot-file for Sigma
 +
 +Execute the gnuplot to produce plots and convert the output into a pdf-file
 +FINISHED
 +</file>
 +###
 +
 +###
 +For completeness, the actual code is:
 +<code Quanty XAS.Quanty>
 +-- script contributed by Yaroslav Kvashnin
 +
 +-- Minimize the output of the program, i.e. for longer calculations
 +-- the program tells the user what it is doing. For these short
 +-- examples the result is instant.
 +Verbosity(0)
 +
 +print("--> This program is the spectra calculator.")
 +print("--> Assuming a specific ionic configuration for Me^{2+} , the L-edge spectra are produced and plotted in a separate file.")
 +print("")
 +print("--> Chosen parameters are as follows:")
 +
 +dofile("conf.in")
 +
 +io.write("  Ion   beta       10Dq       Bext(x)      Bext(y)      Bext(z)       Temperature","\n")
 +io.write(string.format("  %s    %2.7f  %2.7f  %2.7f  %2.7f   %2.7f       %2.2f ",ion,beta,tenDq,Bx,By,Bz,T))
 +io.write("\n")
 +
 +print("")
 +print("---> Setting up the Slater Integrals ...")
 +print("---> Initial values are taken from M.W.Haverkort's PhD Thesis (p.156)")
 +print("")
 +
 +if ion=="Cu" then
 +Nelec=9 
 +-- initial state parameters
 +zeta_3d=0.102; F2dd=12.854; F4dd=7.980   
 +--  final  state parameters                
 +zeta_2p=13.498; F2pd=8.177; G1pd=6.169; G3pd=3.510  
 +--  final  state parameters   
 +Xzeta_3d=0.124; XF2dd=13.611; XF4dd=8.457                  
 +
 +elseif ion=="Ni" then
 +Nelec=8
 +zeta_3d=0.083; F2dd=12.233; F4dd=7.597
 +zeta_2p=11.507; F2pd=7.720; G1pd=5.783; G3pd=3.290
 +Xzeta_3d=0.102; XF2dd=13.005; XF4dd=8.084
 +
 +elseif ion=="Co" then
 +Nelec=7
 +zeta_3d=0.066; F2dd=11.604; F4dd=7.209
 +zeta_2p=9.748; F2pd=7.259; G1pd=5.394; G3pd=3.068
 +Xzeta_3d=0.083; XF2dd=12.395; XF4dd=7.707
 +
 +elseif ion=="Fe" then
 +Nelec=6
 +zeta_3d=0.052; F2dd=10.965; F4dd=6.815
 +zeta_2p=8.200; F2pd=6.792; G1pd=5.000; G3pd=2.843
 +Xzeta_3d=0.067; XF2dd=11.778; XF4dd=7.327
 +
 +elseif ion=="Mn" then
 +Nelec=5
 +zeta_3d=0.040; F2dd=10.315 ; F4dd=6.413
 +zeta_2p=6.846; F2pd=6.320  ; G1pd=4.603; G3pd=2.617
 +Xzeta_3d=0.053; XF2dd=11.154 ; XF4dd=6.942
 +
 +elseif ion=="Cr" then
 +Nelec=4
 +zeta_3d=0.030; F2dd=9.648; F4dd=6.001
 +zeta_2p=5.668; F2pd=5.840; G1pd=4.201; G3pd=2.387
 +Xzeta_3d=0.041; XF2dd=10.521; XF4dd=6.551
 +
 +elseif ion=="V" then
 +Nelec=3
 +zeta_3d=0.022; F2dd=8.961; F4dd=5.576
 +zeta_2p=4.650; F2pd=5.351; G1pd=3.792; G3pd=2.154
 +Xzeta_3d=0.031; XF2dd=9.875; XF4dd=6.152
 +
 +elseif ion=="Ti" then
 +Nelec=2
 +zeta_3d=0.016; F2dd=8.243; F4dd=5.132
 +zeta_2p=3.776; F2pd=4.849; G1pd=3.376; G3pd=1.917
 +Xzeta_3d=0.023; XF2dd=9.21; XF4dd=5.744
 +
 +elseif ion=="Sc" then
 +Nelec=1
 +zeta_3d=0.010; F2dd=0; F4dd=0
 +zeta_2p=3.032; F2pd=4.332; G1pd=2.950; G3pd=1.674
 +Xzeta_3d=0.017; XF2dd=8.530; XF4dd=5.321
 +
 +else 
 +print("Could not recognize the ion name...")
 +os.exit()
 +end
 +
 +-- scaling with beta factor
 +F2dd=beta*F2dd; F4dd=beta*F4dd 
 +F2pd=beta*F2pd; G1pd=beta*G1pd; G3pd=beta*G3pd 
 +XF2dd=beta*XF2dd; XF4dd=beta*XF4dd
 +
 +-- number of possible many-body states in the initial configuration
 +Npsi = fact(10) / (fact(Nelec) * fact(10-Nelec))
 +
 +-- Bringing everything to the same units (eV)
 +T = T * EnergyUnits.Kelvin.value
 +Bx = Bx * EnergyUnits.Tesla.value
 +By = By * EnergyUnits.Tesla.value
 +Bz = Bz * EnergyUnits.Tesla.value
 +
 +-- 10 d-electrons + 6 p-electrons
 +NFermion=16                  
 +NBoson=0
 +-- p-shell [dn]
 +IndexDn_2p={0,2,4}
 +-- p-shell [up]        
 +IndexUp_2p={1,3,5}
 +-- d-shell [dn]         
 +IndexDn_3d={6,8,10,12,14}
 +-- d-shell [up]   
 +IndexUp_3d={7,9,11,13,15}   
 +
 +-- define the operators
 +OppSx   =NewOperator("Sx"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppSy   =NewOperator("Sy"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppSz   =NewOperator("Sz"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppSsqr =NewOperator("Ssqr" ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppSplus=NewOperator("Splus",NFermion, IndexUp_3d, IndexDn_3d)
 +OppSmin =NewOperator("Smin" ,NFermion, IndexUp_3d, IndexDn_3d)
 +
 +OppLx   =NewOperator("Lx"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppLy   =NewOperator("Ly"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppLz   =NewOperator("Lz"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppLsqr =NewOperator("Lsqr" ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppLplus=NewOperator("Lplus",NFermion, IndexUp_3d, IndexDn_3d)
 +OppLmin =NewOperator("Lmin" ,NFermion, IndexUp_3d, IndexDn_3d)
 +
 +OppJx   =NewOperator("Jx"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppJy   =NewOperator("Jy"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppJz   =NewOperator("Jz"   ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppJsqr =NewOperator("Jsqr" ,NFermion, IndexUp_3d, IndexDn_3d)
 +OppJplus=NewOperator("Jplus",NFermion, IndexUp_3d, IndexDn_3d)
 +OppJmin =NewOperator("Jmin" ,NFermion, IndexUp_3d, IndexDn_3d)
 +
 +Oppldots=NewOperator("ldots",NFermion, IndexUp_3d, IndexDn_3d)
 +
 +-- 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 =NewOperator("U", NFermion, IndexUp_3d, IndexDn_3d, {1,0,0})
 +OppF2 =NewOperator("U", NFermion, IndexUp_3d, IndexDn_3d, {0,1,0})
 +OppF4 =NewOperator("U", NFermion, IndexUp_3d, IndexDn_3d, {0,0,1})
 +
 +-- Crystal field operator for the d-shell
 +Akm = PotentialExpandedOnClm("Oh", 2, {0.6,-0.4})
 +OpptenDq = NewOperator("CF", NFermion, IndexUp_3d, IndexDn_3d, Akm)
 +Akm = PotentialExpandedOnClm("Oh",2,{1,0});
 +OppNeg = NewOperator("CF", NFermion, IndexUp_3d, IndexDn_3d, Akm)
 +Akm = PotentialExpandedOnClm("Oh",2,{0,1});
 +OppNt2g = NewOperator("CF", NFermion, IndexUp_3d, IndexDn_3d, Akm)
 +
 +-- Spin-orbit coupling in 2p-shell 
 +Oppcldots= NewOperator("ldots", NFermion, IndexUp_2p, IndexDn_2p)
 +-- core hole potentials:
 +-- direct 
 +OppUpdF0 = NewOperator("U", NFermion, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {1,0}, {0,0})
 +OppUpdF2 = NewOperator("U", NFermion, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,1}, {0,0})
 +-- exchange
 +OppUpdG1 = NewOperator("U", NFermion, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,0}, {1,0})
 +OppUpdG3 = NewOperator("U", NFermion, IndexUp_2p, IndexDn_2p, IndexUp_3d, IndexDn_3d, {0,0}, {0,1})
 +
 +t=math.sqrt(1/2);
 +
 +-- setting up the transition operators with various polarisations
 +-- Dipole X
 +Akm = {{1,-1,t},{1, 1,-t}} ;                
 +TXASx = NewOperator("CF", NFermion, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)
 +-- Dipole Y
 +Akm = {{1,-1,t*I},{1, 1,t*I}} ;             
 +TXASy = NewOperator("CF", NFermion, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)  
 +-- Dipole Z
 +Akm = {{1,0,1}} ;                           
 +TXASz = NewOperator("CF", NFermion, IndexUp_3d, IndexDn_3d, IndexUp_2p, IndexDn_2p, Akm)  
 +
 +TXASxpy  =  t*(TXASx + TXASy) 
 +TXASypz  =  t*(TXASy + TXASz) 
 +TXASzpx  =  t*(TXASz + TXASx) 
 +TXASxpiy =  t*(TXASx + I * TXASy) 
 +TXASypiz =  t*(TXASy + I * TXASz) 
 +TXASzpix =  t*(TXASz + I * TXASx) 
 +
 +--- right circular polarisation
 +TXASr = t*(TXASx - I * TXASy)
 +---  left circular polarisation
 +TXASl =-t*(TXASx + I * TXASy)
 +
 +------------------------------------------------
 +----- Input parameters for the Hamiltonian -----
 +-- in crystal field theory U drops out of the equation
 +U       =  0.000 
 +F0dd    = U+(F2dd+F4dd)*2/63
 +XF0dd   = U+(XF2dd+XF4dd)*2/63
 +-- in crystal field theory U drops out of the equation
 +Upd     =  0.000 
 +F0pd    =  Upd + G1pd*1/15 + G3pd*3/70
 +
 +print("")
 +io.write(string.format("--> Assumed initial configuration: 2p6d%i", Nelec))
 +llist = {F0dd, F2dd, F4dd, zeta_3d}
 +print("")
 +print("--> INITIAL STATE PARAMETERS")
 +io.write("  F0dd    F2dd    F4dd    zeta_3d\n")
 +for i=1,4 do
 +io.write(string.format("%7.4f ",llist[i]))
 +end
 +print("")
 +print("")
 +io.write(string.format("--> Assumed final configuration: 2p5d%i",Nelec+1))
 +llist2 = {XF0dd, XF2dd, XF4dd, Xzeta_3d, zeta_2p, F2pd, G1pd, G3pd}
 +print("")
 +print("--> FINAL STATE PARAMETERS")
 +io.write("  F0dd    F2dd    F4dd    zeta_3d zeta_2p F2pd    G1pd    G3pd \n")
 +for i=1,8 do
 +io.write(string.format("%7.4f ",llist2[i]))
 +end
 +print("")
 +
 +-- initial state Hamiltonian 
 +Hamiltonian =  F0dd*OppF0 + F2dd*OppF2 + F4dd*OppF4 + tenDq*OpptenDq + zeta_3d*Oppldots 
 +Hamiltonian = Hamiltonian + Bx * (2*OppSx + OppLx) + By * (2*OppSy + OppLy) + Bz * (2*OppSz + OppLz)
 +
 +-- final state Hamiltonian 
 +XASHamiltonian = XF0dd*OppF0 + XF2dd*OppF2 + XF4dd*OppF4 + tenDq*OpptenDq + Xzeta_3d*Oppldots + zeta_2p * Oppcldots + F2pd * OppUpdF2 + G1pd * OppUpdG1 + G3pd * OppUpdG3
 +XASHamiltonian = XASHamiltonian + Bx * (2*OppSx + OppLx) + By * (2*OppSy + OppLy) + Bz * (2*OppSz + OppLz)
 +
 +-- in order to make sure we have a filling of 2 electrons we need to define some restrictions
 +StartRestrictions = {NFermion, NBoson, {"111111 0000000000",6,6}, {"000000 1111111111",Nelec,Nelec}}
 +
 +print("")
 +print("--> Calculating the lowest", Npsi ,"eigenstates...")
 +psiList = Eigensystem(Hamiltonian, StartRestrictions, Npsi)
 +
 +
 +oppList={Hamiltonian, OppSsqr, OppLsqr, OppJsqr, OppSz, OppLz, Oppldots, OppF2, OppF4, OppNeg, OppNt2g};
 +
 +
 +print("")
 +print("--> Ground state properties")
 +print("")
 +
 +-- Energy of the lowest multiplet. To be shifted to 0 to get a proper partition function
 +
 +Egrd=psiList[1] * Hamiltonian * psiList[1]
 +
 +print("  <E>    <S^2>  <L^2>  <J^2>  <S_z>  <L_z>  <l.s>  <F[2]> <F[4]> <Neg> <Nt2g>");
 +for key,value in pairs(psiList) do
 +  expvalue = value * oppList * value 
 +  for k,v in pairs(expvalue) do
 +    io.write(string.format("%6.3f ",Complex.Re(v)))
 +  end
 +  io.write("\n")
 +end
 +
 +
 +print("")
 +print("Now the finite temperature properties will be calculated.") 
 +print("Boltzmann averaging of the M=2Ms+Ml operator and of the X-ray absorption spectra.")
 +print("")
 +
 +Z=0
 +Mx =0
 +MSx=0
 +MLx=0
 +My =0
 +MSy=0
 +MLy=0
 +Mz =0
 +MSz=0
 +MLz=0
 +Spectra_z=0
 +Spectra_r=0
 +Spectra_l=0
 +
 +--- Conductivity tensor
 +Sigma_X=0 
 +Sigma_Y=0 
 +Sigma_Z=0 
 +Sigma_X_p_Y=0 
 +Sigma_Y_p_Z=0 
 +Sigma_Z_p_X=0 
 +Sigma_X_p_IY=0 
 +Sigma_Y_p_IZ=0 
 +Sigma_Z_p_IX=0 
 +
 +Z = 0 ; 
 +
 + for j=1, Npsi do
 +  dZ = math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  if (dZ < 0.000000001) then
 +    print("Partition function is sufficiently converged after Npsi=",j-1)
 + print("Z,dZ",Z,dZ)
 + break
 +  end
 +  Z  = Z  + dZ
 +  Mx  = Mx  - psiList[j] * (2 * OppSx + OppLx) * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  MSx = MSx - psiList[j] * (OppSx)             * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  MLx = MLx - psiList[j] * (OppLx)             * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  My  = My  - psiList[j] * (2 * OppSy + OppLy) * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  MSy = MSy - psiList[j] * (OppSy)             * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  MLy = MLy - psiList[j] * (OppLy)             * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Mz  = Mz  - psiList[j] * (2 * OppSz + OppLz) * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  MSz = MSz - psiList[j] * (OppSz)             * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  MLz = MLz - psiList[j] * (OppLz)             * psiList[j] * math.exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Spectra_z = Spectra_z + CreateSpectra(XASHamiltonian, TXASz, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Spectra_r = Spectra_r + CreateSpectra(XASHamiltonian, TXASr, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Spectra_l = Spectra_l + CreateSpectra(XASHamiltonian, TXASl, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_X = Sigma_X + CreateSpectra(XASHamiltonian, TXASx, psiList[j],{{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_Y = Sigma_Y + CreateSpectra(XASHamiltonian, TXASy, psiList[j],{{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_Z = Sigma_Z + CreateSpectra(XASHamiltonian, TXASz, psiList[j],{{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_X_p_Y = Sigma_X_p_Y + CreateSpectra(XASHamiltonian, TXASxpy, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_Y_p_Z = Sigma_Y_p_Z + CreateSpectra(XASHamiltonian, TXASypz, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_Z_p_X = Sigma_Z_p_X + CreateSpectra(XASHamiltonian, TXASzpx, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_X_p_IY = Sigma_X_p_IY + CreateSpectra(XASHamiltonian, TXASxpiy, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_Y_p_IZ = Sigma_Y_p_IZ + CreateSpectra(XASHamiltonian, TXASypiz, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 +  Sigma_Z_p_IX = Sigma_Z_p_IX + CreateSpectra(XASHamiltonian, TXASzpix, psiList[j], {{"Emin",-10}, {"Emax",20}, {"NE",3500}, {"Gamma",1.0}}) * exp(-(psiList[j] * Hamiltonian * psiList[j] - Egrd)/T)
 + end;
 +Mx  = Mx  / Z
 +MSx = MSx / Z
 +MLx = MLx / Z
 +My  = My  / Z
 +MSy = MSy / Z
 +MLy = MLy / Z
 +Mz  = Mz  / Z
 +MSz = MSz / Z
 +MLz = MLz / Z
 +Spectra_z = Spectra_z/Z
 +Spectra_r = Spectra_r/Z
 +Spectra_l = Spectra_l/Z
 +Sigma_X = Sigma_X/Z
 +Sigma_Y = Sigma_Y/Z
 +Sigma_Z = Sigma_Z/Z
 +Sigma_X_p_Y = Sigma_X_p_Y/Z
 +Sigma_Y_p_Z = Sigma_Y_p_Z/Z
 +Sigma_Z_p_X = Sigma_Z_p_X/Z
 +Sigma_X_p_IY = Sigma_X_p_IY/Z
 +Sigma_Y_p_IZ = Sigma_Y_p_IZ/Z
 +Sigma_Z_p_IX = Sigma_Z_p_IX/Z
 +
 +print("")
 +
 +io.write(string.format("For a magnetic field of ( %f %f %f ) Tesla\n", Bx/EnergyUnits.Tesla.value, By/EnergyUnits.Tesla.value, Bz/EnergyUnits.Tesla.value))
 +io.write(string.format("At temperature T=%f Kelvin, the magnetic moment is:\n", T/EnergyUnits.Kelvin.value))
 +print("     Spin      Orbital   Total")
 +io.write(string.format("x  %9.5f %9.5f %9.5f\n", Complex.Re(MSx),Complex.Re(MLx),Complex.Re(Mx)))
 +io.write(string.format("y  %9.5f %9.5f %9.5f\n", Complex.Re(MSy),Complex.Re(MLy),Complex.Re(My)))
 +io.write(string.format("z  %9.5f %9.5f %9.5f\n", Complex.Re(MSz),Complex.Re(MLz),Complex.Re(Mz)))
 +
 +Spectra_ave  = (Spectra_z + Spectra_l + Spectra_r)/3
 +Spectra_XMCD = (Spectra_r - Spectra_l)
 +
 +Sigma_11 = Sigma_X
 +Sigma_22 = Sigma_Y
 +Sigma_33 = Sigma_Z
 +
 +Sigma_12 = Sigma_X_p_Y - I * Sigma_X_p_IY - 0.5 * (1 - I) * (Sigma_X + Sigma_Y)
 +Sigma_21 = Sigma_X_p_Y + I * Sigma_X_p_IY - 0.5 * (1 + I) * (Sigma_Y + Sigma_X) 
 +
 +Sigma_13 = Sigma_Z_p_X + I * Sigma_Z_p_IX - 0.5 * (1 + I) * (Sigma_X + Sigma_Z)
 +Sigma_31 = Sigma_Z_p_X - I * Sigma_Z_p_IX - 0.5 * (1 - I) * (Sigma_Z + Sigma_X)
 +
 +Sigma_23 = Sigma_Y_p_Z - I * Sigma_Y_p_IZ - 0.5 * (1 - I) * (Sigma_Y + Sigma_Z)
 +Sigma_32 = Sigma_Y_p_Z + I * Sigma_Y_p_IZ - 0.5 * (1 + I) * (Sigma_Z + Sigma_Y)
 +
 +
 +print("Printing the components of the conductivity tensor (sigma[i,j])")
 +Sigma_11.Print({{"file","Sigma11.dat"}})
 +Sigma_22.Print({{"file","Sigma22.dat"}})
 +Sigma_33.Print({{"file","Sigma33.dat"}})
 +Sigma_12.Print({{"file","Sigma12.dat"}})
 +Sigma_21.Print({{"file","Sigma21.dat"}})
 +Sigma_13.Print({{"file","Sigma13.dat"}})
 +Sigma_31.Print({{"file","Sigma31.dat"}})
 +Sigma_23.Print({{"file","Sigma23.dat"}})
 +Sigma_32.Print({{"file","Sigma32.dat"}})
 +
 +print("Printing the XAS spectra")
 +Spectra_z.Print({{"file","XASSpecZ.dat"}})
 +Spectra_r.Print({{"file","XASSpecR.dat"}})
 +Spectra_l.Print({{"file","XASSpecL.dat"}})
 +Spectra_ave.Print({{"file","XASSpecAVER.dat"}})
 +Spectra_XMCD.Print({{"file","XASSpecXMCD.dat"}})
 +
 +--  prepare the gnuplot output for XAS 
 +gnuplotInput = [[
 +set autoscale   # scale axes automatically
 +set xtic auto   # set xtics automatically
 +set ytic auto   # set ytics automatically
 +
 +set style line  1 lt 1 lw 2 lc rgb "#FF0000"
 +set style line  2 lt 1 lw 2 lc rgb "#00FF00"
 +set style line  3 lt 1 lw 2 lc rgb "#0000FF"
 +set style line  4 lt 1 lw 2 lc rgb "#000000"
 +
 +set xlabel "E (eV)" font "Times,12"
 +set ylabel "Intensity (arb. units)" font "Times,12"
 +
 +set out 'XASSpec.ps'
 +set size 1.0, 1.0
 +set terminal postscript portrait enhanced color "Times" 12
 +
 +set multiplot layout 4, 1
 +plot "XASSpecZ.dat"    using 1:(-$3) title 'z-polarized' with lines ls 1,\
 +     "XASSpecR.dat"    using 1:(-$3) title 'R-polarized' with lines ls 2,\
 +     "XASSpecL.dat"    using 1:(-$3) title 'L-polarized' with lines ls 3
 +plot "XASSpecAVER.dat" using 1:(-$3) title 'Average' with lines ls 4,\
 +     "XASSpecXMCD.dat" using 1:(-$3) title 'XMCD' with lines ls 1
 +unset multiplot
 +
 +]]
 +
 +print("Prepare gnuplot-file for XAS spectra")
 +-- write the gnuplot script to a file
 +file = io.open("XASSpec.gnuplot", "w")
 +file:write(gnuplotInput)
 +file:close()
 +
 +-- prepare the gnuplot output for Sigma 
 +gnuplotInput = [[
 +set autoscale   # scale axes automatically
 +set xtic auto   # set xtics automatically
 +set ytic auto   # set ytics automatically
 +set style line  1 lt 1 lw 2 lc 1
 +set style line  2 lt 1 lw 2 lc 3
 +
 +set xlabel "E (eV)" font "Times,10"
 +set ylabel "Intensity (arb. units)" font "Times,10"
 +
 +set yrange [-0.3:0.3]
 +
 +set out 'SigmaTensor.ps'
 +set size 1.0, 1.0
 +set terminal postscript portrait enhanced color  "Times" 8
 +
 +set multiplot layout 6, 3
 +
 +plot "Sigma11.dat" u 1:2 title 'Re[{/Symbol s}_{11}]' with lines ls 1,\
 +     "Sigma11.dat" u 1:3 title 'Im[{/Symbol s}_{11}]' with lines ls 2
 +plot "Sigma12.dat" u 1:2 title 'Re[{/Symbol s}_{12}]' with lines ls 1,\
 +     "Sigma12.dat" u 1:3 title 'Im[{/Symbol s}_{12}]' with lines ls 2
 +plot "Sigma13.dat" u 1:2 title 'Re[{/Symbol s}_{13}]' with lines ls 1,\
 +     "Sigma13.dat" u 1:3 title 'Im[{/Symbol s}_{13}]' with lines ls 2
 +plot "Sigma21.dat" u 1:2 title 'Re[{/Symbol s}_{21}]' with lines ls 1,\
 +     "Sigma21.dat" u 1:3 title 'Im[{/Symbol s}_{21}]' with lines ls 2
 +plot "Sigma22.dat" u 1:2 title 'Re[{/Symbol s}_{22}]' with lines ls 1,\
 +     "Sigma22.dat" u 1:3 title 'Im[{/Symbol s}_{22}]' with lines ls 2
 +plot "Sigma23.dat" u 1:2 title 'Re[{/Symbol s}_{23}]' with lines ls 1,\
 +     "Sigma23.dat" u 1:3 title 'Im[{/Symbol s}_{23}]' with lines ls 2
 +plot "Sigma31.dat" u 1:2 title 'Re[{/Symbol s}_{31}]' with lines ls 1,\
 +     "Sigma31.dat" u 1:3 title 'Im[{/Symbol s}_{31}]' with lines ls 2
 +plot "Sigma32.dat" u 1:2 title 'Re[{/Symbol s}_{32}]' with lines ls 1,\
 +     "Sigma32.dat" u 1:3 title 'Im[{/Symbol s}_{32}]' with lines ls 2
 +plot "Sigma33.dat" u 1:2 title 'Re[{/Symbol s}_{33}]' with lines ls 1,\
 +     "Sigma33.dat" u 1:3 title 'Im[{/Symbol s}_{33}]' with lines ls 2
 +
 +unset multiplot
 +]]
 +
 +print("Prepare gnuplot-file for Sigma")
 +
 +-- write the gnuplot script to a file
 +file = io.open("SigmaTensor.gnuplot", "w")
 +file:write(gnuplotInput)
 +file:close()
 +
 +print("")
 +print("Execute the gnuplot to produce plots and convert the output into a pdf-file")
 +
 +-- call gnuplot to execute the script
 +os.execute("gnuplot XASSpec.gnuplot ; ps2pdf XASSpec.ps ; ps2eps XASSpec.ps ; mv XASSpec.eps temp.eps ; eps2eps temp.eps XASSpec.eps ; rm temp.eps")
 +os.execute("gnuplot SigmaTensor.gnuplot ; ps2pdf SigmaTensor.ps ; ps2eps SigmaTensor.ps ;  mv SigmaTensor.eps temp.eps ; eps2eps temp.eps SigmaTensor.eps ; rm temp.eps")
 +
 +print("FINISHED")
 +</code>
 +###
 +
 +
 +===== Table of contents =====
 +{{indexmenu>.#1|msort}}

You are here: Quanty - a quantum many body script language » Documentation » Tutorials » Small programs: A quick start » XAS

Print/export