diff --git a/libhreels/README.md b/libhreels/README.md new file mode 100755 index 0000000000000000000000000000000000000000..1c6ae5c58f71c022695d132b76f816eb67301755 --- /dev/null +++ b/libhreels/README.md @@ -0,0 +1,58 @@ +# HREELS data handling and simulation + +Set of data analysis routines for surface vibrational spectroscopy based on a Delta05 spectrometer. Simple data reading and plotting routines are provided, but also advanced routines for in depth analysis. + +A fast data browser as graphical user interface is included as **ViewHREELS.py**. It can also take command line arguments. + +Additional a python interface to a Fortran-based calculation of full HREEL spectra is provided in **class lambin** (in **calcHREELS20.py**). Note that this part requires a local compilation (via f2py3) from the Fortran90 files. Only for the Linux-based python version 3.8, there are complied files provided. For details see below. + +Most of the routines are within the **class HREELS** (HREELS.py). A simple command line program is "showHREELS.py", which reads one data file and plots the spectrum with a second amplified trace. + + Usage: >python showHREELS.py filename [factor wavenumber] + + E.g.: python showHREELS.py H9H03 100 53 + +More general usage: + +Read dataset by calling the HREELS class: + **data1 = HREELS('filename', datapath='datapath')** # Here you can omit the extension '.gph' + # The second argument is optional + + This will assign/calculate the following properties: + + data1.filename + data1.datapath + data1.startTime + data1.stopTime + data1.totalTime + data1.timePerChannel + data1.numberOfSegments + data1.energy # Electron kinetic energy + data1.filament # Filament current in Ampere + data1.segments # list of segment info + data1.data # [(-100.1021616, 259.5), (-98.5374352, 264.5), ...] + data1.xdata + data1.ydata + data1.maxIntensity # Count rate of elastic peak + data1.resolution # width of elastic peak + +The following methods are defined within the class: + + info() : Print information about the spectrum + plot() : Draws the spectrum. Optional arguments are: + (xmin=None, xmax=None, factor=1, label='x', normalized=False, color="b-",marker=True) + + plotInfoAmp() Draws spectrum together with amplified trace. + pickPeak() : Select peak by mouse cursor. The call of figure() is required before. + selectData(): Returns the data between wavenumbers x1 and x2 + findIndex(lossenergy): Returns the data array index for a given energy loss + setMarker(x, y, ymin=0, size=None): Sets vertical marker with text label + plotWaterFall(...): + + +# calcHREELS20 + +These routines are based on the publication "Computation of the surface electron-energy-loss spectrum in specular geometry for an arbitrary plane-stratified medium" by P. Lambin, J.-P. Vigneron, and A. A. Lucas, in the Journal "Computer Physics Communications 60, 351-64(1990)". + + + diff --git a/libhreels/analysizeData.py b/libhreels/analysizeData.py deleted file mode 100644 index 5184c1ef186f696163dc21f2f6cb7b9467e04ffa..0000000000000000000000000000000000000000 --- a/libhreels/analysizeData.py +++ /dev/null @@ -1,45 +0,0 @@ -import numpy as np -from matplotlib import pyplot as plt - -def findPeriod(list, period): - '''For a list of times (relative to a sync pulse) calculate the average time deviation squared. - ''' - residuum = [] - sum_residuum = 0 - for each in list: - r = each % period - residuum.append(r) - sum_residuum += r - av_res = sum_residuum / len(list) - deviation = 0 - for res in residuum: - deviation += (res - av_res)**2 - deviation /= (len(residuum)*period**2) - return deviation, av_res - -def makeTimeRelative(times, durations): - list = [] - synced = False - t_zero = 0 - for each in pulses: - if synced: - list.append(pulses[each][0]-t_zero) - if 9000 < pulses[each][2] < 9800: - print("synced") - t_zero = pulses[each][0] - synced = True - return list - -data = [] # read dat from file [[time in microsec, low duration, high duration], [...], ...] -with open("data433_Signal2.txt", 'r') as f: - for line in f: - data.append([int(each) for each in line.split(',')]) - -sumOn = 0; sumOff = 0 -for pulse in data[1:]: - sumOn += pulse[2] - sumOff += pulse[1] -sumOn /= len(data) -sumOff /= len(data) - -print("Average on and off times:", sumOn, sumOff) diff --git a/libhreels/calcHREELS.py b/libhreels/calcHREELS.py deleted file mode 100644 index b3000e031383a2a427659344e567760b5a5c88b2..0000000000000000000000000000000000000000 --- a/libhreels/calcHREELS.py +++ /dev/null @@ -1,212 +0,0 @@ -#!/usr/bin/env python3 -import numpy as np -import json -import sys, re, os -from libhreels.HREELS import myPath -from copy import deepcopy - -libDir = os.path.dirname(os.path.realpath(__file__)) - -if not (sys.version.startswith('3.6') or sys.version.startswith('3.8')): - print('''Make sure the Fortran routines 'myEels2' and 'myBoson' - have been complied with the proper f2py3 for the right - python version!!''') -try: - from libhreels import myEels2 as LambinEELS # wrapper for myEels2.f90 - from libhreels import myBoson as LambinBoson # wrapper for myBoson.f90 - from libhreels import myEels3 as LambinEELS3 # wrapper for myEels3.f90 -except: - print('myEels2 and MyBoson are not available here (Check your version)') - -# Experimental setup as dictionary: -setup = { - "e0": 4.0, - "theta": 60., - "phia": 0.33, - "phib": 2.0, - "temperature": 298., - "debug": False -} -# Instrumental function describing elastic peak shape: -instrument = { - "width": 18., - "intensity": 100000., - "asym": 0.01, - "gauss": 0.88 -} - - -def importMaterials(string='', path=libDir): - ''' Returns a dictionary with all phonon parameters for the material provided - as string argument. If the string is empty or not matching, a list of all available - materials is printed. - ''' - file = os.path.join(myPath(path),'materials.json') - with open(file) as json_file: - materials = json.load(json_file) - try: - mat = materials[string] - except: - print('No data for material >>{}<< found in {} materials.json!!'.format(string, path)) - print('Available materials:\n{}\n'.format(materials.keys())) - mat = 'None' - return mat - - -def addDrude(wPlasma, gPlasma, material): - ''' Adds a simple Drude response to the materials properties (which are provided - as 3rd argument) and returns a new materials dictionary with all phonon parameters. Note - that at least the eps_infinity has to given before. - ''' - newMaterial = deepcopy(material) - # To select the Drude response, the Lambin Fortran code requires a negative first oscillator frequency - # which is then interpreted as plasma frequency. - # The values of the former LO parameter are irrelevant (but need to be provided). - if not wPlasma > 0: - print('Cannot add Drude to material',material) - return material - try: - if len(newMaterial['wTO']) > 0: - newMaterial['wTO'] += [-1*wPlasma] - newMaterial['gTO'] += [-1*gPlasma] - newMaterial['wLO'] += [0] - newMaterial['gLO'] += [0] - return newMaterial - except: - print('Cannot add Drude to material',material) - return material - -################################################################################ -################################################################################ -class lambin: - def __init__(self, film, setup=setup, instrument=instrument): - self.e0 = setup['e0'] - self.theta = setup['theta'] - self.phia = setup['phia'] - self.phib = setup['phib'] - self.temperature = setup['temperature'] - self.debug = setup['debug'] - self.width = instrument['width'] - self.gauss = instrument['gauss'] - self.intensity = instrument['intensity'] - self.asym = instrument['asym'] - self.layers = len(film) # number of layers - self.neps = self.layers - name_size = self.layers - self.name = []; self.thick=[]; self.listNOsci=[]; self.epsinf =[]; Q = [] - allTO=[]; allgTO=[]; allgLO=[]; nDrude=0 - name2 = [] - for layer in film: - try: - a = layer[0]['name'] - except: - a = 'None' - self.name.append('{:<10}'.format(a[:10])) # film name and material - name2.append(a) - try: - a = layer[1] - except: - a = 10000. - self.thick.append(a) - self.epsinf.append(layer[0]['eps']) - nTO = 2 * len(layer[0]['wTO']) - allTO.extend(layer[0]['wTO']) - allgTO.extend(layer[0]['gTO']) - allTO.extend(layer[0]['wLO']) - allgTO.extend(layer[0]['gLO']) - Q.extend(2*len(layer[0]['wTO'])*[10.]) - self.listNOsci.append(nTO) - - if len(allTO)!=sum(self.listNOsci) or len(allgTO)!=sum(self.listNOsci): - print('Error in materials: ', layer[0]) - self.wOsc = np.array(allTO) - self.gOsc = np.array(allgTO) - self.osc = np.array([self.wOsc, np.array(Q), self.gOsc]) - return - - def calcSurfaceLoss(self,x): - ''' Calculate the surface loss spectrum for the array of x, which needs to be an equidistant array. - All parameters are defined in the class __init__() call.''' - wmin = min(x) - wmax = max(x)-0.001 - dw = (wmax-wmin)/(len(x)-1) # assumes that x is an equidistant array - wn_array_size = len(x) # size of array for x and epsilon (wn_array, loss_array) - nper = 1. - contrl = '{:<10}'.format('None'[:10]) # Can be 'image' to include image charge - mode = '{:<10}'.format('kurosawa'[:10]) - wn_array,loss_array = LambinEELS.mod_doeels.doeels(self.e0,self.theta,self.phia,self.phib, - wmin,wmax,dw,self.layers,self.neps,nper,self.name, - self.thick,self.epsinf,self.listNOsci,self.osc,contrl,mode,wn_array_size) - i=0 - for item in wn_array: - if item > 0: break - i += 1 - return wn_array[i-1:], loss_array[i-1:] - - def calcSurfaceLoss3(self,x): - ''' Calculate the surface loss spectrum for the array of x, which needs to be an equidistant array. - All parameters are defined in the class __init__() call.''' - wmin = min(x) - wmax = max(x)-0.001 - dw = (wmax-wmin)/(len(x)-1) # assumes that x is an equidistant array - wn_array_size = len(x) # size of array for x and epsilon (wn_array, loss_array) - nper = 1. - contrl = '{:<10}'.format('None'[:10]) # Can be 'image' to include image charge - mode = '{:<10}'.format('kurosawa'[:10]) - wn_array,loss_array = LambinEELS3.mod_doeels.doeels(self.e0,self.theta,self.phia,self.phib, - wmin,wmax,dw,self.layers,self.neps,nper,self.name, - self.thick,self.epsinf,self.listNOsci,self.osc,contrl,mode,wn_array_size) - i=0 - for item in wn_array: - if item > 0: break - i += 1 - return wn_array[i-1:], loss_array[i-1:] - - def calcHREELS(self,x, normalized=True): - emin = min(x) - emax = max(x)-0.001 - norm = 1 - xLoss,loss_array = self.calcSurfaceLoss(x) - wmin = min(xLoss) - wmax = max(xLoss) - xOut,spectrum,n = LambinBoson.doboson3(self.temperature,self.width,self.gauss,self.asym, - emin,emax,wmin,wmax,loss_array,self.debug,len(loss_array)) - if normalized: - norm = max(spectrum[:n]) - # return xOut[:n], spectrum[:n]/norm - return xOut[:len(x)], spectrum[:len(x)]/norm - - def calcEps(self, x): - epsArray = [] - nOsci = len(self.wOsc) - for wn in x: - yn = LambinEELS.mod_doeels.seteps(self.listNOsci,nOsci,self.osc,self.epsinf,wn,self.layers) - epsArray.append(yn) - return np.transpose(np.array(epsArray)) - - def calcEps3(self, x): - epsArray = [] - nOsci = len(self.wOsc) - for wn in x: - yn = LambinEELS3.mod_doeels.seteps(self.listNOsci,nOsci,self.osc,self.epsinf,wn,self.layers) - epsArray.append(yn) - return np.transpose(np.array(epsArray)) - -#################################################################################### -def myMain(): - import matplotlib.pyplot as plt - oxide = importMaterials('NiO') - substrate = importMaterials('Ag') - x = np.linspace(-100.,1000.,1000) - for d in [4., 20., 50., 200., 10000.]: - film1 = lambin(film=[[oxide,d],[substrate,1000.]]) - # xs, spectrum = film1.calcSurfaceLoss(x) - xs, spectrum = film1.calcHREELS(x, normalized=True) - plt.plot(xs,spectrum, label='{:.0f}'.format(d)) - plt.ylabel('HREELS Intensity') - plt.xlabel('Energy Loss (cm$^{-1}$)') - plt.legend(title='Thickness in ($\AA$)') - plt.show() - -if __name__ == '__main__': - myMain() \ No newline at end of file diff --git a/libhreels/calcHREELS3.py b/libhreels/calcHREELS20.py old mode 100644 new mode 100755 similarity index 63% rename from libhreels/calcHREELS3.py rename to libhreels/calcHREELS20.py index 785ecd6f7713f0bbbcb1ebd9344726cb0b59d2ab..2c0715aa067994b4ae05973c786e334a9689f8cd --- a/libhreels/calcHREELS3.py +++ b/libhreels/calcHREELS20.py @@ -5,17 +5,19 @@ import sys, re, os from libhreels.HREELS import myPath from copy import deepcopy +import scipy.integrate as integrate + libDir = os.path.dirname(os.path.realpath(__file__)) if not (sys.version.startswith('3.6') or sys.version.startswith('3.8')): - print('''Make sure the Fortran routines 'myEels2' and 'myBoson' + print('''Make sure the Fortran routines 'myEels20' and 'myBoson' have been complied with the proper f2py3 for the right python version!!''') try: - from . import myEels3 as LambinEELS # wrapper for myEels3.f90 + from libhreels import myEels20 as LambinEELS # wrapper for myEels20.f90 from libhreels import myBoson as LambinBoson # wrapper for myBoson.f90 except: - print('myEels2 and MyBoson are not available here (Check your version)') + print('myEels20 and MyBoson are not available here (Check your version)') # Experimental setup as dictionary: setup = { @@ -45,31 +47,41 @@ def importMaterials(string='', path=libDir): materials = json.load(json_file) try: mat = materials[string] + + #Following if-case for preventing old material data to be used in wrong terminology + if mat["wLO"][0]==1 and mat["gLO"][0]==1 and mat["wTO"][0]<0 and mat["gTO"][0]<0: + print('It seems you are trying to load a Material from an older version. Parameter will be altered to fit the current Version.') + wTO = 0 + gTO = -mat["gTO"][0] + wLO = -mat["wTO"][0] + gLO = -mat["gTO"][0] + mat["wTO"][0] = wTO + mat["gTO"][0] = gTO + mat["wLO"][0] = wLO + mat["gLO"][0] = gLO except: print('No data for material >>{}<< found in {} materials.json!!'.format(string, path)) print('Available materials:\n{}\n'.format(materials.keys())) mat = 'None' return mat - -def addDrude(wPlasma, gPlasma, material): - ''' Adds a simple Drude response to the materials properties (which are provided - as 3rd argument) and returns a new materials dictionary with all phonon parameters. Note +def addDrude(wLOPlasma, gLOPlasma, material, gTOPlasma='None'): + ''' Adds a generalized Drude response to the materials properties (which are provided + as last argument) and returns a new materials dictionary with all phonon parameters. Note that at least the eps_infinity has to given before. ''' + if gTOPlasma == 'none': + gTOPlasma = gLOPlasma newMaterial = deepcopy(material) # To select the Drude response, the Lambin Fortran code requires a negative first oscillator frequency # which is then interpreted as plasma frequency. # The values of the former LO parameter are irrelevant (but need to be provided). - if not wPlasma > 0: - print('Cannot add Drude to material',material) - return material try: if len(newMaterial['wTO']) > 0: - newMaterial['wTO'] += [-1*wPlasma] - newMaterial['gTO'] += [-1*gPlasma] - newMaterial['wLO'] += [0] - newMaterial['gLO'] += [0] + newMaterial['wTO'] += [0.] + newMaterial['gTO'] += [gTOPlasma] + newMaterial['wLO'] += [wLOPlasma] + newMaterial['gLO'] += [gLOPlasma] return newMaterial except: print('Cannot add Drude to material',material) @@ -91,7 +103,7 @@ class lambin: self.asym = instrument['asym'] self.layers = len(film) # number of layers self.neps = self.layers - name_size = self.layers + # name_size = self.layers self.name = []; self.thick=[]; self.listNOsci=[]; self.epsinf =[]; Q = [] allTO=[]; allgTO=[]; allgLO=[]; nDrude=0 name2 = [] @@ -142,7 +154,7 @@ class lambin: i += 1 return wn_array[i-1:], loss_array[i-1:] - def calcHREELS(self,x, normalized=True): + def calcHREELS(self,x, normalized=True, areanormalized=False): emin = min(x) emax = max(x)-0.001 norm = 1 @@ -153,7 +165,22 @@ class lambin: emin,emax,wmin,wmax,loss_array,self.debug,len(loss_array)) if normalized: norm = max(spectrum[:n]) - # return xOut[:n], spectrum[:n]/norm + if areanormalized: #edit by HHE + try: + areanormalize_xstart = np.argmin(abs(x+100.)) #seems to be oddly complicated, but is way more stable than x.index(-100.) or where() + except: + areanormalize_xstart = 0 + try: + areanormalize_xend = np.argmin(abs(x-1000.)) + except: + areanormalize_xend = len(x) + cropped_spectra=spectrum[areanormalize_xstart:areanormalize_xend] + cropped_x=x[areanormalize_xstart:areanormalize_xend] + + norm=integrate.simps(cropped_spectra, dx=x[areanormalize_xstart+1]-x[areanormalize_xstart]) + + else: + print("not normalized") return xOut[:len(x)], spectrum[:len(x)]/norm def calcEps(self, x): @@ -171,16 +198,40 @@ def myMain(): from copy import deepcopy import os - x = np.linspace(-100.,900,1100) + x = np.linspace(-100.,1000,400) + x2 = np.linspace(-100.,1000,2400) + material = {'eps': 4., 'wTO': [0], 'gTO': [20], 'wLO': [598.7], 'gLO': [20], } + + print("Ag imported as: ", importMaterials('Ag')) film1 = lambin(film=[[material,10000.]]) film1.temperature = 100 - xs, spectrum = film1.calcSurfaceLoss(x) - plt.plot(xs,spectrum) + #xs, spectrum = film1.calcSurfaceLoss(x) + xs, spectrum = film1.calcHREELS(x,normalized=False,areanormalized=False) + xs1, spectrum1 = film1.calcHREELS(x2,normalized=False,areanormalized=False) + xs2, spectrum2 = film1.calcHREELS(x,areanormalized=True) + xs3, spectrum3 = film1.calcHREELS(x2,areanormalized=True) + + norm_test=integrate.simps(spectrum, dx=xs[2]-xs[1])/(max(xs)-min(xs)) + norm_test1=integrate.simps(spectrum1, dx=xs1[2]-xs1[1])/(max(xs1)-min(xs1)) + norm_test2=integrate.simps(spectrum2, dx=xs2[2]-xs2[1])/(max(xs2)-min(xs2)) + norm_test3=integrate.simps(spectrum3, dx=xs3[2]-xs3[1])/(max(xs3)-min(xs3)) + print(norm_test,norm_test1,norm_test2,norm_test3) + print(xs[2]-xs[1],xs1[2]-xs1[1],xs2[2]-xs2[1],xs3[2]-xs3[1]) + + plt.plot(xs[:-1],spectrum[:-1], label='normalized=Flase, '+str(len(x))+' points, E0 at '+str(max(spectrum))) + plt.plot(xs1[:-2],spectrum1[:-2], label='normalized=Flase, '+str(len(x2))+' points, E0 at '+str(max(spectrum1))) + plt.plot(xs2[:-10],spectrum2[:-10],label='area normalized=True, '+str(len(x))+' points, E0 at '+str(max(spectrum2))) + plt.plot(xs3[:-20],spectrum3[:-20],label='area normalized=True, '+str(len(x2))+' points, E0 at '+str(max(spectrum3))) + + print('spec2/spec0=', max(spectrum2)/max(spectrum)) + print('spec3/spec0=', max(spectrum3)/max(spectrum)) + print('spec0/spec2=', max(spectrum)/max(spectrum2)) + print('spec0/spec3=', max(spectrum)/max(spectrum3)) # pureDrude = {'eps': 4., # 'wTO': [-400], 'gTO': [-20], 'wLO': [0], 'gLO': [0], diff --git a/libhreels/dielectrics.py b/libhreels/dielectrics.py deleted file mode 100644 index 5a98d6798f5fcbf630e2f2ff7f77963a98afbc1e..0000000000000000000000000000000000000000 --- a/libhreels/dielectrics.py +++ /dev/null @@ -1,140 +0,0 @@ -# Testing dielectric models and their surface loss function -import matplotlib.pyplot as plt -import numpy as np -from numpy import real, imag, sqrt -#from scipy.constants import physical_constants - -# Thz = 100*physical_constants['speed of light in vacuum'][0] # conversion from cm^-1 to Hz - -def doping2plasmaFrequency(doping,epsInfinity=1.): - '''Returns the bulk plasma frequency in cm-1 for a doping - given by the argument in cm-3. - Check if eps_Infinity handling is correct.''' - return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(epsInfinity) - -def doping2surfacePlasma(doping,epsInfinity=1.): - '''Returns the surface plasma frequency in cm-1 for a doping - given by the argument in cm-3''' - return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(1+epsInfinity) - -def loss(eps): - return np.imag(1/eps) - -def surfaceLoss(eps): - return np.imag(1/(1+eps)) - -def reflectivity(eps): # IR reflectivity - sq = np.sqrt(eps) - a =np.abs((sq-1)/(sq+1)) - return a*a - -def sigma(eps,w,eps_Infinity=1): # Complex conductivity - return np.imag((eps-eps_Infinity)*w) - -def plotDielectrics(x,eps): - fig, axs = plt.subplots(3, 1, sharex=True, figsize=(6,6)) - axs[0].plot(x, -np.imag(eps), label='-Im $\epsilon (\omega )$') - axs[0].plot(x, np.real(eps), label='Re $\epsilon (\omega )$') - axs[0].legend() - axs[0].set_ylabel('Dielectric Function') - - axs[1].plot(x, surfaceLoss(eps), linestyle='-') - axs[1].set_ylabel('Surface Loss Function') - axs[1].set_ylim([0,None]) - - axs[2].plot(x, reflectivity(eps)) - axs[2].set_ylabel('Reflectivity') - axs[2].set_xlabel('Frequency') - axs[2].set_ylim([0,1]) - axs[0].set_xlim(left=0, right=max(x)) - - plt.show() -class oscillator: - def __init__(self, wTO, wLO, gammaTO=20, gammaLO=20): - self.wTO = wTO - self.wLO = wLO - self.gammaTO = gammaTO - self.gammaLO = gammaLO - - def eps(self,w): - nom = self.wLO*self.wLO - w*w + 1j*self.gammaLO*w - denom = self.wTO*self.wTO - w*w + 1j*self.gammaTO*w - return nom/denom - - def __call__(self, w): - return self.eps(w) - -class simpleOsci(oscillator): - def __init__(self, wTO, Q, gammaTO=20): - self.wTO = wTO - self.Q = Q - self.gammaTO = gammaTO - - def eps(self,w): - nom = self.wTO*self.wTO*self.Q - denom = self.wTO*self.wTO - w*w + 1j*self.gammaTO*w - return nom/denom - -class drude: - '''Returns the additive Drude response with plasma frequency and damping as parameters. Note that for the - full dielctric response, this Drude contribution has to be added to eps_infinity and any phonon - contribution. - ''' - def __init__(self, wPL, gamma): - self.wPL = wPL - self.gamma = gamma - - def eps(self,w): - # Starting from libhreels version 1.1.4, it has been corrected: - return -(self.wPL*self.wPL)/(w*w - 1j*self.gamma*w) - - - def __call__(self, w): - return self.eps(w) - -class drude2: - '''Returns the additive extended Drude response with three arguments: plasma frequency and two damping - parameters. The first damping parameter is the plasma damping and the second the damping at frequency - zero. Note that for the full dielectric response, this Drude contribution has to be added to eps_infinity - and any phonon contribution. - - eps = drude2(omega_p, gamma_p, gamma_0) - ''' - def __init__(self, omega_p, gamma_p, gamma_0=None): - self.omega_p = omega_p - self.gamma_p = gamma_p - if gamma_0: - self.gamma_0 = gamma_0 - else: - self.gamma_0 = gamma_p - - def eps(self,w): - newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero - w = newW - return -((self.omega_p**2 + 1j*(self.gamma_p-self.gamma_0)*w)/(w*(w - 1j *self.gamma_0))) - - def __call__(self, w): - return self.eps(w) - -def myMain(): - # BaTiO3 - data = [ - [178.6, 8.02, 3.09800E-02], - [270.6, 24.00, 10.0800E-02], - [522.9, 1.20, 6.97200E-02]] - # Width (3rd parameter) is given in fraction of TO frequency - - oscis = [simpleOsci(TO, Q, f*TO) for (TO, Q, f) in data] - - x = np.linspace(5,1000,num=1200) - epsInfinity = 5.25 - eps = epsInfinity - for each in oscis: - eps += each(x) - - # Now add Drude contribution: - eps += drude2(1200,60)(x) - plotDielectrics(x, eps) - -if __name__ == '__main__': - myMain() diff --git a/libhreels/dielectrics20.py b/libhreels/dielectrics20.py new file mode 100755 index 0000000000000000000000000000000000000000..39b8ff1527386f4a275faba7b0b3605edaa0007f --- /dev/null +++ b/libhreels/dielectrics20.py @@ -0,0 +1,242 @@ +# Testing dielectric models and their surface loss function +import matplotlib.pyplot as plt +import numpy as np +from numpy import real, imag, sqrt +from scipy import constants + + +def doping2chargecarrierdensity(doping): #use doping in percentage by mass, not volume + '''Returns the charge carrier density for a given doping in percentage of mass.''' + return float(doping*3.31664067935517E+020*1E6) # doping[%] * charge carrier/cm^3 * 1E6 = n_e [m^-3] + +def doping2plasmaFrequency(doping,epsInfinity=1.): #use doping in percentage by mass, not volume + '''Returns the plasma frequency for a given doping in percentage of mass.''' + eps_0 = constants.value("vacuum electric permittivity") + e = constants.value('elementary charge') + m_e = constants.value('electron mass') + c = constants.value('speed of light in vacuum') + eps_r = epsInfinity + n_e = doping2chargecarrierdensity(doping) + + w_P = np.sqrt(n_e*e*e/(eps_r*eps_0*m_e)) #Hz + w_P = w_P /1E+12 * 33.35641 #Hz -> THz -> cm^-1 + return w_P + + +def doping2surfacePlasma(doping,epsInfinity=1.): + '''Returns the surface plasma frequency in cm-1 for a doping + given by the argument in cm-3''' + return np.sqrt(doping)*1.8817885780819758e-06/np.sqrt(1+epsInfinity) + +def loss(eps): + '''Returns the loss function for a given eps''' + return np.imag(eps) #HHE changed former sign + +def surfaceLoss(eps): + '''Returns the surface loss function for a given eps''' + return np.imag(-1/(1+eps)) #HHE changed former sign, according to Phd of FSc + +def reflectivity(eps): # IR reflectivity + '''Returns the reflectivity for a given eps''' + sq = np.sqrt(eps) + a = np.abs((sq-1)/(sq+1)) + return a*a + +def sigma(eps,w): + '''Returns the complex conductivity for a given eps''' + return (eps-1)*w/1j + +def plotDielectrics(x,eps, title=" ", plot_show=True): + '''This method plots a given (x,eps)-dataset as Re/Im(eps), SurfaceLoss and Reflectivity + ''' + fig, axs = plt.subplots(3, 1, sharex=True) + fig.suptitle(str(title), fontsize=16) + axs[0].plot(x, np.imag(eps), label='Im( $\epsilon (\omega )$ )') + axs[0].plot(x, np.real(eps), label='Re( $\epsilon (\omega )$ )') + axs[0].legend() + axs[0].set_ylabel('Dielectric Function') + + axs[1].plot(x, surfaceLoss(eps), lineStyle='-', label='Loss function') + axs[1].set_ylabel('Surface Loss Function') + axs[1].set_ylim([-0.10,None]) + + axs[2].plot(x, reflectivity(eps)) + axs[2].set_ylabel('Reflectivity') + axs[2].set_xlabel('Frequency') + axs[2].set_ylim([0.0,1]) + axs[0].set_xlim(left=0, right=max(x)) + + if plot_show==True: + plt.show() + +class oscillator: + '''This defines epsilon of an oscillator with a given (wTO, gTO, wLO, gLO) + ''' + def __init__(self, wTO=177., gTO=20., wLO=184., gLO=20.): + self.wTO = wTO + self.wLO = wLO + self.gammaTO = gTO + self.gammaLO = gLO + + def eps(self,w): + nom = self.wLO*self.wLO - w*w - 1j*w*self.gammaLO #sign changed according to gervais paper, formular 2 -> neg im(eps), therefore changed back to lambins kurosawa + denom = self.wTO*self.wTO - w*w - 1j*w*self.gammaTO #sign changed according to gervais paper, formular 2 -> neg im(eps), therefore changed back to lambins kurosawa + return nom/denom + + def __call__(self, w): + return self.eps(w) + +class drude: + '''This defines an additive Drude component to a dielectric function: + eps_drude = drude2(omega_p, gamma_p, gamma_0) + ... eps(omega) + eps_drude(omega) + ''' + def __init__(self, omega_p, gamma_p, gamma_0): + + self.omega_p = omega_p + self.gamma_p = gamma_p + self.gamma_0 = gamma_0 + + print("called drude with w_P, g_P, g_0:", omega_p, gamma_p, gamma_0) + + + def eps(self,w): + newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero + w = newW + #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0))) + + nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) + denom = w*w + 1j*w*self.gamma_0 #sign changed according to Lambin notation + + return -(nom/denom) #HHE changed + + def __call__(self, w): + return self.eps(w) + + +class cole_davidson: + '''This defines an additive Cole-Davidson component to a dielectric function: + eps_cole_davidson = cole_davidson(omega_p, gamma_p, gamma_0) + ... eps(omega) + eps_drude(omega) + ''' + def __init__(self, omega_p, gamma_p, gamma_0, beta): + + self.omega_p = omega_p + self.gamma_p = gamma_p + self.gamma_0 = gamma_0 + self.beta = beta + + print("called coledavidson with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta) + + + def eps(self,w): + newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero + w = newW + #print("Drude (x,y): ", w, -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0))) + #return -((self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p))/(w*w + 1j*w*self.gamma_0)) #HHE changed + + #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P) + #denom = (w*w - 1j*w*self.gamma_0)**self.beta #sign changed according to gervais paper, formular 7 + + hilf = 1j * w * self.gamma_0 * (1-1j*w*(1/self.gamma_0))**(self.beta) + + nom = self.omega_p**2 - 1j*w*(self.gamma_p) - w*w + hilf #unfortunately also the nominator changes due to cole-davidson + denom = hilf #sign changed according to gervais paper, formular 7 -> did not match, changed back to lambin formularism + + + return -(nom/denom) + + def __call__(self, w): + return self.eps(w) + +class cole_cole: + '''This defines an additive Cole-Davidson component to a dielectric function: + eps_cole_cole = cole_cole(omega_p, gamma_p, gamma_0, beta) + ... + ''' + def __init__(self, omega_p, gamma_p, gamma_0, beta): + + self.omega_p = omega_p + self.gamma_p = gamma_p + self.gamma_0 = gamma_0 + self.beta = beta + + print("called colecole with w_P, g_P, g_0, beta:", omega_p, gamma_p, gamma_0, beta) + + def eps(self,w): + newW = np.where(w==0, 0.000567894, w) # Avoid zeros to avoid devision by zero + w = newW + #nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) #very simple approach by just putting an exponent to denominator (works only for gamma_0=gamma_P) + #denom = w*w + w*(1j*self.gamma_0)**self.beta + + hilf = 1j**(self.beta+1) * w**(self.beta+1) * self.gamma_0**(1-self.beta) + + nom = self.omega_p**2 + 1j*w*(self.gamma_0-self.gamma_p) - w*w - hilf #unfortunately also the nominator changes due to cole-cole + denom = 1j*w*(self.gamma_0) - hilf #sign changed according to lambin formularism + + return -(nom/denom) + + def __call__(self, w): + return self.eps(w) + + +class simpleOscillator(oscillator): + '''This defines epsilon of an oscillator with a given (wPl, Q and gTO) + ''' + def __init__(self, wPl, Q, gTO=20): + self.wPl = wPl + self.Q = Q + self.gammaTO = gTO + + def eps(self,w): + nom = self.wPl*self.wPl*self.Q + denom = self.wPl*self.wPl - w*w - 1j*self.gammaTO*w #HHE changed, wrong sign of Im + return nom/denom + +class simpleDrude: + '''This defines epsilon of an oscillator with a given (wPl, gamma). gamma = gammaTO = gammaLO + ''' + def __init__(self, wPL, gamma): + self.wPL = wPL + self.gamma = gamma + + def eps(self,w): + return -(self.wPL*self.wPL)/(w*w-1j*w*self.gamma) #HHE changed + + def __call__(self, w): + return self.eps(w) + + + +def myMain(): + # BaTiO3 + simple_data = [ + [178.6, 8.02, 3.09800E-02], + [270.6, 24.00, 10.0800E-02], + [522.9, 1.20, 6.97200E-02]] + + data = [ + [177.4 , 1.9262 , 184.003 , 9.718], + [272.84 , 93.3091 , 470.972 , 14.31 ], + [506.379 , 42.60 , 739.651 , 33.044], + ] + + x = np.linspace(0,1000,num=1200) + epsInfinity = 5.25 + + simpleoscis = [simpleOscillator(TO, Q, f*TO) for (TO, Q, f) in simple_data] + eps = epsInfinity + for each in simpleoscis: + eps += each(x) #sum + plotDielectrics(x, eps, title="Eps of BaTiO calced with simpleOscillator", plot_show=False) + + oscis = [oscillator(wTO=a, wLO=c, gTO=b, gLO=d) for (a, b, c, d) in data] + eps2 = epsInfinity + for each in oscis: + eps2 *= each(x) #kurosawa + plotDielectrics(x, eps2, title="Eps of BaTiO calced with oscillator", plot_show=False) + + plt.show() + +if __name__ == '__main__': + myMain() diff --git a/libhreels/materials.json b/libhreels/materials.json old mode 100644 new mode 100755 index 02518d7fed00d079c21a93c412d12297a6dcecab..bc95115f29de261cc59ee045862b0d8d2cba7320 --- a/libhreels/materials.json +++ b/libhreels/materials.json @@ -79,7 +79,7 @@ "gLO": [ 3.09999, 4.62163, 25.0 ], "reference": "Gervais'93", "comment": "300 K", - "name": "SrTiO_3" + "name": "Gervais_SrTiO_3" }, "Gervais_doped": { "eps": 5.01, @@ -91,7 +91,7 @@ "comment": "300 K", "name": "SrTiO_3" }, - "J0B04": { + "J0B04": { "eps": 5.2, "wTO": [60, 171.8, 531.97], "gTO": [80, 30, 17.706], @@ -100,6 +100,35 @@ "reference": "Fit to J0B04", "name": "SrTiO_3" }, + "Fitted_EBE02" : { + "eps": 5.2, + "wTO": [94.85, 178.28, 528.921], + "gTO": [2.00, 9.14, 18.73], + "wLO": [176.49, 476.26, 833.7236], + "gLO": [10.43, 18.26 , 52.63], + "reference": "Fit by HHE", + "name": "Fitted_EBE02_300K" + }, + "Fitted_EBE02_210914" : { + "eps": 3.024434, + "wTO": [147.877698, 175.0514, 531.3065], + "gTO": [5.02880219, 13.73305, 18.76442], + "wLO": [168.951344, 475.621265, 833.761677], + "gLO": [10.0004693, 19.75355, 60.70802], + "reference": "Fit by HHE, 0.05% doped, 300K, fitted with calcHREELS3, last updated: 210916", + "name": "Fitted_EBE02_300K" + }, + "Fitted_K7701_210916" : { + "eps": 4.57655596, + "wTO": [50.0211, 174.060581, 534.354046], + "gTO": [1.56852054, 8.17815492, 12.6232742], + "wLO": [171.56, 482.400957, 820.342514], + "gLO": [8.52832484, 21.0070097, 58.6519913], + "reference": "Fit by HHE, 0.05% doped, 114K, fitted with calcHREELS3, last updated: 210917", + "name": "Fitted_K7701_100K" + }, + + "STO_DEG1": { "eps": 5.14285, "wTO": [ 98.0 , 174.367 , 548.821, 0.0 ], @@ -109,15 +138,6 @@ "reference": "not working", "name": "SrTiO_3" }, - "BTO": { - "eps": 5.1295, - "wTO": [177.4 , 272.84 , 506.379], - "gTO": [ 1.9262 , 93.3091 , 42.60 ], - "wLO": [184.003 , 470.972 , 739.651], - "gLO": [ 9.718 , 14.31 , 33.044], - "reference": "BTO323", - "name": "BaTiO_3" - }, "KTO": { "eps": 5.12878, "wTO": [91.0, 201.0, 548.0, 710.0, 752.0], @@ -172,6 +192,15 @@ "reference": "Irani'71, Ehrenreich'62, Furtak'75", "name": "Silver" }, + "Ag4": { + "eps": 1.0, + "wTO": [0], + "gTO": [161.0], + "wLO": [31460], + "gLO": [161.0], + "reference": "Irani'71, Ehrenreich'62, Furtak'75", + "name": "Silver new" + }, "Pt": { "eps": 1.0, "wTO": [-185541.0], @@ -181,6 +210,16 @@ "reference": "Weaver'75, Seignac'72", "name": "Platinum" }, + "Pt4": { + "eps": 1.0, + "wTO": [0], + "gTO": [16134.0], + "wLO": [185541.0], + "gLO": [16134.0], + "reference": "Weaver'75, Seignac'72", + "name": "Platinum new" + }, + "test": { "eps": 5.25, "wTO": [493.7, -10.0], diff --git a/libhreels/materials20.json b/libhreels/materials20.json new file mode 100755 index 0000000000000000000000000000000000000000..e9dee5f28379396c96c37cbe7c2e5e8bf95a207a --- /dev/null +++ b/libhreels/materials20.json @@ -0,0 +1,224 @@ +{ + "NiO": { + "eps": 5.25, + "wTO": [393.7], + "gTO": [10.8], + "wLO": [584.7], + "gLO": [10.8], + "reference": "Schumann'22", + "name": "NiO" + }, + "BaO": { + "eps": 3.7, + "wTO": [145.0], + "gTO": [46.0], + "wLO": [422.0], + "gLO": [56.0], + "reference": "Goian'18", + "name": "BaO" + }, + "BaO_2": { + "eps": 3.7, + "wTO": [145.0], + "gTO": [46.0], + "wLO": [414.0], + "gLO": [56.0], + "reference": "...", + "name": "BaO" + }, + "Willet": { + "eps": 4.12, + "wTO": [188.0, 427.0, 495.72, 650.79, 708.2], + "gTO": [0.4, 5.0, 3.8, 22.5, 55.3], + "wLO": [276.4, 596.1, 495.5, 744.1, 702.2], + "gLO": [3.7, 7.2, 3.8, 12.1, 66.0], + "reference": "Willet-Gies'14", + "name": "LaAlO_3" + }, + "STO": { + "eps": 5.14285, + "wTO": [ 98.0 , 174.367 , 548.821 ], + "gTO": [ 62.488 , 9.89216, 25.0 ], + "wLO": [171.679 , 475.624 , 817.369 ], + "gLO": [ 6.11124, 29.9772 , 60.0 ], + "reference": "Crandles'99", + "name": "SrTiO_3" + }, + "Crandles99_2e18_80K": { + "eps": 5.16, + "wTO": [ 50.5 , 174.2 , 548.9 ], + "gTO": [ 16 , 6.0 , 10.0 ], + "wLO": [171.56 , 481.0 , 805.7 ], + "gLO": [ 4.8 , 5.35 , 29.0 ], + "reference": "Crandles'99 2e18 80K", + "name": "SrTiO_3" + }, + "Crandles99_2e18_300K": { + "eps": 5.14, + "wTO": [93.9831, 173.175, 543.346 ], + "gTO": [25.8272, 8.42197, 17.5385 ], + "wLO": [169.788, 474.634, 795.385 ], + "gLO": [ 4.9, 7.07977, 28.3571 ], + "reference": "Crandles'99 2e18 300K", + "name": "SrTiO_3" + }, + "Crandles99_2e19_80K": { + "eps": 5.1, + "wTO": [50.0211, 173.071, 548.723], + "gTO": [34.8738, 5.00001, 10.0], + "wLO": [ 170.94, 493.881, 823.411], + "gLO": [4.69999, 15.3783, 52.6903], + "reference": "Crandles'99 2e19 80K", + "name": "SrTiO_3" + }, + "Gervais": { + "eps": 5.2, + "wTO": [ 88.5533 , 175.999 , 543.0 ], + "gTO": [ 27.9999 , 5.99999, 18.7496 ], + "wLO": [171.967 , 475.999 , 796.556 ], + "gLO": [ 3.09999, 4.62163, 25.0 ], + "reference": "Gervais'93", + "comment": "300 K", + "name": "Gervais_SrTiO_3" + }, + "Gervais_doped": { + "eps": 5.01, + "wTO": [ 93.0 , 175.99 , 547.99 ], + "gTO": [ 42.874 , 8.9935, 18.999 ], + "wLO": [172.79 , 472.53 , 795.0 ], + "gLO": [ 4.1918, 9.9999, 26.0 ], + "reference": "Gervais'93", + "comment": "300 K", + "name": "SrTiO_3" + }, + "J0B04": { + "eps": 5.2, + "wTO": [60, 171.8, 531.97], + "gTO": [80, 30, 17.706], + "wLO": [170.9, 481, 814.92], + "gLO": [53.6, 13.3 , 36.497], + "reference": "Fit to J0B04", + "name": "SrTiO_3" + }, + "Fitted_EBE02" : { + "eps": 5.2, + "wTO": [94.85, 178.28, 528.921], + "gTO": [2.00, 9.14, 18.73], + "wLO": [176.49, 476.26, 833.7236], + "gLO": [10.43, 18.26 , 52.63], + "reference": "Fit by HHE", + "name": "Fitted_EBE02_300K" + }, + "Fitted_EBE02_210914" : { + "eps": 3.024434, + "wTO": [147.877698, 175.0514, 531.3065], + "gTO": [5.02880219, 13.73305, 18.76442], + "wLO": [168.951344, 475.621265, 833.761677], + "gLO": [10.0004693, 19.75355, 60.70802], + "reference": "Fit by HHE, 0.05% doped, 300K, last updated: 210916", + "name": "Fitted_EBE02_300K" + }, + "Fitted_K7701_210916" : { + "eps": 4.57655596, + "wTO": [50.0211, 174.060581, 534.354046], + "gTO": [1.56852054, 8.17815492, 12.6232742], + "wLO": [171.56, 482.400957, 820.342514], + "gLO": [8.52832484, 21.0070097, 58.6519913], + "reference": "Fit by HHE, 0.05% doped, 114K, last updated: 210917", + "name": "Fitted_K7701_100K" + }, + + + "STO_DEG1": { + "eps": 5.14285, + "wTO": [ 98.0 , 174.367 , 548.821, 0.0 ], + "gTO": [ 62.488 , 9.89216, 25.0 , -100 ], + "wLO": [171.679 , 475.624 , 817.369, 200 ], + "gLO": [ 6.11124, 29.9772 , 60.0 , 0.0 ], + "reference": "not working", + "name": "SrTiO_3" + }, + "BTO": { + "eps": 5.1295, + "wTO": [177.4 , 272.84 , 506.379], + "gTO": [ 1.9262 , 93.3091 , 42.60 ], + "wLO": [184.003 , 470.972 , 739.651], + "gLO": [ 9.718 , 14.31 , 33.044], + "reference": "BTO323", + "name": "BaTiO_3" + }, + "KTO": { + "eps": 5.12878, + "wTO": [91.0, 201.0, 548.0, 710.0, 752.0], + "gTO": [6.5, 6.0, 13.0, 130.0, 20.0], + "wLO": [186.0, 423.0,701.0, 750.0, 821.0], + "gLO": [6.0, 6.5, 130.0, 20.0, 13.0], + "reference": "Jandl'91", + "name": "KTaO_3" + }, + "Massa97": { + "eps": 5.12878, + "wTO": [1317.9, 0.0, 181.8, 209.9, 330.5], + "gTO": [1500.5, 1170.0, 185.1, 321.8, 391.5], + "wLO": [1618.5, 900.3, 12.8, 829.7, 49.3], + "gLO": [3967.7, 689.4, 5.0, 41.9, 85.9], + "reference": "Massa'97", + "name": "NdNiO_3" + }, + "Massa97_77K": { + "eps": 5.12878, + "wTO": [ 398.0, 430.4, 470.0, 561.1, 729.9, 0.0], + "gTO": [ 428.8, 435.8, 542.1, 567.8, 1108.5, 146.6], + "wLO": [ 149.2, 31.5, 254.5, 107.3, 1155.6, 1705.3], + "gLO": [ 887.6, 28.0, 290.8, 110.8, 1504.6, 2087.1], + "reference": "Massa'97", + "name": "NdNiO_3" + }, + "S2R1O4": { + "eps": 5.5, + "wTO": [195, 369, 412.5, 466 ], + "gTO": [9.2, 18, 17.1, 18 ], + "wLO": [250, 400, 450, 600], + "gLO": [9.2, 18, 17.1, 18 ], + "reference": "Farbe der Ruthenate", + "name": "SRO-214" + }, + "Ag": { + "eps": 1.0, + "wTO": [0], + "gTO": [161.0], + "wLO": [31460.0], + "gLO": [161.0], + "reference": "Irani'71, Ehrenreich'62, Furtak'75", + "name": "Silver" + }, + "Pt": { + "eps": 1.0, + "wTO": [0], + "gTO": [16134.0], + "wLO": [185541.0], + "gLO": [16134.0], + "reference": "Weaver'75, Seignac'72", + "name": "Platinum" + }, + "test": { + "eps": 5.25, + "wTO": [493.7, -10.0], + "gTO": [10.8, -5.0], + "wLO": [584.7, 1], + "gLO": [10.8, 1], + "reference": "for playing", + "name": "test" + }, + "Vacuum": { + "eps": 1.0, + "wTO": [1000], + "gTO": [100], + "wLO": [1001], + "gLO": [100], + "reference": "...", + "name": "vacuum" + } +} + + diff --git a/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so b/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so new file mode 100755 index 0000000000000000000000000000000000000000..59bf896c1432e147d62c4476225db64664fc9923 Binary files /dev/null and b/libhreels/myEels20.cpython-38-x86_64-linux-gnu.so differ diff --git a/pyproject.toml b/pyproject.toml index 853c05196388c2cb8a949a8c3f4160b562baad7e..bea482e9c73419dc3812a20053ce1225d9561ded 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [tool.poetry] name = "libhreels" -version = "1.4.0" +version = "2.0.0" description = "Handling, simulating, and plotting HREELS and Auger spectroscopy data" authors = ["Wolf Widdra <wolf.widdra@physik.uni-halle.de>"] include = ["*./libhreels/*"] @@ -27,12 +27,12 @@ pytest = "^5.2" [tool.poetry.scripts] dielectrics = "libhreels.dielectrics:myMain" -calchreels = "libhreels.calcHREELS:myMain" +calchreels = "libhreels.calcHREELS20:myMain" viewhreelstest = "libhreels.ViewHREELS:myMain" [tool.poetry.plugins."gui_scripts"] viewauger = "libhreels.ViewAuger:myMain" -viewhreels = "libhreels.ViewHREELS:myMain" +viewhreels = "libhreels.ViewHREELS20:myMain" [build-system] requires = ["poetry>=0.12"]