[ase2-svncheckins] [svn commit /var/www/svn/CamposASE2] r1084 - trunk/ASE/Transport

svn at fysik.dtu.dk svn at fysik.dtu.dk
Tue Oct 28 14:57:32 CET 2008


Author: carstenr
Date: 2008-10-28 14:57:32 +0100 (Tue, 28 Oct 2008)
New Revision: 1084

Modified:
   trunk/ASE/Transport/LOE.py
Log:
Updated the Lowest Order Expansion in the el-ph coupling


Modified: trunk/ASE/Transport/LOE.py
===================================================================
--- trunk/ASE/Transport/LOE.py	2008-09-02 15:49:40 UTC (rev 1083)
+++ trunk/ASE/Transport/LOE.py	2008-10-28 13:57:32 UTC (rev 1084)
@@ -1,5 +1,5 @@
 import Numeric as num
-from ASE.Transport.Hilbert  import Hilbert, HilbertKernelInterpolate
+from ASE.Transport.Hilbert import Hilbert, HilbertKernelInterpolate
 
 def Dagger(matrix):
     return num.conjugate(num.swapaxes(matrix, 0, 1))
@@ -47,9 +47,12 @@
                  phononcoupling,
                  frequencies,
                  energyrange,
-                 fermienergy):
-        print 'Remember to align fermi levels, and'
-        print 'check principal layer convergence'
+                 fermienergy=0.0):
+        """Initialize necessary attributes.
+
+        scatteringgreenfunction is any object with the methods SetEnergy,
+        GetRepresentation, GetLeftLeadSelfEnergy, and GetRightLeadSelfEnergy.
+        """
         self.frequencies = frequencies
         self.phononcoupling = phononcoupling
         self.temperature = temperature
@@ -62,32 +65,29 @@
         Gd = Dagger(G)
         sigmal = gf.GetLeftLeadSelfEnergy()
         sigmar = gf.GetRightLeadSelfEnergy()
+        del gf
+        
         Gam1 = 1.j * (sigmal - Dagger(sigmal))
         Gam2 = 1.j * (sigmar - Dagger(sigmar))
-        del gf, sigmal, sigmar
-        
+        del sigmal, sigmar
+
         GGam2 = Multiply(G, Gam2)
         Gam2Gd = Multiply(Gam2, Gd)
+        self.ballistic = num.trace(Multiply(G, Gam1, Gd, Gam2)).real #1<->2
+        del Gam2
+        
         GdGam1G = Multiply(Gd, Gam1, G)
         A1 = Multiply(G, Gam1, Gd)
         A2 = Multiply(G, Gam2Gd)
-        ballistic = num.trace(Multiply(G, Gam1, Gd, Gam2))
-        del G, Gd, Gam1, Gam2
+        del G, Gd, Gam1
 
-        print 'Ballistic =', ballistic
-        if abs(ballistic.imag) > 1e-7:
-            print 'Too large imag. part in ballistic prefactor'
-
-        self.ballistic = ballistic.real
-        self.I_sym = []
-        self.I_asym = []
+        self.symmetric = []
+        self.asymmetric = []
         self.hilbert = []
         ker = HilbertKernelInterpolate(2 * len(energyrange))
-
         for w, M in zip(self.frequencies, self.phononcoupling):
             box = Box(w, -w, energyrange, temperature)
             hil = Hilbert(box, ker=ker)
-            assert max(num.absolute(hil.imag)) < 1e-7
             self.hilbert.append(hil.real)
 
             MA1M = Multiply(M, A1, M)
@@ -96,35 +96,44 @@
             MA2mA1M = MA2M - MA1M
             del MA1M
 
-            sym = num.trace(Multiply(-2 * MA2M + 1.j * (
-                Multiply(MAM, GGam2) - Multiply(Gam2Gd, MAM)), GdGam1G))
+            sym = num.trace(Multiply(GdGam1G, MA2M + .5j * (
+                Multiply(Gam2Gd, MAM) - Multiply(MAM, GGam2))))
 
-            asym = num.trace(Multiply((Multiply(Gam2Gd, MA2mA1M) +
-                                       Multiply(MA2mA1M, GGam2)), GdGam1G))
+            asym = num.trace(Multiply(GdGam1G, (Multiply(Gam2Gd, MA2mA1M) +
+                                                Multiply(MA2mA1M, GGam2))))
 
-            print 'sym =', sym, 'asym =',asym
-            if abs(sym.imag) > 1e-7:
-                print 'Too large imag. part in symmetric current, for mode', i
-            if abs(asym.imag) > 1e-7:
-                print 'Too large imag. part in asymmetric current, for mode', i
+            self.symmetric.append(sym.real)
+            self.asymmetric.append(asym.real)
 
-            self.I_sym.append(sym.real)
-            self.I_asym.append(asym.real)
-
-    def GetCurrent(self, V):
-        I = self.ballistic * V + self.GetSymmetric(V) + self.GetAsymmetric(V)
+    def Current(self, V):
+        I = self.Ballistic(V) + self.Symmetric(V) + self.Asymmetric(V)
         return I
 
-    def GetSymmetric(self, V):
+    def Ballistic(self, V):
+        return self.ballistic * V
+
+    def Symmetric(self, V):
         T = self.temperature
-        I_sym = 0.0
-        for w, sym in zip(self.frequencies, self.I_sym):
-            I_sym += .5 * (EBose(w + V, T) - EBose(w - V, T)) * sym
-        return I_sym
+        I = 0.0
+        for w, sym in zip(self.frequencies, self.symmetric):
+            I += (EBose(w - V, T) - EBose(w + V, T)) * sym
+        return I
 
-    def GetAsymmetric(self, V):
-        I_asym = 0.0
-        for asym, hil in zip(self.I_asym, self.hilbert):
-            box = Box(V, 0, self.energyrange, self.temperature)
-            I_asym += -.5 * asym * Trapez(self.energyrange, box * hil)
-        return I_asym
+    def Asymmetric(self, V):
+        I = 0.0
+        for asym, hil in zip(self.asymmetric, self.hilbert):
+            box = Box(0, -V, self.energyrange, self.temperature)
+            I += -.5 * asym * Trapez(self.energyrange, box * hil)
+        return I
+
+    def DifferentialConductance(self, bias):
+
+        I = num.array([self.Current(V) for V in bias])
+
+        # Differential conductance -- central difference
+        dG = num.empty(len(I), num.Float)
+        dG[1:-1] = (I[2:] - I[:-2]) / (bias[2:] - bias[:-2])
+        dG[0] = (I[1] - I[0]) / (bias[1] - bias[0])
+        dG[-1] = (I[-1] - I[-2]) / (bias[-1] - bias[-2])
+
+        return dG



More information about the ase2-svncheckins mailing list