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