V 2.1.1 (Additive Lorentzian corrected)

libhreels/calcHREELS20.py

@@ -79,6 +79,8 @@ def addDrude(wLOPlasma, gLOPlasma, material, gTOPlasma='None'):
             newMaterial['gTO'] += [gTOPlasma]
             newMaterial['wLO'] += [wLOPlasma]
             newMaterial['gLO'] += [gLOPlasma]
+            if newMaterial.get('Q'):
+                newMaterial['Q'] += [0]
             return newMaterial
     except:
         print('Cannot add Drude to material',material)
@@ -102,7 +104,7 @@ class lambin:
         self.neps = self.layers
         # name_size = self.layers
         self.name = []; self.thick=[]; self.listNOsci=[]; self.epsinf =[]; Q = []
-        allTO=[]; allgTO=[];  allgLO=[]; nDrude=0
+        allTO=[]; allgTO=[];  allgLO=[]; nDrude=0; Qdummy = []
         name2 = []
         for layer in film:
             try:
@@ -122,14 +124,22 @@ class lambin:
             allgTO.extend(layer[0]['gTO'])
             allTO.extend(layer[0]['wLO'])
             allgTO.extend(layer[0]['gLO'])
-            Q.extend(2*len(layer[0]['wTO'])*[10.])
+            qList = layer[0].get('Q')
+            if qList:
+                Q.extend(layer[0]['Q'])
+                Q.extend(len(layer[0]['Q'])*[0.])
+            else:
+                Q.extend(2* len(layer[0]['wTO'])*[0.])
             self.listNOsci.append(nTO)
-        
+
         if len(allTO)!=sum(self.listNOsci) or len(allgTO)!=sum(self.listNOsci):
             print('Error in materials: ', layer[0])
+        if len(Q)!=sum(self.listNOsci) :
+            print('Error in materials (Check Q): ', layer[0])
         self.wOsc = np.array(allTO)
         self.gOsc = np.array(allgTO)
         self.osc = np.array([self.wOsc, np.array(Q), self.gOsc])
+        # print('[self.wOsc, np.array(Q), self.gOsc]: \n',self.osc)
         return
 
     def calcSurfaceLoss(self,x):
@@ -193,21 +203,56 @@ def myMain():
     import matplotlib.pyplot as plt
     import numpy as np
     import os
+    from libhreels import dielectrics20 as dielectrics
 
     x = np.linspace(-100.,1000,400)
-    material = {'eps': 4.,
-                'wTO': [0],  'gTO': [20], 'wLO': [598.7], 'gLO': [20]}
+    material = {'eps': 1.,
+                'wTO': [-200, -750],  
+                'gTO': [  12,   8], 
+                'wLO': [   1,   1],     # this parameter is irrelvant if wTO is negativ
+                'gLO': [   1,   1],     # this parameter is irrelvant if wTO is negativ
+                'Q'  : [  10,  15]}
+    material2 = {'eps': 1.,
+                'wTO': [-200, -750],  
+                'gTO': [  12,   8], 
+                'wLO': [   1,   1],     # this parameter is irrelvant if wTO is negativ
+                'gLO': [   1,   1],     # this parameter is irrelvant if wTO is negativ
+                'Q'  : [  0.1,  15]}
+    material3 = {'eps': 1.,
+                'wTO': [-200, -750,   0],  
+                'gTO': [  12,   8,   50], 
+                'wLO': [   1,   1, 2000],     # this parameter is irrelvant if wTO is negativ
+                'gLO': [   1,   1,   50],     # this parameter is irrelvant if wTO is negativ
+                'Q'  : [  10,  15,    0]}
 
     film1 = lambin(film=[[material,10000.]])
-    film1.temperature = 100
-    xs, spectrum = film1.calcHREELS(x,normalized=True,areanormalized=False)
+    film2 = lambin(film=[[material2,10000.]])
+    film3 = lambin(film=[[material3,10000.]])
+    eps3= film3.calcEps(x)[0]
+    eps2= film2.calcEps(x)[0]
+    eps = film1.calcEps(x)[0]
+    plt.plot(x,np.real(eps))
+    plt.plot(x,np.real(eps2))
+    plt.plot(x,np.real(eps3))
 
