# 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:

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

$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.

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 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:

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

For completeness, the actual code is:

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")