-    plt.plot(xs[:-1],spectrum[:-1], label='normalized=False')
-    plt.ylabel('Surface Loss')
+    ############# Comparison with dielectrics ################
+    osci1  = dielectrics.simpleOscillator(material3['wTO'][0], 
+                                          material3['Q'  ][0],
+                                      gTO=material3['gTO'][0])
+    osci2  = dielectrics.simpleOscillator(material3['wTO'][1], 
+                                          material3['Q'  ][1],
+                                      gTO=material3['gTO'][1])
+    drude = dielectrics.drude(material3['wLO'][2],material3['gTO'][2],material3['gLO'][2])
+    epsInfinity = material3['eps']
+    eps_dielectrics = epsInfinity * (osci1(x) + osci2(x)  + (1 + drude(x)))
+    plt.plot(x,np.real(eps_dielectrics),linestyle='dotted')
+
+
+
+    plt.ylabel(r'$Re(\epsilon)$')
     plt.xlabel('Energy Loss (cm$^{-1}$)')
-    plt.ylim(bottom=0)
-    plt.title('Bulk plasmon at 600 cm-1')
-    # plt.legend(title=r'Plasma frequency')
+    plt.xlim(left=5)
+    # plt.ylim(-6000,6000)
 
     plt.text(0.99, 0.01,os.path.basename(__file__), fontsize=10, ha='right', va='bottom', transform=plt.gcf().transFigure)
     output_filename = os.path.splitext(__file__)[0] + '.png'
@@ -215,6 +260,43 @@ def myMain():
 
     plt.show()
 
+    plt.plot(x,np.imag(eps))
+    plt.plot(x,np.imag(eps2))
+    plt.plot(x,np.imag(eps3))
+    plt.plot(x,np.imag(eps_dielectrics))
+    plt.xlim(left=0)
+    plt.ylim(-1500,1500)
+    plt.show()
+
+    # plt.plot(x,np.imag(dielectrics.sigma(eps,x)))
+    # plt.plot(x,np.imag(dielectrics.sigma(eps2,x)))
+    # plt.plot(x,np.imag(dielectrics.sigma(eps3,x)))
+    # plt.xlim(left=0)
+    # plt.ylim(-1500,1500)
+    # plt.show()
+
+    xs, spectrum = film3.calcHREELS(x,normalized=True,areanormalized=False)
+    plt.plot(xs[:-1],spectrum[:-1], label='normalized=False')
+    plt.show()
+
+    material4 = {'eps': 1.,
+                'wTO': [ 200, 750,   0],  
+                'gTO': [  12,   8,   50], 
+                'wLO': [ 600, 950, 2000],     
+                'gLO': [  10,  10,   50],     
+                'Q'  : [  10,  15,    0]}
+
+    # material4 = {'eps': 1.,
+    #             'wTO': [ 200, 750,   0],  
+    #             'gTO': [  12,   8,   50], 
+    #             'wLO': [ 300, 950, 2000], 
+    #             'gLO': [  10,  10,   50]}
+
+
+    film4 = lambin(film=[[material4,10000.]])
+    xs, spectrum = film4.calcHREELS(x,normalized=True,areanormalized=False)
+    plt.plot(xs[:-1],spectrum[:-1], label='normalized=False')
+    plt.show()
 
 if __name__ == '__main__':
 	myMain()
\ No newline at end of file

libhreels/myEels20.f90

@@ -982,7 +982,7 @@ subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps)
         deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) )
       
       else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q
-        addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, wn*osc(3,j))
+        addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -1*wn*osc(3,j))  ! Sign of imaginary part changed (WFW)
       
       else                      ! osc(1,j) = 0   -> it is a Drude term
         addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j))

pyproject.toml

 [tool.poetry]
 name = "libhreels"
-version = "2.1.0"
+version = "2.1.1"
 description = "Handling, simulating, and plotting HREELS and Auger spectroscopy data"
 authors = ["Wolf Widdra <wolf.widdra@physik.uni-halle.de>"]
 include = ["*./libhreels/*"]